/******************************************************************* * Solve a Complex Linear System By Gauss Method with full pivoting * * and correction process. * * ---------------------------------------------------------------- * * SAMPLE RUN: * * Solve Complex Linear System AX = B with: * * * * ( 47,-15) ( 62,5) ( 0,-72) (61, 20) ( 629,988) * * A = ( 6, 14) (-17,3) (-102,91) ( 7,-12) and B = (-180,825) * * ( 13, 55) ( 32,8) ( 41, 7) (25, 1) ( 877,441) * * (111,25) ( 40,0) ( 12,-82) (58,-30) ( 734,-88) * * * * System solutions: * * X1 = -1.000000 + 3.000000 I * * X2 = 2.000000 + 10.000000 I * * X3 = 7.000000 - 5.000000 I * * X4 = 17.000000 + 6.000000 I * * * * ---------------------------------------------------------------- * * Reference: "Algebre, 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,M,N; double dta, eps; Matc A; Vecc B, X; 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, Veci KP, Veci 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 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 * * (modified in the process). * * 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() /************************************************************ * Solve a Complex Linear System AX = B By Gauss Method with * * full pivoting and a correction process. * * --------------------------------------------------------- * * INPUTS: * * eps, dta : absolute and relative precisions * * M : maximum number of iterations * * N : size of linear system * * A : complex matrix * * B : right-hand side (complex vector) * * OUTPUTS: * * it: flag, =-1 if no convergence, =0 if matrix A * * is singular, =1 convergence ok. * * X : contains system solution (X1,X2,...,Xn) * * * ************************************************************/ void RSLCGTC(double eps, double dta, int M, int N, Matc A, Vecc B, int *it, Vecc X) { int I,J,L; double phi1,phi2; COMPLEX C0,C1,S,Z0; Matc C; Vecc B1, W, Z; Veci KP, LP; Z0.R=0.0; Z0.I=0.0; TSCGT(eps,N,A,it,C,KP,LP); if (*it==1) { // Save B in B1 before BSCGT for (J=1; J<=N; J++) B1[J] = B[J]; BSCGT(N,C,B,KP,LP, X); // Restore B after BSCGT for (J=1; J<=N; J++) B[J] = B1[J]; *it=-1; L=1; while (*it==-1 && L<=M) { phi1=0.0; for (I=1; I<=N; I++) if (CABS(X[I]) > phi1) phi1=CABS(X[I]); for (I=1; I<=N; I++) { S.R=Z0.R; S.I=Z0.I; for (J=1; J<=N; J++) { C0.R=S.R; C0.I=S.I; CMUL(A[I][J],X[J], &C1); CADD(C0,C1, &S); } CDIF(B[I],S, &W[I]); } BSCGT(N,C,W,KP,LP,Z); for (I=1; I<=N; I++) { C0.R=X[I].R; C0.I=X[I].I; CADD(C0,Z[I], &X[I]); } phi2=0.0; for (I=1; I<=N; I++) if (CABS(Z[I]) > phi2) phi2=CABS(Z[I]); if (phi2/phi1 < dta) *it=1; else L++; } // while } // if *it==1 } // RSLCGTC() void main() { N=4; A[1][1].R= 47.0; A[1][1].I=-15.0; A[1][2].R= 62.0; A[1][2].I= 5.0; A[1][3].R= 0.0; A[1][3].I=-72.0; A[1][4].R= 61.0; A[1][4].I= 20.0; A[2][1].R= 6.0; A[2][1].I= 14.0; A[2][2].R= -17.0; A[2][2].I= 3.0; A[2][3].R=-102.0; A[2][3].I= 91.0; A[2][4].R= 7.0; A[2][4].I=-12.0; A[3][1].R=13.0; A[3][1].I=-55.0; A[3][2].R=32.0; A[3][2].I= 8.0; A[3][3].R=41.0; A[3][3].I= 7.0; A[3][4].R=25.0; A[3][4].I= 1.0; A[4][1].R=111.0; A[4][1].I= 25.0; A[4][2].R= 40.0; A[4][2].I= 0.0; A[4][3].R= 12.0; A[4][3].I=-82.0; A[4][4].R= 58.0; A[4][4].I=-30.0; B[1].R= 629.0; B[1].I=988.0; B[2].R=-180.0; B[2].I=825.0; B[3].R= 877.0; B[3].I=441.0; B[4].R= 734.0; B[4].I=-88.0; eps=1e-10; dta=1e-10; M=1; RSLCGTC(eps,dta,M,N,A,B,&it,X); printf("\n System solutions:\n"); for (I=1; I<=N; I++) { printf(" X%d = %f", I, X[I].R); if (X[I].I >= 0.0) printf(" + "); else printf(" - "); printf("%f I\n", fabs(X[I].I)); } printf("\n"); } // end of file trslcgtc.cpp