/****************************************************** * LU decomposition routines used by test_lu.cpp * * with dynamic allocations * * * * C++ version by J-P Moreau, Paris * * --------------------------------------------------- * * Reference: * * * * "Numerical Recipes by W.H. Press, B. P. Flannery, * * S.A. Teukolsky and W.T. Vetterling, Cambridge * * University Press, 1986". * * --------------------------------------------------- * * Uses: basis_r.cpp and vmblock.cpp * ******************************************************/ #include #include #define TINY 2.5e-16 /************************************************************** * Given an N x N matrix A, this routine replaces it by the LU * * decomposition of a rowwise permutation of itself. A and N * * are input. INDX is an output vector which records the row * * permutation effected by the partial pivoting; D is output * * as -1 or 1, depending on whether the number of row inter- * * changes was even or odd, respectively. This routine is used * * in combination with LUBKSB to solve linear equations or to * * invert a matrix. Return code is 1, if matrix is singular. * **************************************************************/ int LUDCMP(REAL **A, int n, int *INDX, int *d) { REAL AMAX,DUM, SUM; int I,IMAX,J,K; REAL *VV; void *vmblock = NULL; vmblock = vminit(); VV = (REAL *) vmalloc(vmblock, VEKTOR, n, 0); if (! vmcomplete(vmblock)) { LogError ("No Memory", 0, __FILE__, __LINE__); return 1; } *d=1; for (I=1; I AMAX) AMAX=ABS(A[I][J]); if(AMAX < TINY) return 1; VV[I] = 1.0 / AMAX; } // i loop for (J=1; J= AMAX) { IMAX = I; AMAX = DUM; } } // i loop if (J != IMAX) { for (K=1; K0; I--) { SUM = B[I]; if (I < n) { for (J=I+1; J