DECLARE SUB DCGT (eps!, N%, det!()) DECLARE FUNCTION CABS! (Z!()) DECLARE SUB CDIF (Z1!(), Z2!(), Z!()) DECLARE SUB CDIV (Z1!(), Z2!(), Z!()) DECLARE SUB CMUL (Z1!(), Z2!(), Z!()) DECLARE SUB TSCGT (eps!, N%, it%, KP%(), LP%()) '************************************************************* '* 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 = 1.741656E+07 - 1.059832E+07 I * '* * '* --------------------------------------------------------- * '* 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 CDetMat DEFINT I-N OPTION BASE 1 DIM SHARED NMAX NMAX = 20 DIM SHARED A(NMAX, NMAX, 2) 'Complex matrix A DIM SHARED C(NMAX, NMAX, 2) 'Complex matrix C DIM det(2) 'Complex determinant of A ' Example #1 ' N = 3 {size of complex matrix A} ' A(1,1,1)=1.0: A(1,1,2)=0.0 ' A(1,2,1)=0.0: A(1,2,2)=1.0 ' A(1,3,1)=0.0: A(1,3,2)=1.0 ' A(2,1,1)=0.0: A(2,1,2)=1.0 ' A(2,2,1)=1.0: A(2,2,2)=0.0 ' A(2,3,1)=1.0: A(2,3,2)=0.0 ' A(3,1,1)= 0.0: A(3,1,2)=1.0 ' A(3,2,1)= 1.0: A(3,2,2)=0.0 ' A(3,3,1)=-1.0: A(3,3,2)=0.0 ' ( Det = -4.0 + 0 I ) *) ' Example #2 N = 4 A(1, 1, 1) = 47!: A(1, 1, 2) = -15! A(1, 2, 1) = 62!: A(1, 2, 2) = 5! A(1, 3, 1) = 0!: A(1, 3, 2) = -72! A(1, 4, 1) = 61!: A(1, 4, 2) = 20! A(2, 1, 1) = 6!: A(2, 1, 2) = 14! A(2, 2, 1) = -17!: A(2, 2, 2) = 3! A(2, 3, 1) = -102!: A(2, 3, 2) = 91! A(2, 4, 1) = 7!: A(2, 4, 2) = -12! A(3, 1, 1) = 13!: A(3, 1, 2) = -55! A(3, 2, 1) = 32!: A(3, 2, 2) = 8! A(3, 3, 1) = 41!: A(3, 3, 2) = 7! A(3, 4, 1) = 25!: A(3, 4, 2) = 1! A(4, 1, 1) = 111!: A(4, 1, 2) = 25! A(4, 2, 1) = 40!: A(4, 2, 2) = 0! A(4, 3, 1) = 12!: A(4, 3, 2) = -82! A(4, 4, 1) = 58!: A(4, 4, 2) = -30! eps = 1E-10 'required precision DCGT eps, N, det() CLS PRINT PRINT " Det = "; det(1); IF det(2) = 0! THEN PRINT " + "; ELSE PRINT " - "; END IF PRINT ABS(det(2)); " I" PRINT END 'of main program 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 '**************************************************************** '* 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. * '**************************************************************** SUB DCGT (eps, N, det()) DIM C0(2), C1(2), Z0(2), Z1(2) 'Complex numbers 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 = 0 THEN det(1) = Z0(1): det(2) = Z0(2) ELSE det(1) = Z1(1): det(2) = Z1(2) FOR k = 1 TO N C0(1) = det(1): C0(2) = det(2) C1(1) = C(k, k, 1): C1(2) = C(k, k, 2) CMUL C0(), C1(), det() NEXT k l = 0 FOR k = 1 TO N - 1 IF LP(k) <> k THEN l = l + 1 IF KP(k) <> k THEN l = l + 1 NEXT k IF (l MOD 2) <> 0 THEN det(1) = -det(1): det(2) = -det(2) END IF END IF END SUB 'DCGT '***************************************************************** '* 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), T0(2), TMP(2), TMP1(2) 'Complex numbers ' C=A 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 = 1 TO N T0(1) = C(i, j, 1): T0(2) = C(i, j, 2) IF CABS(T0()) > 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 T0(1) = C(k, j, 1): T0(2) = C(k, j, 2) C(k, j, 1) = C(l0, j, 1): C(k, j, 2) = C(l0, j, 2) C(l0, j, 1) = T0(1): C(l0, j, 2) = T0(2) NEXT j END IF IF k0 <> k THEN FOR i = 1 TO N T0(1) = C(i, k, 1): T0(2) = C(i, k, 2) C(i, k, 1) = C(i, k0, 1): C(i, k, 2) = C(i, k0, 2) C(i, k0, 1) = T0(1): C(i, k0, 2) = T0(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(), T0() C(i, k, 1) = T0(1): C(i, k, 2) = T0(2) FOR j = k + 1 TO N C0(1) = C(i, j, 1): C0(2) = C(i, j, 2) 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(), T0() C(i, j, 1) = T0(1): C(i, j, 2) = T0(2) NEXT j NEXT i k = k + 1 END IF WEND T0(1) = C(N, N, 1): T0(2) = C(N, N, 2) IF it = 1 AND CABS(T0()) < eps THEN it = 0 END SUB 'TSCGT