DECLARE SUB ICGT (eps!, N%, it%) DECLARE SUB MATMUL (N%) DECLARE FUNCTION CABS! (Z!()) DECLARE SUB CDIV (Z1!(), Z2!(), Z!()) DECLARE SUB CMUL (Z1!(), Z2!(), Z!()) DECLARE SUB TSCGT (eps!, N%, it%, KP%(), LP%()) DECLARE SUB BSCGT (N%, KP%(), LP%()) '**************************************************************** '* calculate inverse of a complex square matrix by Gauss Method * '* with full pivoting. * '* ------------------------------------------------------------ * '* SAMPLE RUN: * '* Calculate inverse of complex square matrix: * '* * '* | (0, 1) ( 2,0) (-1,-1) | * '* A = | (1, 0) ( 1,1) ( 0,-3) | * '* | (0,-2) (-1,1) ( 3, 0) | * '* * '* Inverse matrix: * '* ( 0.0000, 0.0000) ( 0.3333, 0.0000) ( 0.0000, 0.3333) * '* ( 0.7500, 0.0000) ( -0.1667, -0.0833) ( 0.3333, 0.0833) * '* ( 0.2500, -0.2500) ( -0.0833, 0.2500) ( 0.2500, -0.0833) * '* * '* Product AM1*A: * '* ( 1.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) * '* ( 0.0000, 0.0000) ( 1.0000, 0.0000) ( 0.0000, 0.0000) * '* ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 1.0000, 0.0000) * '* * '* ------------------------------------------------------------ * '* Ref.: "Algebre, Algorithmes et programmes en Pascal * '* By Jean-Louis Jardrin, DUNOD Paris, 1988". * '* * '* Quick Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '**************************************************************** 'Program InvMatC DEFINT I-N OPTION BASE 1 DIM SHARED NMAX NMAX = 20 DIM SHARED A(NMAX, NMAX, 2), AM1(NMAX, NMAX, 2), MI(NMAX, NMAX, 2) DIM SHARED C(NMAX, NMAX, 2), W(NMAX, 2), Z(NMAX, 2) N = 3 'size of complex matrix A A(1, 1, 1) = 0!: A(1, 1, 2) = 1! A(1, 2, 1) = 2!: A(1, 2, 2) = 0! A(1, 3, 1) = -1!: A(1, 3, 2) = -1! A(2, 1, 1) = 1!: A(2, 1, 2) = 0! A(2, 2, 1) = 1!: A(2, 2, 2) = 1! A(2, 3, 1) = 0!: A(2, 3, 2) = -3! A(3, 1, 1) = 0!: A(3, 1, 2) = -2! A(3, 2, 1) = -1!: A(3, 2, 2) = 1! A(3, 3, 1) = 3!: A(3, 3, 2) = 0! eps = 1E-10 ICGT eps, N, it F\$ = " (###.####,###.####)" CLS PRINT PRINT " Inverse matrix:" FOR i = 1 TO N FOR j = 1 TO N PRINT USING F\$; AM1(i, j, 1); AM1(i, j, 2); NEXT j PRINT NEXT i ' Check AM1 * A = I MATMUL N PRINT PRINT " Product AM1*A:" FOR i = 1 TO N FOR j = 1 TO N PRINT USING F\$; MI(i, j, 1); MI(i, j, 2); NEXT j PRINT NEXT i PRINT END 'of main program '************************************************************** '* 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 * '* KP: contains the column exchanges. * '* LP: contains the line exchanges. * '* OUTPUT: * '* Z : system solution complex vector. * '************************************************************** 'Note: C, W, Z are shared with main program and ICGT SUB BSCGT (N, KP(), LP()) DIM C0(2), C1(2), S(2), Z0(2)'Complex numbers DIM TMP(2), TMP1(2) Z0(1) = 0!: Z0(2) = 0! FOR k = 1 TO N - 1 l0 = LP(k) IF l0 <> k THEN 'CSwap(W(k),W(l0)) C0(1) = W(l0, 1): C0(2) = W(l0, 2) W(l0, 1) = W(k, 1): W(l0, 2) = W(k, 2) W(k, 1) = C0(1): W(k, 2) = C0(2) END IF FOR i = k + 1 TO N C0(1) = W(i, 1): C0(2) = W(i, 2) 'CMUL(C(i,k),W(k),C1) TMP(1) = C(i, k, 1): TMP(2) = C(i, k, 2) TMP1(1) = W(k, 1): TMP1(2) = W(k, 2) CMUL TMP(), TMP1(), C1() 'CDIF(C0,C1,W(i)) W(i, 1) = C0(1) - C1(1) W(i, 2) = C0(2) - C1(2) NEXT i NEXT k 'CDIV(W(N),C(N,N),Z(N)) TMP(1) = W(N, 1): TMP(2) = W(N, 2) TMP1(1) = C(N, N, 1): TMP1(2) = C(N, N, 2) CDIV TMP(), TMP1(), C1() Z(N, 1) = C1(1): Z(N, 2) = C1(2) FOR i = N - 1 TO 1 STEP -1 S(1) = Z0(1): S(2) = Z0(2) FOR j = i + 1 TO N C0(1) = S(1): C0(2) = S(2) 'CMUL(C(i,j),Z(j),C1) TMP(1) = C(i, j, 1): TMP(2) = C(i, j, 2) TMP1(1) = Z(j, 1): TMP1(2) = Z(j, 2) CMUL TMP(), TMP1(), C1() 'CADD(C0,C1,S) S(1) = C0(1) + C1(1): S(2) = C0(2) + C1(2) NEXT j 'CDIF(W(i),S,C0) C0(1) = W(i, 1) - S(1): C0(2) = W(i, 2) - S(2) 'CDIV(C0,C(i,i),Z(i)) TMP(1) = C(i, i, 1): TMP(2) = C(i, i, 2) CDIV C0(), TMP(), TMP1() Z(i, 1) = TMP1(1): Z(i, 2) = TMP1(2) NEXT i FOR k = N - 1 TO 1 STEP -1 k0 = KP(k) IF k0 <> k THEN 'CSwap(Z(k),Z(k0)) TMP(1) = Z(k0, 1): TMP(2) = Z(k0, 2) Z(k0, 1) = Z(k, 1): Z(k0, 2) = Z(k, 2) Z(k, 1) = TMP(1): Z(k, 2) = TMP(2) END IF NEXT k END SUB 'BSCGT FUNCTION CABS (Z()) CABS = SQR(Z(1) ^ 2 + Z(2) ^ 2) END FUNCTION SUB CADD (Z1(), Z2(), Z()) Z(1) = Z1(1) + Z2(1) Z(2) = Z1(2) + Z2(2) END SUB SUB CDIF (Z1(), Z2(), Z()) Z(1) = Z1(1) - Z2(1) Z(2) = Z1(2) - Z2(2) END SUB SUB CDIV (Z1(), Z2(), Z()) DIM C1(2) D = Z2(1) ^ 2 + Z2(2) ^ 2 IF D < 2E-12 THEN PRINT #2, " Complex Divide by zero!" ELSE C1(1) = Z2(1): C1(2) = -Z2(2) CMUL Z1(), C1(), Z() Z(1) = Z(1) / D: Z(2) = Z(2) / D END IF END SUB SUB CMUL (Z1(), Z2(), Z()) Z(1) = Z1(1) * Z2(1) - Z1(2) * Z2(2) Z(2) = Z1(1) * Z2(2) + Z1(2) * Z2(1) END SUB '**************************************************************** '* calculate inverse of a complex square matrix by Gauss Method * '* with full pivoting. * '* ------------------------------------------------------------ * '* INPUTS: * '* eps : required precision * '* N : size of complex matrix A * '* A : input complex matrix (NxN) * '* OUTPUTS: * '* it : =0, matrix is singular; =1, matris is regular * '* AM1 : inverse complex matrix (NxN). * '**************************************************************** ' A, AM1, C, W, Z are shared with main program and BSCGT SUB ICGT (eps, N, it) DIM Z0(2), Z1(2) DIM KP(NMAX), LP(NMAX) Z0(1) = 0!: Z0(2) = 0! Z1(1) = 1!: Z1(2) = 0! TSCGT eps, N, it, KP(), LP() IF it = 1 THEN FOR l = 1 TO N FOR i = 1 TO N IF i = l THEN W(i, 1) = Z1(1): W(i, 2) = Z1(2) ELSE W(i, 1) = Z0(1): W(i, 2) = Z0(2) END IF NEXT i BSCGT N, KP(), LP() FOR i = 1 TO N AM1(i, l, 1) = Z(i, 1) AM1(i, l, 2) = Z(i, 2) NEXT i NEXT l END IF END SUB 'ICGT SUB MATMUL (N) '******************************************* '* MULTIPLICATION OF TWO SQUARE COMPLEX * '* MATRICES * '* --------------------------------------- * '* INPUTS: AM1 MATRIX N*N * '* A MATRIX N*N * '* N INTEGER * '* --------------------------------------- * '* OUTPUTS: MI MATRIX N*N PRODUCT AM1*A * '* * '******************************************* DIM SUM(2), PROD(2), TMP(2), TMP1(2)'Complex numbers FOR i = 1 TO N FOR j = 1 TO N SUM(1) = 0!: SUM(2) = 0! FOR k = 1 TO N 'SUM=SUM+AM1(I,K)*A(K,J) TMP(1) = AM1(i, k, 1): TMP(2) = AM1(i, k, 2) TMP1(1) = A(k, j, 1): TMP1(2) = A(k, j, 2) CMUL TMP(), TMP1(), PROD() 'CADD(SUM,PROD,SUM) SUM(1) = SUM(1) + PROD(1) SUM(2) = SUM(2) + PROD(2) NEXT k MI(i, j, 1) = SUM(1) MI(i, j, 2) = SUM(2) NEXT j NEXT i END SUB '***************************************************************** '* 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. * '***************************************************************** SUB TSCGT (eps, N, it, KP(), LP()) DIM C0(2), C1(2), P0(2), TMP(2), TMP1(2) 'Complex numbers FOR i = 1 TO N FOR j = 1 TO N C(i, j, 1) = A(i, j, 1) C(i, j, 2) = A(i, j, 2) NEXT j NEXT i it = 1: k = 1 WHILE it = 1 AND k < N P0(1) = C(k, k, 1): P0(2) = C(k, k, 2) l0 = k: k0 = k FOR i = k TO N FOR j = k TO N C0(1) = C(i, j, 1): C0(2) = C(i, j, 2) IF CABS(C0()) > CABS(P0()) THEN P0(1) = C(i, j, 1): P0(2) = C(i, j, 2) l0 = i: k0 = j END IF NEXT j NEXT i LP(k) = l0: KP(k) = k0 IF CABS(P0()) < eps THEN it = 0 ELSE IF l0 <> k THEN FOR j = k TO N 'CSwap(C(k,j),C(l0,j)) C0(1) = C(l0, j, 1): C0(2) = C(l0, j, 2) C(l0, j, 1) = C(k, j, 1): C(l0, j, 2) = C(k, j, 2) C(k, j, 1) = C0(1): C(k, j, 2) = C0(2) NEXT j END IF IF k0 <> k THEN FOR i = 1 TO N 'CSwap(C(i,k),C(i,k0)) C0(1) = C(i, k0, 1): C0(2) = C(i, k0, 2) C(i, k0, 1) = C(i, k, 1): C(i, k0, 2) = C(i, k, 2) C(i, k, 1) = C0(1): C(i, k, 2) = C0(2) NEXT i END IF FOR i = k + 1 TO N C0(1) = C(i, k, 1): C0(2) = C(i, k, 2) 'CDIV(C0,P0,C(i,k)) CDIV C0(), P0(), TMP() C(i, k, 1) = TMP(1): C(i, k, 2) = TMP(2) FOR j = k + 1 TO N C0(1) = C(i, j, 1): C0(2) = C(i, j, 2) 'CMUL(C(i,k),C(k,j),C1) TMP(1) = C(i, k, 1): TMP(2) = C(i, k, 2) TMP1(1) = C(k, j, 1): TMP1(2) = C(k, j, 2) CMUL TMP(), TMP1(), C1() 'CDIF(C0,C1,C(i,j)) C(i, j, 1) = C0(1) - C1(1) C(i, j, 2) = C0(2) - C1(2) NEXT j NEXT i k = k + 1 END IF WEND C0(1) = C(N, N, 1): C0(2) = C(N, N, 2) IF it = 1 AND CABS(C0()) < eps THEN it = 0 END SUB 'TSCGT