[BACK]Return to imachprax.c CVS log [TXT] [DIR] Up to [Local Repository] / imach / src

Diff for /imach/src/imachprax.c between versions 1.3 and 1.4

version 1.3, 2023/06/22 11:28:07 version 1.4, 2023/06/22 12:50:51
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.4  2023/06/22 12:50:51  brouard
     Summary: stil on going
   
   Revision 1.3  2023/06/22 11:28:07  brouard    Revision 1.3  2023/06/22 11:28:07  brouard
   *** empty log message ***    *** empty log message ***
   
Line 2617  void linmin(double p[], double xi[], int Line 2620  void linmin(double p[], double xi[], int
 }   } 
   
 /**** praxis ****/  /**** 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 )  double pythag( double x, double y )
   
 /******************************************************************************/  /******************************************************************************/
Line 2657  double pythag( double x, double y ) Line 2695  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[], 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- */  /*   Computation of the singular values and complete orthogonal decom- */
 /* position of a real rectangular matrix A, */  /* position of a real rectangular matrix A, */
 /* A = U diag (q) V^T, U^T U = V^T V =I , */  /* A = U diag (q) V^T, U^T U = V^T V =I , */
Line 2807  void svdcmp(double **a, int m, int n, do Line 2852  void svdcmp(double **a, int m, int n, do
         }           } 
         break;           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; */        x=w[l];  /* shift from bottom 2 x 2 minor; */
       nm=k-1;         nm=k-1; 
       y=w[nm];         y=w[nm]; 
Line 3184  void powell(double p[], double **xi, int Line 3229  void powell(double p[], double **xi, int
           xi[j][1]=xi[j][j+1]; /* Standard method of conjugate directions */            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 */            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  #endif
 #ifdef LINMINORIGINAL  #ifdef LINMINORIGINAL
 #else  #else
Line 5301  void mlikeli(FILE *ficres,double p[], in Line 5366  void mlikeli(FILE *ficres,double p[], in
 #else  /* FLATSUP */  #else  /* FLATSUP */
 /*  powell(p,xi,npar,ftol,&iter,&fret,func);*/  /*  powell(p,xi,npar,ftol,&iter,&fret,func);*/
 /*   praxis ( t0, h0, n, prin, x, beale_f ); */  /*   praxis ( t0, h0, n, prin, x, beale_f ); */
   int prin=4;  /*   int prin=4; */
   double h0=0.25;  /*   double h0=0.25; */
 #include "praxis.h"  /* #include "praxis.h" */
   /* Be careful that praxis start at x[0] and powell start at p[1] */    /* Be careful that praxis start at x[0] and powell start at p[1] */
    /* praxis ( ftol, h0, npar, prin, p, func ); */     /* praxis ( ftol, h0, npar, prin, p, func ); */
 p1= (p+1); /*  p *(p+1)@8 and p *(p1)@8 are equal p1[0]=p[1] */  /* p1= (p+1); /\*  p *(p+1)@8 and p *(p1)@8 are equal p1[0]=p[1] *\/ */
 printf("Praxis \n");  /* printf("Praxis \n"); */
 fprintf(ficlog, "Praxis \n");fflush(ficlog);  /* fprintf(ficlog, "Praxis \n");fflush(ficlog); */
 praxis ( ftol, h0, npar, prin, p1, func );  /* praxis ( ftol, h0, npar, prin, p1, func ); */
 printf("End Praxis\n");  /* printf("End Praxis\n"); */
 #endif  /* FLATSUP */  #endif  /* FLATSUP */
   
 #ifdef LINMINORIGINAL  #ifdef LINMINORIGINAL
Removed from v.1.3  
changed lines
  Added in v.1.4

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>