/********************************************************************************* * Test Program of units Ucmat and Clu concerning complex matrices * * ------------------------------------------------------------------------------ * * SAMPLE RUNS: * * (Complex matrix A(2,3) is defined as: * * (3,1) (1,0.5) (-7,0) * * (2,0) (4,0.0) ( 1,0) * * Complex matrix B(3,2) is defined as: * * (0.25,0.5) (2,2) * * (6, 0.0) ( 0,0) * * (1, 0.0) (-3,0) * * Find C(2,2) = A * B ) * * * * A * B = * * -0.75 4.75 25.00 8.00 * * 25.50 1.00 1.00 4.00 * * * * (Complex matrix A(4,4) is defined as: * * (1,1) ( 0,0.5) (2,0) ( 1,0) * * (0,0) (-3,0.0) (0,0) (-5,0) * * (2,0) (-5,0.0) (0,0) ( 2,0) * * (0,0) ( 9,0.0) (0,0) ( 0,0) * * Find A power 5 ) * * * * A ^ p = * * 12.00 76.00 -1543.50 488.00 24.00 48.00 141.00 -623.50 * * 0.00 0.00 -13608.00 0.00 0.00 0.00 -4455.00 0.00 * * 24.00 48.00 -123.00 275.00 24.00 40.00 8974.00 221.00 * * 0.00 0.00 8019.00 -0.00 0.00 0.00 -10935.00 0.00 * * * * (Complex matrix A(4,4) is defined as: * * (2,0.5) (1,0) (1,0) (4,0) * * (0,0.0) (1,0) (3,0) (2,0) * * (3,0.0) (1,0) (2,0) (1,0) * * (1,0.0) (1,0) (0,0) (3,7) * * Find complex determinant of A ) * * * * Det(A) = * * -12.50 26.00 * * * * (Complex matrix A(4,4) is defined as: * * ( 8,0) (2.00,0) (3,0) (12.0,0) * * ( 2,0) (4.00,0) (7,0) (0.25,0) * * ( 3,0) (7,00,0) (3,0) ( 5.0,0) * * (12,0) (0.25,0) (5,0) ( 2.0,0) * * Find Gauss decomposition of A ) * * * * A decomposed = * * 12.000 0.000 0.250 0.000 5.000 0.000 2.000 0.000 * * -0.167 0.000 6.938 0.000 1.750 0.000 4.500 0.000 * * -0.250 0.000 -0.571 0.000 5.168 0.000 -2.651 0.000 * * -0.667 0.000 -0.264 0.000 0.154 0.000 9.069 0.000 * * * * (Solve linear system A.X = B with B = ((25,0) (13.25,0) (18,0) (19,25,0)) ) * * * * Solution: * * 1.000000 0.000000 1.000000 0.000000 1.000000 0.000000 1.000000 0.00000 * * * * (Find inverse of complex matrix A ) * * * * Inverse of A = * * -0.040155 0.000000 -0.085788 0.000000 0.056557 0.000000 0.110262 0.00000 * * -0.085788 0.000000 -0.062018 0.000000 0.202201 0.000000 0.016978 0.00000 * * 0.056557 0.000000 0.202201 0.000000 -0.130316 0.000000 -0.038826 0.00000 * * 0.110262 0.000000 0.016978 0.000000 -0.038826 0.000000 -0.066631 0.00000 * * * * C++ Release 1.1 By J-P Moreau Paris. * * * * Release 1.1 (08/16/07): Added function INVC() in Module clu.pas. * * Release 1.2 (08/20/07): Added verify A * Inv(A) = Unity matrix. * *********************************************************************************/ #include #include "ucmat.h" // For COMPLEX, CMat CMat A,B,C; // three complex matrices int IP[SIZE]; COMPLEX B1[SIZE]; int Error,i,j,m,n,p; COMPLEX det; void CDet(int, CMat, COMPLEX); void DECC (int, CMat, int *, int *); void INVC (int, CMat, CMat); void MulCMat(int, int, CMat,CMat, CMat); void PowCMat(int, int, CMat,CMat); void SOLC (int, CMat, int, COMPLEX *, int *); void main() { printf(" Adjust window size (then any key + Enter)..."); getchar(); printf("\n"); m=2; n=3; A[1][1][0]= 3.0; A[1][1][1]=1.0; A[1][2][0]= 1.0; A[1][2][1]=0.5; A[1][3][0]=-7.0; A[1][3][1]=0.0; A[2][1][0]= 2.0; A[2][1][1]=0.0; A[2][2][0]= 4.0; A[2][2][1]=0.0; A[2][3][0]= 1.0; A[2][3][1]=0.0; B[1][1][0]= 0.25; B[1][1][1]=0.5; B[1][2][0]= 2.0; B[1][2][1]=2.0; B[2][1][0]= 6.0; B[2][1][1]=0.0; B[2][2][0]= 0.0; B[2][2][1]=0.0; B[3][1][0]= 1.0; B[3][1][1]=0.0; B[3][2][0]=-3.0; B[3][2][1]=0.0; MulCMat(m,n, A, B, C); printf("\n A * B =\n"); for (i=1; i<=m; i++) { for (j=1; j<=m; j++) printf(" %9.2f %9.2f", C[i][j][0], C[i][j][1]); printf("\n"); } m=4; p=5; A[1][1][0]= 1.0; A[1][1][1]=1.0; A[1][2][0]= 0.0; A[1][2][1]=0.5; A[1][3][0]= 2.0; A[1][3][1]=0.0; A[1][4][0]= 1.0; A[1][4][1]=0.0; A[2][1][0]= 0.0; A[2][1][1]=0.0; A[2][2][0]=-3.0; A[2][2][1]=0.0; A[2][3][0]= 0.0; A[2][3][1]=0.0; A[2][4][0]=-5.0; A[2][4][1]=0.0; A[3][1][0]= 2.0; A[3][1][1]=0.0; A[3][2][0]=-5.0; A[3][2][1]=0.0; A[3][3][0]= 0.0; A[3][3][1]=0.0; A[3][4][0]= 2.0; A[3][4][1]=0.0; A[4][1][0]= 0.0; A[4][1][1]=0.0; A[4][2][0]= 9.0; A[4][2][1]=0.0; A[4][3][0]= 0.0; A[4][3][1]=0.0; A[4][4][0]= 0.0; A[4][4][1]=0.0; PowCMat(m,p,A,B); printf("\n A ^ p =\n"); for (i=1; i<=m; i++) { for (j=1; j<=m; j++) printf(" %9.2f %9.2f", B[i][j][0], B[i][j][1]); printf("\n"); } m=4; A[1][1][0]= 2.0; A[1][1][1]=0.5; A[1][2][0]= 1.0; A[1][2][1]=0.0; A[1][3][0]= 1.0; A[1][3][1]=0.0; A[1][4][0]= 4.0; A[1][4][1]=0.0; A[2][1][0]= 0.0; A[2][1][1]=0.0; A[2][2][0]= 1.0; A[2][2][1]=0.0; A[2][3][0]= 3.0; A[2][3][1]=0.0; A[2][4][0]= 2.0; A[2][4][1]=0.0; A[3][1][0]= 3.0; A[3][1][1]=0.0; A[3][2][0]= 1.0; A[3][2][1]=0.0; A[3][3][0]= 2.0; A[3][3][1]=0.0; A[3][4][0]= 1.0; A[3][4][1]=0.0; A[4][1][0]= 1.0; A[4][1][1]=0.0; A[4][2][0]= 1.0; A[4][2][1]=0.0; A[4][3][0]= 0.0; A[4][3][1]=0.0; A[4][4][0]= 3.0; A[4][4][1]=7.0; CDet(m,A,det); printf("\n Det(A) =\n"); printf(" %9.2f %9.2f\n", det[0], det[1]); m=4; A[1][1][0]= 8.0; A[1][1][1]=0.0; A[1][2][0]= 2.0; A[1][2][1]=0.0; A[1][3][0]= 3.0; A[1][3][1]=0.0; A[1][4][0]=12.0; A[1][4][1]=0.0; A[2][1][0]= 2.0; A[2][1][1]=0.0; A[2][2][0]= 4.0; A[2][2][1]=0.0; A[2][3][0]= 7.0; A[2][3][1]=0.0; A[2][4][0]=0.25; A[2][4][1]=0.0; A[3][1][0]= 3.0; A[3][1][1]=0.0; A[3][2][0]= 7.0; A[3][2][1]=0.0; A[3][3][0]= 3.0; A[3][3][1]=0.0; A[3][4][0]= 5.0; A[3][4][1]=0.0; A[4][1][0]=12.0; A[4][1][1]=0.0; A[4][2][0]=0.25; A[4][2][1]=0.0; A[4][3][0]= 5.0; A[4][3][1]=0.0; A[4][4][0]= 2.0; A[4][4][1]=0.0; // call decomposition procedure (only once) DECC (m, A, IP, &Error); printf("\n A decomposed =\n"); for (i=1; i<=m; i++) { for (j=1; j<=m; j++) printf(" %9.3f %9.3f", A[i][j][0], A[i][j][1]); printf("\n"); } // complex right side B1[1][0] = 25.0; B1[1][1] = 0.0; B1[2][0] = 13.25; B1[2][1] = 0.0; B1[3][0] = 18.0; B1[3][1] = 0.0; B1[4][0] = 19.25; B1[4][1] = 0.0; // Solve complex linear system SOLC (m, A, 3, B1, IP); printf("\n Solution:\n"); for (i=1; i<=m; i++) printf(" %9.6f %9.6f", B1[i][0], B1[i][1]); printf("\n"); // Restore matrix A to call inversion routine m=4; A[1][1][0]= 8.0; A[1][1][1]=0.0; A[1][2][0]= 2.0; A[1][2][1]=0.0; A[1][3][0]= 3.0; A[1][3][1]=0.0; A[1][4][0]=12.0; A[1][4][1]=0.0; A[2][1][0]= 2.0; A[2][1][1]=0.0; A[2][2][0]= 4.0; A[2][2][1]=0.0; A[2][3][0]= 7.0; A[2][3][1]=0.0; A[2][4][0]=0.25; A[2][4][1]=0.0; A[3][1][0]= 3.0; A[3][1][1]=0.0; A[3][2][0]= 7.0; A[3][2][1]=0.0; A[3][3][0]= 3.0; A[3][3][1]=0.0; A[3][4][0]= 5.0; A[3][4][1]=0.0; A[4][1][0]=12.0; A[4][1][1]=0.0; A[4][2][0]=0.25; A[4][2][1]=0.0; A[4][3][0]= 5.0; A[4][3][1]=0.0; A[4][4][0]= 2.0; A[4][4][1]=0.0; // call inversion routine (in clu.cpp) INVC(m, A, B); // inverse of matrix A is now in matrix B, // matrix A is destroyed. printf("\n Inverse of A =\n"); for (i=1; i<=m; i++) { for (j=1; j<=m; j++) printf(" %9.6f %9.6f", B[i][j][0], B[i][j][1]); printf("\n"); } // Restore matrix A to verify A * Inv(A) m=4; A[1][1][0]= 8.0; A[1][1][1]=0.0; A[1][2][0]= 2.0; A[1][2][1]=0.0; A[1][3][0]= 3.0; A[1][3][1]=0.0; A[1][4][0]=12.0; A[1][4][1]=0.0; A[2][1][0]= 2.0; A[2][1][1]=0.0; A[2][2][0]= 4.0; A[2][2][1]=0.0; A[2][3][0]= 7.0; A[2][3][1]=0.0; A[2][4][0]=0.25; A[2][4][1]=0.0; A[3][1][0]= 3.0; A[3][1][1]=0.0; A[3][2][0]= 7.0; A[3][2][1]=0.0; A[3][3][0]= 3.0; A[3][3][1]=0.0; A[3][4][0]= 5.0; A[3][4][1]=0.0; A[4][1][0]=12.0; A[4][1][1]=0.0; A[4][2][0]=0.25; A[4][2][1]=0.0; A[4][3][0]= 5.0; A[4][3][1]=0.0; A[4][4][0]= 2.0; A[4][4][1]=0.0; // call inversion routine (in clu.cpp) MulCMat(m, m, A, B, C); // C now contains product A * Inv(A) printf("\n A * Inv(A)=\n"); for (i=1; i<=m; i++) { for (j=1; j<=m; j++) printf(" %9.6f %9.6f", C[i][j][0], C[i][j][1]); printf("\n"); } } // end of file tucmat.cpp