DECLARE SUB RSLCGTC (eps!, dta!, M%, N%, it%) DECLARE FUNCTION CABS! (Z!()) DECLARE SUB CADD (Z1!(), Z2!(), Z!()) DECLARE SUB CDIV (Z1!(), Z2!(), Z!()) DECLARE SUB CMUL (Z1!(), Z2!(), Z!()) DECLARE SUB BSCGT (N%, KP%(), LP%()) DECLARE SUB TSCGT (eps!, N%, it%, KP%(), LP%()) '******************************************************************** '* Solve a Complex Linear System By Gauss Method with full pivoting * '* and correction process. * '* ---------------------------------------------------------------- * '* SAMPLE RUN: * '* Solve Complex Linear System AX = B with: * '* * '* ( 47,-15) ( 62,5) ( 0,-72) (61, 20) ( 629,988) * '* A = ( 6, 14) (-17,3) (-102,91) ( 7,-12) and B = (-180,825) * '* ( 13, 55) ( 32,8) ( 41, 7) (25, 1) ( 877,441) * '* (111,25) ( 40,0) ( 12,-82) (58,-30) ( 734,-88) * '* * '* System solutions: * '* ( -1.0000, 3.0000) * '* ( 2.0000, 10.0000) * '* ( 7.0000, -5.0000) * '* ( 17.0000, 6.0000) * '* * '* ---------------------------------------------------------------- * '* Reference: "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 CSysLin DEFINT I-N OPTION BASE 1 DIM SHARED NMAX NMAX = 20 ' Matc = Array(1..NMAX,1..NMAX) of Complex ' Vecc = Array(1..NMAX) of Complex ' Veci = Array(1..NMAX) of Integer ' I,it,M,N: Integer ' dta, eps: Real ' Common complex matrices and vectors DIM SHARED A(NMAX, NMAX, 2) DIM SHARED B(NMAX, 2), X(NMAX, 2) DIM SHARED C(NMAX, NMAX, 2), W(NMAX, 2), Z(NMAX, 2) N = 4 'Size of linear system 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! B(1, 1) = 629!: B(1, 2) = 988! B(2, 1) = -180!: B(2, 2) = 825! B(3, 1) = 877!: B(3, 2) = 441! B(4, 1) = 734!: B(4, 2) = -88! eps = 1E-10: dta = 1E-10: M = 1 RSLCGTC eps, dta, M, N, it F\$ = " (###.####,###.####)" CLS PRINT PRINT " System solutions:" FOR J = 1 TO N PRINT USING F\$; X(J, 1); X(J, 2) NEXT J 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 other subroutines. 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 '************************************************************* '* Solve a Complex Linear System AX = B By Gauss Method with * '* full pivoting and a correction process. * '* --------------------------------------------------------- * '* INPUTS: * '* eps, dta : absolute and relative precisions * '* M : maximum number of iterations * '* N : size of linear system * '* A : complex matrix * '* B : right-hand side (complex vector) * '* OUTPUTS: * '* it: flag, =-1 if no convergence, =0 if matrix A * '* is singular, =1 convergence ok. * '* X : contains system solution (X1,X2,...,Xn) * '* * '************************************************************* 'Note: A, B, X are shared with main program SUB RSLCGTC (eps, dta, M, N, it) DIM C0(2), C1(2), S(2), TMP(2), TMP1(2), Z0(2) 'Complex numbers DIM KP(NMAX), LP(NMAX) Z0(1) = 0!: Z0(2) = 0! TSCGT eps, N, it, KP(), LP() IF it = 1 THEN FOR I = 1 TO N W(I, 1) = B(I, 1): W(I, 2) = B(I, 2) NEXT I BSCGT N, KP(), LP() FOR I = 1 TO N X(I, 1) = Z(I, 1) X(I, 2) = Z(I, 2) NEXT I PRINT " Vector X after first BSCGT:" FOR I = 1 TO N PRINT " "; X(I, 1); " "; X(I, 2) NEXT I INPUT " ", R\$ it = -1: L = 1 WHILE it = -1 AND L <= M phi1 = 0! FOR I = 1 TO N TMP(1) = X(I, 1): TMP(2) = X(I, 2) IF CABS(TMP()) > phi1 THEN phi1 = CABS(TMP()) NEXT I FOR I = 1 TO N S(1) = Z0(1): S(2) = Z0(2) FOR J = 1 TO N C0(1) = S(1): C0(2) = S(2) 'CMUL(A(I,J),X(J),C1) TMP(1) = A(I, J, 1): TMP(2) = A(I, J, 2) TMP1(1) = X(J, 1): TMP1(2) = X(J, 2) CMUL TMP(), TMP1(), C1() CADD C0(), C1(), S() NEXT J 'CDIF(B(I),S,W(I)) W(I, 1) = B(I, 1) - S(1) W(I, 2) = B(I, 2) - S(2) NEXT I BSCGT N, KP(), LP() FOR I = 1 TO N C0(1) = X(I, 1): C0(2) = X(I, 2) 'CADD(C0,Z(I),X(I)) X(I, 1) = C0(1) + Z(I, 1) X(I, 2) = C0(2) + Z(I, 2) NEXT I phi2 = 0! FOR I = 1 TO N TMP(1) = Z(I, 1): TMP(2) = Z(I, 2) IF CABS(TMP()) > phi2 THEN phi2 = CABS(TMP()) NEXT I IF (phi2 / phi1) < dta THEN it = 1 ELSE L = L + 1 END IF WEND END IF END SUB 'SLCGTC '***************************************************************** '* 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. * '***************************************************************** 'NOTE: C is shared with main program and other Subroutines. 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