Marsyas
0.6.0-alpha
|
00001 /****************************************************** 00002 * LU decomposition routines used by test_lu.cpp * 00003 * with dynamic allocations * 00004 * * 00005 * C++ version by J-P Moreau, Paris * 00006 * --------------------------------------------------- * 00007 * Reference: * 00008 * * 00009 * "Numerical Recipes by W.H. Press, B. P. Flannery, * 00010 * S.A. Teukolsky and W.T. Vetterling, Cambridge * 00011 * University Press, 1986". * 00012 * --------------------------------------------------- * 00013 * Uses: basis_r.cpp and vmblock.cpp * 00014 ******************************************************/ 00015 00016 #include "lu.h" 00017 00018 /************************************************************** 00019 * Given an N x N matrix A, this routine replaces it by the LU * 00020 * decomposition of a rowwise permutation of itself. A and N * 00021 * are input. INDX is an output vector which records the row * 00022 * permutation effected by the partial pivoting; D is output * 00023 * as -1 or 1, depending on whether the number of row inter- * 00024 * changes was even or odd, respectively. This routine is used * 00025 * in combination with LUBKSB to solve linear equations or to * 00026 * invert a matrix. Return code is 1, if matrix is singular. * 00027 **************************************************************/ 00028 int LUDCMP(REAL **A, int n, int *INDX, int *d) { 00029 00030 REAL AMAX,DUM, SUM; 00031 int I,IMAX,J,K; 00032 IMAX = 0; 00033 00034 REAL *VV; 00035 void *vmblock = NULL; 00036 00037 vmblock = vminit(); 00038 VV = (REAL *) vmalloc(vmblock, VEKTOR, NMAX, 0); 00039 00040 if (! vmcomplete(vmblock)) 00041 { 00042 // LogError ("No Memory", 0, __FILE__, __LINE__); 00043 return -1; 00044 } 00045 00046 *d=1; 00047 00048 for (I=1; I<n+1; I++) { 00049 AMAX=0.0; 00050 for (J=1; J<n+1; J++) 00051 if (ABS(A[I][J]) > AMAX) AMAX=ABS(A[I][J]); 00052 00053 if(AMAX < TINY) 00054 return 1; //singular Matrix! 00055 VV[I] = 1.0 / AMAX; 00056 } // i loop 00057 00058 for (J=1; J<n+1; J++) { 00059 00060 for (I=1; I<J; I++) { 00061 SUM = A[I][J]; 00062 for (K=1; K<I; K++) 00063 SUM = SUM - A[I][K]*A[K][J]; 00064 A[I][J] = SUM; 00065 } // i loop 00066 AMAX = 0.0; 00067 00068 for (I=J; I<n+1; I++) { 00069 SUM = A[I][J]; 00070 for (K=1; K<J; K++) 00071 SUM = SUM - A[I][K]*A[K][J]; 00072 A[I][J] = SUM; 00073 DUM = VV[I]*ABS(SUM); 00074 if (DUM >= AMAX) { 00075 IMAX = I; 00076 AMAX = DUM; 00077 } 00078 } // i loop 00079 00080 if (J != IMAX) { 00081 for (K=1; K<n+1; K++) { 00082 DUM = A[IMAX][K]; 00083 A[IMAX][K] = A[J][K]; 00084 A[J][K] = DUM; 00085 } // k loop 00086 *d = -*d; 00087 VV[IMAX] = VV[J]; 00088 } 00089 00090 INDX[J] = IMAX; 00091 00092 if (ABS(A[J][J]) < TINY) A[J][J] = TINY; 00093 00094 if (J != n) { 00095 DUM = 1.0 / A[J][J]; 00096 for (I=J+1; I<n+1; I++) 00097 A[I][J] *= DUM; 00098 } 00099 } // j loop 00100 00101 free(vmblock); 00102 return 0; 00103 00104 } // subroutine LUDCMP 00105 00106 00107 /***************************************************************** 00108 * Solves the set of N linear equations A . X = B. Here A is * 00109 * input, not as the matrix A but rather as its LU decomposition, * 00110 * determined by the routine LUDCMP. INDX is input as the permuta-* 00111 * tion vector returned by LUDCMP. B is input as the right-hand * 00112 * side vector B, and returns with the solution vector X. A, N and* 00113 * INDX are not modified by this routine and can be used for suc- * 00114 * cessive calls with different right-hand sides. This routine is * 00115 * also efficient for plain matrix inversion. * 00116 *****************************************************************/ 00117 void LUBKSB(REAL **A, int n, int *INDX, REAL *B) { 00118 REAL SUM; 00119 int I,II,J,LL; 00120 00121 II = 0; 00122 00123 for (I=1; I<n+1; I++) { 00124 LL = INDX[I]; 00125 SUM = B[LL]; 00126 B[LL] = B[I]; 00127 if (II != 0) 00128 for (J=II; J<I; J++) 00129 SUM = SUM - A[I][J]*B[J]; 00130 else if (SUM != 0.0) II = I; 00131 B[I] = SUM; 00132 } // i loop 00133 00134 for (I=n; I>0; I--) { 00135 SUM = B[I]; 00136 if (I < n) { 00137 for (J=I+1; J<n+1; J++) 00138 SUM = SUM - A[I][J]*B[J]; 00139 } 00140 B[I] = SUM / A[I][I]; 00141 } // i loop 00142 00143 00144 } // LUBKSB 00145 00146 // end of file lu.cpp 00147 00148