/**************************************************** * Solving a complex linear system A Z = B, * * where A is complex matrix of size N*N, * * B is a complex matrix of size N*M, * * Z is a complex solution matrix of size N*M. * * ------------------------------------------------- * * Ref.: "Mathématiques en Turbo-Pascal By M. Ducamp * * and A. Reverchon (vol 2), Eyrolles, Paris, 1988" * * [BIBLI 05]. * * ------------------------------------------------- * * SAMPLE RUN: * * * * Solve complex linear system: * * * * z1 + z2 + iz3 = 5 * * z1 + iz2 + z3 = -2i * * iz1 + z2 + z3 = 1 * * * * Here: N=3, M=1 * * * * System solution: * * * * z1 = 1.500000 i -0.500000 * * z2 = 1.000000 i 1.000000 * * z3 = -0.500000 i -2.500000 * * * * DET = 2.00000000000000E+0001 * * * * C++ version by J-P Moreau, Paris. * * (www.jpmoreau.fr) * ***************************************************** Explanations: The original complex linear system to solve is: [ a11 a12 ... a1N ] [ z1 ] [ b1 ] [ a21 a22 ... a2N ] [ z2 ] [ b2 ] [ ... ... ... ... ] [... ] = [... ] [ aN1 aN2 ... aNN ] [ zN ] [ bN ] a k][l = c k][l + i d k][l z k = x k + i y k b k = e k + i f k where c, d, e, f, x, y are real numbers. The system is replaced by the following REAL system of size 2*N: [ c11 -d11 c12 -d12 ... c1N -d1N ] [ x1 ] [ e1 ] [ d11 c11 d12 c12 ... d1N c1N ] [ y1 ] [ f1 ] [ c21 -d21 c22 -d22 ... c2N -d2N ] [ x2 ] [ e2 ] [ d21 c21 d22 c22 ... d2N c2N ] [ y2 ] = [ f2 ] [ ... .... ... .... ... ... .... ] [ ...] [ ...] [ cN1 -dN1 cN2 -dN2 ... cNN -dNN ] [ xN ] [ eN ] [ dN1 cN1 dN2 cN2 ... dNN cNN ] [ yN ] [ fN ] The real system is then solved by a classic Gauss-Jordan method. */ #include #include #define MACH_EPS 2e-16 #define NMAX 3 #define MMAX 1 typedef double MAT[NMAX][NMAX]; typedef double VEC[NMAX]; typedef double MAT2[2*NMAX][2*NMAX]; typedef double MAT3[2*NMAX][MMAX]; MAT Xm, Ym; VEC Xv, Yv; MAT2 M; MAT3 V; int i,j,li,co; double det; /****************************************** * SOLVING A LINEAR MATRIX SYSTEM AX = B * * with Gauss-Jordan method using full * * pivoting at each step. During the pro- * * cess, original AA and BB matrices are * * destroyed to spare storage location. * * --------------------------------------- * * INPUTS: AA MATRIX N*N * * BB MATRIX N*M * * --------------------------------------- * * OUTPUS: AA INVERSE OF AA N*N * * DET DETERMINANT OF AA * * BB SOLUTION MATRIX N*M * * --------------------------------------- * * NOTA - If M=0 inversion of AA matrix * * only. * ******************************************/ void MATINV(int N, int M, MAT2 AA, MAT3 BB, double *DET) { double PC[2*NMAX],PL[2*NMAX],CS[2*NMAX]; double PV,PAV,temp,TT; int I,IK,J,JK,K; // Initializations : *DET= 1.0; for (I=0; I PAV) { PV=AA[I][J]; PAV= fabs(PV); IK=I; JK=J; } } // Search terminated, the pivot is in location I=IK, J=JK. // Memorizing pivot location: PC[K]=JK; PL[K]=IK; // DETERMINANT DET is actualised // If DET=0, ERROR MESSAGE and STOP // Machine dependant EPSMACH equals here 2e-16 if (IK!=K) *DET=-*DET; if (JK!=K) *DET=-*DET; *DET=*DET*PV; temp= fabs(*DET); if (temp < MACH_EPS) { printf("\n The determinant equals ZERO !!!\n"); return; } // POSITIONNING PIVOT IN K,K: if(IK!=K) for (I=0; I=0; I--) { IK=(int) PC[I]; if (IK==I) goto fin1; // EXCHANGE LINES I and PC(I) of matrix AA: for (J=0; J=0; J--) { JK=(int) PL[J]; if (JK==J) goto fin2; // EXCHANGE COLUMNS J ET PL(J) of matrix AA: for (I=0; I