'***************************************************** '* Solving a complex linear system A Z = B, * '* where A is complex matrix of size N*N, * '* B is a complex matrix of size N*M, * '* Z is a complex solution matrix of size N*M. * '* ------------------------------------------------- * '* Ref.: "Mathematiques en Turbo-Pascal By M. Ducamp * '* and A. Reverchon (vol 2), Eyrolles, Paris, 1988" * '* [BIBLI 05]. * '* ------------------------------------------------- * '* SAMPLE RUN: * '* * '* Solve complex linear system: * '* * '* z1 + z2 + iz3 = 5 * '* z1 + iz2 + z3 = -2i * '* iz1 + z2 + z3 = 1 * '* * '* Here: N=3, M=1 * '* * '* System solution: * '* * '* z1 = 1.500000 i -0.500000 * '* z2 = 1.000000 i 1.000000 * '* z3 = -0.500000 i -2.500000 * '* * '* DET = 20 * '* * '* Basic version by J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '***************************************************** ' Explanations: ' The original complex linear system to solve is: ' ( a11 a12 ... a1N ) ( z1 ) ( b1 ) ' ( a21 a22 ... a2N ) ( z2 ) ( b2 ) ' ( ... ... ... ... ) (... ) = (... ) ' ( aN1 aN2 ... aNN ) ( zN ) ( bN ) ' a k,l = c k,l + i d k,l ' z k = x k + i y k ' b k = e k + i f k ' where c, d, e, f, x, y are real numbers. ' The system is replaced by the following REAL system of size 2*N: ' ( c11 -d11 c12 -d12 ... c1N -d1N ) ( x1 ) ( e1 ) ' ( d11 c11 d12 c12 ... d1N c1N ) ( y1 ) ( f1 ) ' ( c21 -d21 c22 -d22 ... c2N -d2N ) ( x2 ) ( e2 ) ' ( d21 c21 d22 c22 ... d2N c2N ) ( y2 ) = ( f2 ) ' ( ... .... ... .... ... ... .... ) ( ...) ( ...) ' ( cN1 -dN1 cN2 -dN2 ... cNN -dNN ) ( xN ) ( eN ) ' ( dN1 cN1 dN2 cN2 ... dNN cNN ) ( yN ) ( fN ) ' The real system is then solved by a classic Gauss-Jordan method. ' '----------------------------------------------------------------- defdbl a-h,o-z defint i-n NMAX=3 : MMAX=1 Dim Xm(NMAX,NMAX), Ym(NMAX,NMAX) Dim Xv(NMAX), Yv(NMAX) 'used by subroutine 1000 Dim A(2*NMAX,2*NMAX) Dim B(2*NMAX,MMAX) Dim PC(2*NMAX), PL(2*NMAX), CS(2*NMAX) li=NMAX: ico=MMAX 'read Xm matrix real part of original A matrix Xm(1,1)= 1: Xm(1,2)=1: Xm(1,3)=0 ' 1 1 0 Xm(2,1)= 1: Xm(2,2)=0: Xm(2,3)=1 ' 1 0 1 Xm(3,1)= 0: Xm(3,2)=1: Xm(3,3)=1 ' 0 1 1 'read Ym matrix imaginary part of original A matrix Ym(1,1)= 0: Ym(1,2)=0: Ym(1,3)=1 ' 0 0 1 Ym(2,1)= 0: Ym(2,2)=1: Ym(2,3)=0 ' 0 1 0 Ym(3,1)= 1: Ym(3,2)=0: Ym(3,3)=0 ' 1 0 0 'read Xv vector real part of 1st column of original B matrix Xv(1)=5: Xv(2)=0: Xv(3)=1 ' 5 0 1 'read Yv vector imaginary part of 1st column of original B matrix Yv(1)=0: Yv(2)=-2: Yv(3)=0 ' 0 -2 0 'The complex system AZ=B is replaced by real system MX=V of size 2*NMAX 'create M matrix (name A for subroutine 1000) for i=1 to NMAX for j=1 to NMAX A(2*i-1,2*j-1)=Xm(i,j): A(2*i,2*j)=Xm(i,j) A(2*i,2*j-1)=Ym(i,j): A(2*i-1,2*j)=-Ym(i,j) next j next i 'create V matrix (name B for subroutine 1000) for i=1 to NMAX B(2*i-1,1)=Xv(i): B(2*i,1)=Yv(i) next ' solve real linear system MX=V of zize 2*NMAX ' by Gauss-Jordan method - See program sysmat.bas gosub 1000 'call MATINV(li*2, co, M, V, det) ' write complex solutions stored in vector V (name B) in the form: ' x1,y1,x2,y2,...,xn,yn, with zk = xk + i yk... cls print print " System solution:" print for i=1 to NMAX print using " z# = ####.###### i ####.######"; i; B(2*i-1,1); B(2*i,1) next print print " DET = "; DET print print END 'of main program '******************************************* '* SOLVING A LINEAR MATRIX SYSTEM AX = B * '* with Gauss-Jordan method using full * '* pivoting at each step. During the pro- * '* cess, original A and B matrices are * '* destroyed to spare storage location. * '* --------------------------------------- * '* INPUTS: A MATRIX N*N * '* B MATRIX N*M * '* --------------------------------------- * '* OUTPUS: A INVERSE OF A N*N * '* DET DETERMINANT OF A * '* B SOLUTION MATRIX N*M * '* --------------------------------------- * '* NOTA - If M=0, inversion of A matrix * '* only. * '******************************************* 1000 'Labels: 100, 200, 300 'Initializations N=2*li : M=ico DET = 1# FOR I = 1 TO N PC(I) = 0# PL(I) = 0# CS(I) = 0# NEXT I 'Main loop FOR K = 1 TO N 'Searching greatest pivot PV = A(K, K) IK = K: JK = K PAV = ABS(PV) FOR I = K TO N FOR J = K TO N temp = ABS(A(I, J)) IF temp > PAV THEN PV = A(I, J) PAV = ABS(PV) IK = I: JK = J END IF NEXT J NEXT I 'Search terminated, the pivot is in location I=IK, J=JK. 'Memorizing pivot location: PC(K) = JK PL(K) = IK 'DETERMINANT DET is actualised 'If DET=0, ERROR MESSAGE and STOP 'Machine dependant EPSMACH equals here 2E-16 IF IK <> K THEN DET = -DET END IF IF JK <> K THEN DET = -DET END IF DET = DET * PV temp = ABS(DET) IF temp < EPSMACH THEN print #2, print #2, " The determinant equals ZERO !!!" END 'stop program END IF 'POSITIONNING PIVOT IN K,K:} IF IK <> K THEN for I=1 to N 'EXCHANGE LINES IK and K of matrix A:} TT = A(IK, I) A(IK, I) = A(K, I) A(K, I) = TT NEXT I END IF IF M <> 0 THEN FOR I = 1 TO M TT = B(IK, I) B(IK, I) = B(K, I) B(K, I) = TT NEXT I END IF 'PIVOT is at correct line IF JK <> K THEN FOR I = 1 TO N 'EXCHANGE COLUMNS JK and K of matrix A: TT = A(I, JK) A(I, JK) = A(I, K) A(I, K)=TT NEXT I END IF 'The PIVOT is at correct column. 'and is located in K,K. 'Column K of matrix A is stored in CS vector 'then column K is set to zero. for I=1 to N CS(I) = A(I, K) A(I, K) = 0# NEXT I CS(K) = 0# A(K, K) = 1# 'Line K of matrix A is modified: temp = ABS(PV) IF temp < EPSMACH THEN print #2, print #2, " PIVOT TOO SMALL - STOP" END 'stop program END IF FOR I = 1 TO N A(K, I) = A(K, I) / PV NEXT I IF M <> 0 THEN FOR I = 1 TO M B(K, I) = B(K, I) / PV NEXT I END IF 'Other lines of matrix A are modified: FOR J = 1 TO N IF J = K THEN GOTO 100 END IF FOR I = 1 TO N 'Line J of matrix A is modified: A(J, I) = A(J, I) - CS(J) * A(K, I) NEXT I IF M <> 0 THEN FOR I = 1 TO M B(J, I) = B(J, I) - CS(J) * B(K, I) NEXT I END IF 100 NEXT J 'Line K is ready NEXT K 'end K loop 'MATRIX A INVERSION IS DONE - REARRANGEMENT OF MATRIX A 'EXCHANGE LINES FOR I = N TO 1 STEP -1 IK = INT(PC(I)) IF IK = I THEN GOTO 200 END IF 'EXCHANGE LINES I and PC(I) of matrix A: FOR J = 1 TO N TT = A(I, J) A(I, J) = A(IK, J) A(IK, J) = TT NEXT J IF M <> 0 THEN FOR J = 1 TO M TT = B(I, J) B(I, J) = B(IK, J) B(IK, J) = TT NEXT J END IF 'NO MORE EXCHANGE is NECESSARY 'Go to next line. 200 NEXT I 'End of loop i=N to 1 step -1 'EXCHANGE COLUMNS FOR J = N TO 1 STEP -1 JK = INT(PL(J)) IF JK = J THEN GOTO 300 END IF 'EXCHANGE COLUMNS J ET PL(J) of matrix A: FOR I = 1 TO N TT = A(I, J) A(I, J) = A(I, JK) A(I, JK) = TT NEXT I 'NO MORE EXCHANGE is NECESSARY 'Go to next column.} 300 NEXT J 'REARRANGEMENT TERMINATED RETURN 'End of subroutine 1000 ' end of file zmatsys.bas