/************************************************************ * Determinant of a complex square matrix * * By Gauss Method with full pivoting * * --------------------------------------------------------- * * SAMPLE RUN: * * Calculate the determinant of complex matrix: * * ( 47,-15) ( 62,5) ( 0,-72) (61, 20) * * ( 6, 14) (-17,3) (-102,91) ( 7,-12) * * ( 13, 55) ( 32,8) ( 41, 7) (25, 1) * * (111,25) ( 40,0) ( 12,-82) (58,-30) * * * * Det = 17416564.000000 - 10598320.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 N; double eps; Matc A; COMPLEX det; 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<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; } } /**************************************************************** * 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=1; 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() /*************************************************************** * The DCGT procedure calculates the complex determinant of a * * complex square matrix by the Gauss method with full pivoting * * ------------------------------------------------------------ * * INPUTS: * * eps: required precision * * N : size of matrix A * * A : complex matrix of size N x N * * OUTPUT: * * det: complex determinant. * ***************************************************************/ void DCGT(double eps, int N, Matc A, COMPLEX *det) { int it,k,l; COMPLEX C0,Z0,Z1; Matc C; int KP[NMAX], LP[NMAX]; 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==0) { det->R=Z0.R; det->I=Z0.I; } else { det->R=Z1.R; det->I=Z1.I; for (k=1; k<=N; k++) { C0.R=det->R; C0.I=det->I; CMUL(C0,C[k][k],det); } l=0; for (k=1; kR=-det->R; det->I=-det->I; } } } // DCGT() void main() { /* Example #1 N=3; //size of complex matrix A A[1][1].R=1.0; A[1][1].I=0.0; A[1][2].R=0.0; A[1][2].I=1.0; A[1][3].R=0.0; A[1][3].I=1.0; A[2][1].R=0.0; A[2][1].I=1.0; A[2][2].R=1.0; A[2][2].I=0.0; A[2][3].R=1.0; A[2][3].I=0.0; A[3][1].R= 0.0; A[3][1].I=1.0; A[3][2].R= 1.0; A[3][2].I=0.0; A[3][3].R=-1.0; A[3][3].I=0.0; ( Det = -4.0 + 0 I ) Example #2 */ 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; eps=1e-10; DCGT(eps, N, A, &det); printf("\n Det = %f", det.R); if (det.I>=0.0) printf(" + "); else printf(" - "); printf("%f I\n\n", fabs(det.I)); } // end of file cdetmat.cpp