'*********************************************************
'*          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