/************************************************************* * Solving a complex homogeneous linear system by Gauss * * Method with full pivoting. * * ---------------------------------------------------------- * * SAMPLE RUN: * * | -1 i -i | * * Solve AX = 0 with A = | -i -1 1 | * * | i 1 -1 | * * * * Basic Family of system solutions: * * Solution 1 * * X1 = 0.000000 + 1.000000 I * * X2 = 1.000000 + 0.000000 I * * X3 = 0.000000 + 0.000000 I * * * * Solution 2 * * X1 = 0.000000 - 1.000000 I * * X2 = 0.000000 + 0.000000 I * * X3 = 1.000000 + 0.000000 I * * * * ---------------------------------------------------------- * * 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]; int i,l,N,M0,R0; double eps; Matc A, VX; 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; } /********************************************************* * Procedure RSHCGT solves the complex homogeneous linear * * system AX = 0 by Gauss method with full pivoting. * * ------------------------------------------------------ * * INPUTS: * * eps: required precision * * N : size of complex matrix A (NxN) * * OUTPUTS: * * R0: rank of A in output * * M0: dimension of system solutions space * * (if R0 <> N) * * VX: contains in output the M0 solution * * vectors (stored in columns 1..M0) when * * R0<>N, or the unique solution vector * * (stored in column 1) when R0=N. * *********************************************************/ void RSHCGT(double eps, int N, Matc A, int *R0, int *M0, Matc VX) { int i,j,k, k0,l,l0,n0; COMPLEX C0,C1, P0,Q0,S, Z0,Z1; int I0[NMAX]; Z0.R=0.0; Z0.I=0.0; Z1.R=1.0; Z1.I=0.0; for (k=1; k<=N; k++) I0[k]=k; *R0=0; k=1; while (*R0==0 && k<=N) { P0.R=A[k][k].R; P0.I=A[k][k].I; l0=k; k0=k; for (i=k; i<=N; i++) for (j=k; j<=N; j++) if (CABS(A[i][j]) > CABS(P0)) { P0.R=A[i][j].R; P0.I=A[i][j].I; l0=i; k0=j; } if (CABS(P0)=*R0+1; i--) if (i==*R0+l) { VX[i][l].R=Z1.R; VX[i][l].I=Z1.I; } else { VX[i][l].R=Z0.R; VX[i][l].I=Z0.I; } for (i=*R0; i>0; 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(A[i][j],VX[j][l],&C1); CADD(C0,C1,&S); } CDIV(S,A[i][i],&C0); VX[i][l].R=-C0.R; VX[i][l].I=-C0.I; } } for (i=1; i= 0.0) printf(" + "); else printf(" - "); printf("%f I\n", fabs(VX[i][l].I)); } } } printf("\n"); } // end of file rshcgt.cpp