From 3c1d71f8fea190564c88d6eff3171f9b8b31c6c7 Mon Sep 17 00:00:00 2001
From: "N. Brouard" <brouard@ined.fr>
Date: Thu, 22 Jun 2023 12:50:51 +0000
Subject: [PATCH] Summary: stil on going

---
 src/imachprax.c | 87 ++++++++++++++++++++++++++++++++++++++++++-------
 1 file changed, 76 insertions(+), 11 deletions(-)

diff --git a/src/imachprax.c b/src/imachprax.c
index b20b69f..ae5bade 100644
--- a/src/imachprax.c
+++ b/src/imachprax.c
@@ -1,6 +1,9 @@
 /* $Id$
   $State$
   $Log$
+  Revision 1.3  2023/06/22 11:28:07  brouard
+  *** empty log message ***
+
   Revision 1.2  2023/06/22 11:22:40  brouard
   Summary: with svd but not working yet
 
@@ -2614,6 +2617,41 @@ void linmin(double p[], double xi[], int n, double *fret,double (*func)(double [
 } 
 
 /**** praxis ****/
+
+void transpose_in_place ( int n, double **a )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+   TRANSPOSE_IN_PLACE transposes a square matrix in place.
+  Licensing:
+    This code is distributed under the GNU LGPL license.
+
+    Input, int N, the number of rows and columns of the matrix A.
+
+    Input/output, double A[N*N], the matrix to be transposed.
+*/
+{
+  int i;
+  int j;
+  double t;
+
+  /* for ( j = 0; j < n; j++ ){ */
+  /*   for ( i = 0; i < j; i++ ) { */
+  for ( j = 1; j <= n; j++ ){
+    for ( i = 1; i < j; i++ ) {
+      /* t        = a[i+j*n]; */
+      /* a[i+j*n] = a[j+i*n]; */
+      /* a[j+i*n] = t; */
+      t        = a[i][j];
+      a[i][j] = a[j][i];
+      a[j][i] = t;
+    }
+  }
+  return;
+}
+
 double pythag( double x, double y )
 
 /******************************************************************************/
@@ -2654,9 +2692,16 @@ double pythag( double x, double y )
 }
 
 /* void svdcmp(double **a, int m, int n, double w[], double **v)  */
-void svdcmp(double **a, int m, int n, double w[]) 
+void svdminfit(double **a, int m, int n, double w[]) 
 {
-  /* Golub 1970 Algol60 */
+  /* From numerical recipes */
+  /* Given a matrix a[1..m][1..n], this routine computes its singular value */
+  /*   decomposition, A = U ·W ·V T . The matrix U replaces a on output. */
+  /*   The diagonal matrix of singular values W is out- put as a vector w[1..n]. */
+  /*   The matrix V (not the transpose V T ) is output as v[1..n][1..n]. */
+  
+  /* But in fact from Golub 1970 Algol60 */
+
 /*   Computation of the singular values and complete orthogonal decom- */
 /* position of a real rectangular matrix A, */
 /* A = U diag (q) V^T, U^T U = V^T V =I , */
@@ -2804,7 +2849,7 @@ void svdcmp(double **a, int m, int n, double w[])
         } 
         break; 
       } 
-      if (its == 30) nrerror("no convergence in 30 dsvdcmp iterations"); 
+      if (its == 30) nrerror("no convergence in 30 svdcmp iterations"); 
       x=w[l];  /* shift from bottom 2 x 2 minor; */
       nm=k-1; 
       y=w[nm]; 
@@ -3181,6 +3226,26 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret
 	  xi[j][1]=xi[j][j+1]; /* Standard method of conjugate directions */
 	  xi[j][n]=xit[j];      /* and this nth direction by the by the average p_0 p_n */
 	}
+/*
+  Calculate a new set of orthogonal directions before repeating
+  the main loop.
+  Transpose V for SVD (minfit) (because minfit returns the right V in ULV=A):
+*/
+    printf(" Calculate a new set of orthogonal directions before repeating  the main loop.\n  Transpose V for MINFIT:...\n");
+    transpose_in_place ( n, xi );
+/*
+  SVD/MINFIT finds the singular value decomposition of V.
+
+  This gives the principal values and principal directions of the
+  approximating quadratic form without squaring the condition number.
+*/
+    printf(" SVDMINFIT finds the singular value decomposition of V. \n This gives the principal values and principal directions of the\n  approximating quadratic form without squaring the condition number...\n");
+    double *w; /* eigenvalues of principal directions */
+    w=vector(1,n);
+    svdminfit (xi, n, n, w ); /* In Brent's notation find d such that V=Q Diagonal(d) R, and Lambda=d^(-1/2) */
+    /* minfit ( n, vsmall, v, d ); */
+    /* v is overwritten with R. */
+    free_vector(w,1,n); 	
 #endif
 #ifdef LINMINORIGINAL
 #else
@@ -5298,16 +5363,16 @@ void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, doub
 #else  /* FLATSUP */
 /*  powell(p,xi,npar,ftol,&iter,&fret,func);*/
 /*   praxis ( t0, h0, n, prin, x, beale_f ); */
-  int prin=4;
-  double h0=0.25;
-#include "praxis.h"
+/*   int prin=4; */
+/*   double h0=0.25; */
+/* #include "praxis.h" */
   /* Be careful that praxis start at x[0] and powell start at p[1] */
    /* praxis ( ftol, h0, npar, prin, p, func ); */
-p1= (p+1); /*  p *(p+1)@8 and p *(p1)@8 are equal p1[0]=p[1] */
-printf("Praxis \n");
-fprintf(ficlog, "Praxis \n");fflush(ficlog);
-praxis ( ftol, h0, npar, prin, p1, func );
-printf("End Praxis\n");
+/* p1= (p+1); /\*  p *(p+1)@8 and p *(p1)@8 are equal p1[0]=p[1] *\/ */
+/* printf("Praxis \n"); */
+/* fprintf(ficlog, "Praxis \n");fflush(ficlog); */
+/* praxis ( ftol, h0, npar, prin, p1, func ); */
+/* printf("End Praxis\n"); */
 #endif  /* FLATSUP */
 
 #ifdef LINMINORIGINAL
-- 
2.43.0