/****************************************************** * LU decomposition routines used by test_lu.cpp * * with static allocations * * * * C++ version by J-P Moreau, Paris * * for a complex linear system. * * (www.jpmoreau.fr) * * --------------------------------------------------- * * Reference: * * * * "Numerical Recipes by W.H. Press, B. P. Flannery, * * S.A. Teukolsky and W.T. Vetterling, Cambridge * * University Press, 1986". * * * ******************************************************/ #include #include #define SIZE 25 #define TINY 2.5e-16 //complex number typedef struct { double R,I; //algebraic form } COMPLEX; typedef COMPLEX CMAT[SIZE][SIZE]; typedef COMPLEX CVEC[SIZE]; typedef int IMAT[SIZE]; double CABS(COMPLEX Z) { return sqrt(Z.R*Z.R + Z.I*Z.I); } void CMUL(COMPLEX Z1, COMPLEX Z2, COMPLEX *Z) { Z->R = Z1.R*Z2.R - Z1.I*Z2.I; Z->I = Z1.R*Z2.I + Z1.I*Z2.R; } void CDIV(COMPLEX Z1, COMPLEX Z2, COMPLEX *Z) { double d; COMPLEX C; d = Z2.R*Z2.R+Z2.I*Z2.I; if (d<1e-12) printf(" Complex Divide by zero!\n"); else { C.R=Z2.R; C.I=-Z2.I; CMUL(Z1,C,Z); Z->R=Z->R/d; Z->I=Z->I/d; } } /************************************************************** * 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(CMAT A, int n, IMAT INDX, int *d) { COMPLEX AMAX,DUM, SUM, TMP; int I,IMAX,J,K; CVEC VV; *d=1; for (I=1; I CABS(AMAX)) { AMAX.R=A[I][J].R; AMAX.I=A[I][J].I; } if(CABS(AMAX) < TINY) return 1; // VV[I] = 1.0 / AMAX; TMP.R=1.0; TMP.I=0.0; CDIV(TMP, AMAX, &VV[I]); } // i loop for (J=1; J= CABS(AMAX)) { IMAX = I; AMAX.R = DUM.R; AMAX.I = DUM.I; } } // i loop if (J != IMAX) { for (K=1; K0; I--) { SUM.R = B[I].R; SUM.I = B[I].I; if (I < n) { for (J=I+1; J