/************************************************************ !* This program calculates the determinant of a complex * !* square matrix using Function CFindDet (see ref. below). * !* -------------------------------------------------------- * !* SAMPLE RUN: * !* * !* Calculate the deteminant of square matrix: * !* * !* ! 0 -1 -1 -1 ! * !* ! 1 0 -1 -1 ! * !* ! 1 1 0 -1 ! * !* ! 1 1 1 0 ! * !* * !* Determinant = (1.000000, 0.000000) * !* * !* -------------------------------------------------------- * !* * !* C++ Version By J-P Moreau, Paris. * !***********************************************************/ #include #include #define NMAX 25 #define FALSE 0 #define TRUE 1 //complex number typedef struct { double R,I; //algebraic form } COMPLEX; typedef COMPLEX Matc[NMAX][NMAX]; int N; Matc A; COMPLEX D; 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; } } /*--------------------------------------------------------------------------------------------------- !Function to find the determinant of a square matrix !Author : Louisda16th a.k.a Ashwith J. Rego !Description: The subroutine is based on two key points: !1) A determinant is unaltered when row operations are performed: Hence, using this principle, !row operations (column operations would work as well) are used !to convert the matrix into upper traingular form !2) The determinant of a triangular matrix is obtained by finding the product of the diagonal elements !----------------------------------------------------------------------------------------------------*/ void CFindDet(Matc matrix, int n, COMPLEX *cdet) { COMPLEX m, temp, CONE, CZERO; int i, j, k, l; bool DetExists = TRUE; l = 1; CZERO.R=0,0; CZERO.I=0.0; CONE.R=1,0; CONE.I=0.0; // Convert to upper triangular form for (k = 1; kR = CZERO.R; cdet->I = CZERO.I; return; } } for (j = k+1; j<=n; j++) { // m = matrix(j,k)/matrix(k,k) CDIV(matrix[j][k],matrix[k][k],&m); for (i = k+1; i<=n; i++) { // matrix(j,i) = matrix(j,i) - m*matrix(k,i) CMUL(m,matrix[k][i],&temp); matrix[j][i].R = matrix[j][i].R - temp.R; matrix[j][i].I = matrix[j][i].I - temp.I; } } } // Calculate determinant by finding product of diagonal elements cdet->R = 1.0; cdet->I = 0.0; for (i = 1; i<=n; i++) { // CDet = CDet * matrix(i,i) CMUL(*cdet,matrix[i][i],&temp); cdet->R = temp.R; cdet->I = temp.I; } } void main() { int N=4; A[1][1].R = 0.0; A[1][1].I = 0.0; A[2][1].R = 1.0; A[2][1].I = 0.0; A[3][1].R = 1.0; A[3][1].I = 0.0; A[4][1].R = 1.0; A[4][1].I = 0.0; A[1][2].R = -1.0; A[1][2].I = 0.0; A[2][2].R = 0.0; A[2][2].I = 0.0; A[3][2].R = 1.0; A[3][2].I = 0.0; A[4][2].R = 1.0; A[4][2].I = 0.0; A[1][3].R = -1.0; A[1][3].I = 0.0; A[2][3].R = -1.0; A[2][3].I = 0.0; A[3][3].R = 0.0; A[3][3].I = 0.0; A[4][3].R = 1.0; A[4][3].I = 0.0; A[1][4].R = -1.0; A[1][4].I = 0.0; A[2][4].R = -1.0; A[2][4].I = 0.0; A[3][4].R = -1.0; A[3][4].I = 0.0; A[4][4].R = 0.0; A[4][4].I = 0.0; CFindDet(A, N, &D); printf("\n Determinant = (%f,%f)\n\n", D.R, D.I); } // end of file cfinddet.cpp