--- imach/src/imachprax.c	2023/06/22 11:28:07	1.3
+++ imach/src/imachprax.c	2023/06/22 12:50:51	1.4
@@ -1,6 +1,9 @@
-/* $Id: imachprax.c,v 1.3 2023/06/22 11:28:07 brouard Exp $
+/* $Id: imachprax.c,v 1.4 2023/06/22 12:50:51 brouard Exp $
   $State: Exp $
   $Log: imachprax.c,v $
+  Revision 1.4  2023/06/22 12:50:51  brouard
+  Summary: stil on going
+
   Revision 1.3  2023/06/22 11:28:07  brouard
   *** empty log message ***
 
@@ -1377,12 +1380,12 @@ double gnuplotversion=GNUPLOTVERSION;
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imachprax.c,v 1.3 2023/06/22 11:28:07 brouard Exp $ */
+/* $Id: imachprax.c,v 1.4 2023/06/22 12:50:51 brouard Exp $ */
 /* $State: Exp $ */
 #include "version.h"
 char version[]=__IMACH_VERSION__;
 char copyright[]="April 2023,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2022";
-char fullversion[]="$Revision: 1.3 $ $Date: 2023/06/22 11:28:07 $"; 
+char fullversion[]="$Revision: 1.4 $ $Date: 2023/06/22 12:50:51 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -2617,6 +2620,41 @@ void linmin(double p[], double xi[], int
 } 
 
 /**** 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 )
 
 /******************************************************************************/
@@ -2657,9 +2695,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 , */
@@ -2807,7 +2852,7 @@ void svdcmp(double **a, int m, int n, do
         } 
         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]; 
@@ -3184,6 +3229,26 @@ void powell(double p[], double **xi, int
 	  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
@@ -5301,16 +5366,16 @@ void mlikeli(FILE *ficres,double p[], in
 #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