/*************************************************************** * calculate inverse of a complex square matrix by Gauss Method * * with full pivoting. * * ------------------------------------------------------------ * * SAMPLE RUN: * * Calculate inverse of complex square matrix: * * * * | (0, 1) ( 2,0) (-1,-1) | * * A = | (1, 0) ( 1,1) ( 0,-3) | * * | (0,-2) (-1,1) ( 3, 0) | * * * * Inverse matrix: * * ( 0.0000, 0.0000) ( 0.3333, 0.0000) ( 0.0000, 0.3333) * * ( 0.7500, 0.0000) ( -0.1667, -0.0833) ( 0.3333, 0.0833) * * ( 0.2500, -0.2500) ( -0.0833, 0.2500) ( 0.2500, -0.0833) * * * * Product AM1*A: * * ( 1.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) * * ( 0.0000, -0.0000) ( 1.0000, 0.0000) ( 0.0000, 0.0000) * * ( -0.0000, -0.0000) ( -0.0000, 0.0000) ( 1.0000, 0.0000) * * * * ------------------------------------------------------------ * * Ref.: "Algèbre, Algorithmes et programmes en Pascal * * By Jean-Louis Jardrin, DUNOD Paris, 1988". * * * * C++ Release By J-P Moreau, Paris. * * (www.jpmoreau.fr) * ***************************************************************/ #include #include #define NMAX 20 //complex number typedef struct { double R,I; //algebraic form } COMPLEX; typedef COMPLEX Matc[NMAX][NMAX]; typedef COMPLEX Vecc[NMAX]; typedef int Veci[NMAX]; int i,it,j,m, N; double eps; Matc A, AM1, MI; double CABS(COMPLEX Z) { return sqrt(Z.R*Z.R + Z.I*Z.I); } void CADD(COMPLEX Z1, COMPLEX Z2, COMPLEX *Z) { Z->R = Z1.R + Z2.R; Z->I = Z1.I + Z2.I; } void CDIF(COMPLEX Z1, COMPLEX Z2, COMPLEX *Z) { Z->R = Z1.R - Z2.R; Z->I = Z1.I - Z2.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>2.2e-16) { C.R=Z2.R; C.I=-Z2.I; CMUL(Z1,C,Z); Z->R=Z->R/d; Z->I=Z->I/d; } } void CSwap(COMPLEX *Z1, COMPLEX *Z2) { COMPLEX C; C.R=Z1->R; C.I=Z1->I; Z1->R=Z2->R; Z1->I=Z2->I; Z2->R=C.R; Z2->I=C.I; } /**************************************************************** * TSCGT procedure implements the triangularization algorithm of * * Gauss with full pivoting at each step for a complex matrix, A * * and saves the made transformations in KP and LP. * * ------------------------------------------------------------- * * INPUTS: * * N: size of complex matrix A * * A: complex matrix of size N x N * * OUTPUTS; * * it: =0 if A is singular, else =1. * * C: contains the upper triangular matrix and the * * multipliers used during the process. * * KP: contains the column exchanges. * * LP: contains the line exchanges. * ****************************************************************/ void TSCGT(double eps, int N, Matc A, int *it, Matc C, int *KP, int *LP) { int i,j,k,k0,l0; COMPLEX C0,C1,P0,T0; for (i=1; i<=N; i++) for(j=1; j<=N; j++) { C[i][j].R=A[i][j].R; C[i][j].I=A[i][j].I; } *it=1; k=1; while (*it==1 && k<=N) { P0.R=C[k][k].R; P0.I=C[k][k].I; l0=k; k0=k; for (i=k; i<=N; i++) for (j=k; j<=N; j++) if (CABS(C[i][j]) > CABS(P0)) { P0.R=C[i][j].R; P0.I=C[i][j].I; l0=i; k0=j; } LP[k]=l0; KP[k]=k0; if (CABS(P0) < eps) *it=0; else { if (l0!=k) for (j=k; j<=N; j++) { T0.R=C[k][j].R; T0.I=C[k][j].I; C[k][j].R=C[l0][j].R; C[k][j].I=C[l0][j].I; C[l0][j].R=T0.R; C[l0][j].I=T0.I; } if (k0!=k) for (i=1; i<=N; i++) { T0.R=C[i][k].R; T0.I=C[i][k].I; C[i][k].R=C[i][k0].R; C[i][k].I=C[i][k0].I; C[i][k0].R=T0.R; C[i][k0].I=T0.I; } for (i=k+1; i<=N; i++) { C0.R=C[i][k].R; C0.I=C[i][k].I; CDIV(C0,P0,&C[i][k]); for (j=k+1; j<=N; j++) { C0.R=C[i][j].R; C0.I=C[i][j].I; CMUL(C[i][k],C[k][j],&C1); CDIF(C0,C1,&C[i][j]); } } k++; } } if (*it==1 && CABS(C[N][N]) < eps) *it=0; } //TSCGT() /************************************************************* * Function BSCGT calculates the solution of upper triangular * * system [S(n-1)] by the back substitution method and deduces* * from it the solution of system [S]. The call must be made * * only after a call to TSCGT and the matrix of [S] must be * * regular. * * ---------------------------------------------------------- * * INPUTS: * * N : size of matrix C * * C: contains the upper triangular matrix and the * * multipliers used during the triangularization * * process (in output of TSCGT). . * * W : contains the right-side coefficients * * KP: contains the column exchanges. * * LP: contains the line exchanges. * * OUTPUT: * * Z : system solution complex vector. * *************************************************************/ void BSCGT(int N, Matc C, Vecc W, Veci KP, Veci LP, Vecc Z) { int i,j,k,k0,l0; COMPLEX C0,C1,S,Z0; Z0.R=0.0; Z0.I=0.0; for (k=1; k0; i--) { S.R=Z0.R; S.I=Z0.I; for (j=i+1; j<=N; j++) { C0.R=S.R; C0.I=S.I; CMUL(C[i][j],Z[j],&C1); CADD(C0,C1,&S); } CDIF(W[i],S,&C0); CDIV(C0,C[i][i],&Z[i]); } for (k=N-1; k>0; k--) { k0=KP[k]; if (k0!=k) CSwap(&Z[k],&Z[k0]); } } // BSCGT() /*************************************************************** * calculate inverse of a complex square matrix by Gauss Method * * with full pivoting. * * ------------------------------------------------------------ * * INPUTS: * * eps : required precision * * N : size of complex matrix A * * A : input complex matrix (NxN) * * OUTPUTS: * * it : =0, matrix is singular; =1, matris is regular * * AM1 : inverse complex matrix (NxN). * ***************************************************************/ void ICGT(double eps, int N, Matc A, int *it, Matc AM1) { int i,l; COMPLEX Z0,Z1; Matc C; Vecc W,Z; Veci KP,LP; Z0.R=0.0; Z0.I=0.0; Z1.R=1.0; Z1.I=0.0; TSCGT(eps,N,A,it,C,KP,LP); if (*it==1) for (l=1; l<=N; l++) { for (i=1; i<=N; i++) if (i==l) { W[i].R=Z1.R; W[i].I=Z1.I; } else { W[i].R=Z0.R; W[i].I=Z0.I; } BSCGT(N,C,W,KP,LP,Z); for (i=1; i<=N; i++) { AM1[i][l].R=Z[i].R; AM1[i][l].I=Z[i].I; } } } // ICGT() void MATMUL(Matc A, Matc B, Matc C, int N) { /****************************************** * MULTIPLICATION OF TWO SQUARE COMPLEX * * MATRICES * * --------------------------------------- * * INPUTS: A MATRIX N*N * * B MATRIX N*N * * N INTEGER * * --------------------------------------- * * OUTPUTS: C MATRIX N*N PRODUCT A*B * * * ******************************************/ COMPLEX SUM,PROD, TMP; int I,J,K; for (I=1; I<=N; I++) for (J=1; J<=N; J++) { SUM.R=0.0; SUM.I=0.0; for (K=1; K<=N; K++) { // SUM=SUM+A[I][K]*B[K][J] CMUL(A[I][K],B[K][J],&PROD); CADD(SUM,PROD,&TMP); SUM.R=TMP.R; SUM.I=TMP.I; } C[I][J].R=SUM.R; C[I][J].I=SUM.I; } } void main() { N=3; // size of complex matrix A A[1][1].R= 0.0; A[1][1].I= 1.0; A[1][2].R= 2.0; A[1][2].I= 0.0; A[1][3].R=-1.0; A[1][3].I=-1.0; A[2][1].R=1.0; A[2][1].I= 0.0; A[2][2].R=1.0; A[2][2].I= 1.0; A[2][3].R=0.0; A[2][3].I=-3.0; A[3][1].R= 0.0; A[3][1].I=-2.0; A[3][2].R=-1.0; A[3][2].I= 1.0; A[3][3].R= 3.0; A[3][3].I= 0.0; eps=1e-10; ICGT(eps, N, A, &it, AM1); printf("\n Inverse matrix:\n"); for (i=1; i<=N; i++) { for (j=1; j<=N; j++) printf(" (%8.4f,%8.4f)", AM1[i][j].R, AM1[i][j].I); printf("\n"); } // Check AM1 * A = I MATMUL(AM1, A, MI, N); printf("\n Product AM1*A:\n"); for (i=1; i<=N; i++) { for (j=1; j<=N; j++) printf(" (%8.4f,%8.4f)", MI[i][j].R, MI[i][j].I); printf("\n"); } printf("\n"); } // end of file Invmatc.cpp