'********************************************************* '* SOLVING A LINEAR MATRIX SYSTEM AX=B * '* by Gauss-Jordan method * '* ----------------------------------------------------- * '* Basic version by J-P Moreau, Paris * '* (www.jpmoreau.fr) * '* * '* INPUTS from example file sysmat.dat, * '* OUTPUTS to sysmat.lst. * '* ----------------------------------------------------- * '* SAMPLE RUN: * '* Input file * '* 3 * '* 2.0 -1.0 1.0 * '* 1.0 5.0 -2.0 * '* 3.0 -2.0 3.0 * '* 4 * '* 5.0 -1.0 1.0 2.0 * '* 1.0 2.0 3.0 4.0 * '* 3.0 7.556 4.0 4.0 * '* * '* Output file * '* ----------------------------------------------------- * '* SOLVING THE MATRIX LINEAR SYSTEM AX = B * '* ----------------------------------------------------- * '* (BY GAUSS-JORDAN METHOD) * '* * '* N=3 * '* * '* MATRIX A: * '* * '* 2.000000 -1.000000 1.000000 * '* 1.000000 5.000000 -2.000000 * '* 3.000000 -2.000000 3.000000 * '* * '* M=4 * '* * '* MATRIX B: * '* * '* 5.000000 -1.000000 1.000000 2.000000 * '* 1.000000 2.000000 3.000000 4.000000 * '* 3.000000 7.556000 4.000000 4.000000 * '* * '* INVERSE OF MATRIX A: * '* * '* 0.785714 0.071429 -0.214286 * '* -0.642857 0.214286 0.357143 * '* -1.214286 0.071429 0.785714 * '* * '* SOLUTION MATRIX X: * '* * '* 3.357143 -2.262000 0.142857 1.000000 * '* -1.928571 3.770000 1.428571 1.000000 * '* -3.642857 7.294000 2.142857 1.000000 * '* * '* DETERMINANT= 14.000000 * '* * '* * '* VERIFICATION OF AX = B: * '* * '* 5.000000 -1.000000 1.000000 2.000000 * '* 1.000000 2.000000 3.000000 4.000000 * '* 3.000000 7.556000 4.000000 4.000000 * '* * '* * '* ----------------------------------------------------- * '* Subroutines used: * '* 1000 MATINV: Solve a linear system AX = B by the * '* Gauss-Jordan method with full pivoting. * '* 2000 Multiply two square matrices of size (n,n). * '* 2100 Multiply a matrix (n,n) by a matrix (n,m). * '* 5000 Write a header with caption to output file. * '* 5100 Write a footer to output file. * '********************************************************* DEFINT I-N DEFDBL A-H, O-Z EPSMACH = 2E-16 NMAX = 10 MMAX = 5 DIM A(NMAX, NMAX), A1(NMAX, NMAX), D(NMAX, NMAX) DIM B(NMAX, MMAX), C(NMAX, MMAX) DIM PC(NMAX), CL(NMAX), CS(NMAX) 'READ LEFT HAND MATRIX FROM INPUT FILE AND ECHO TO OUTPUT FILE OPEN "sysmat.dat" FOR INPUT AS #1 OPEN "sysmat.lst" FOR OUTPUT AS #2 INPUT #1, N 'Print header nom\$ = "SOLVING THE MATRIX LINEAR SYSTEM AX = B" : GOSUB 5000 PRINT #2, " (BY GAUSS-JORDAN METHOD)" PRINT #2, PRINT #2, " N= "; N PRINT #2, FOR I = 1 TO N FOR J = 1 TO N INPUT #1, temp A(I, J) = temp NEXT J NEXT I PRINT #2, " MATRIX A:" PRINT #2, FOR I = 1 TO N FOR J = 1 TO N PRINT #2, USING " ##.######"; A(I, J); IF J = N THEN PRINT #2, END IF NEXT J NEXT I 'READ FROM INPUT FILE RIGHT HAND MATRIX IF M>0 INPUT #1, M PRINT #2, PRINT #2, " M="; M PRINT #2, IF M <> 0 THEN FOR I = 1 TO N FOR J = 1 TO M INPUT #1, temp B(I, J) = temp NEXT J NEXT I PRINT #2, " MATRIX B:" PRINT #2, FOR I = 1 TO N FOR J = 1 TO M PRINT #2, USING " ##.######"; B(I, J); IF J = M THEN PRINT #2, END IF NEXT J NEXT I END IF 'END SECTION READ DATA CLOSE #1 'STORE MATRIX A IN MATRIX A1 FOR I = 1 TO N FOR J = 1 TO N A1(I, J) = A(I, J) NEXT J NEXT I 'CALL MATRIX INVERSION ROUTINE GOSUB 1000 'PRINT RESULTS TO OUTPUT FILE PRINT #2, PRINT #2, " INVERSE OF MATRIX A:" PRINT #2, FOR I = 1 TO N FOR J = 1 TO N PRINT #2, USING " ##.######"; A(I, J); IF J = N THEN PRINT #2, END IF NEXT J NEXT I IF M <> 0 THEN PRINT #2, PRINT #2, " SOLUTION MATRIX X:" PRINT #2, FOR I = 1 TO N FOR J = 1 TO M PRINT #2, USING " ##.######"; B(I, J); IF J = M THEN PRINT #2, END IF NEXT J NEXT I END IF PRINT #2, print #2, USING " DETERMINANT= ##.######"; DET PRINT #2, 'VERIFY THAT NEW PRODUCT A*B EQUALS ORIGINAL B MATRI ' (OPTIONAL) IF M <> 0 THEN PRINT #2, print #2, " VERIFICATION OF AX = B:" PRINT #2, 'CALL MULTIPLICATION ROUTINE 1 GOSUB 2100 FOR I = 1 TO N FOR J = 1 TO M PRINT #2, USING " ##.######"; C(I, J); IF J = M THEN PRINT #2, END IF NEXT J NEXT I ELSE PRINT #2, PRINT #2, " VERIFICATION OF A1 * A = I:" PRINT #2, 'CALL MULTIPLICATION ROUTINE 0 GOSUB 2000 FOR I = 1 TO N FOR J = 1 TO N PRINT #2, USING " ##.######"; D(I, J); IF J = N THEN print #2, END IF NEXT J NEXT I END IF 'Print footer GOSUB 5100 CLOSE #2 PRINT PRINT " Results in file sysmat.lst..." input R\$ 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 'Subroutine MATINV 'Labels: 100, 200, 300 'Initializations 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 updated '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 '******************************************* '* MULTIPLICATION OF TWO SQUARE REAL * '* MATRICES * '* --------------------------------------- * '* INPUTS: A1 MATRIX N*N * '* B MATRIX N*N * '* N INTEGER * '* --------------------------------------- * '* OUTPUTS: C MATRIX N*N PRODUCT A*B * '* * '******************************************* 2000 'Start multiplication FOR I = 1 TO N FOR J = 1 TO N SUM = 0# FOR K = 1 TO N SUM = SUM + A1(I, K) * B(K, J) C(I, J) = SUM NEXT K NEXT J NEXT I RETURN '******************************************* '* MULTIPLICATION OF TWO REAL MATRICES * '* --------------------------------------- * '* INPUTS: A1 MATRIX N*N * '* B MATRIX N*M * '* N INTEGER * '* M INTEGER * '* --------------------------------------- * '* OUTPUTS: C MATRIX N*M PRODUCT A*B * '* * '******************************************* 2100 'Start multiplication FOR I = 1 TO N FOR J = 1 TO M SUM = 0# FOR K = 1 TO N SUM = SUM + A1(I, K) * B(K, J) C(I, J) = SUM NEXT K NEXT J NEXT I RETURN 5000 'Write header PRINT #2, "------------------------------------------------------" PRINT #2, nom\$ PRINT #2, "------------------------------------------------------" RETURN 5100 'Write footer PRINT #2, PRINT #2, "------------------------------------------------------" RETURN 'End of file sysmat.bas j-p moreau april 2000