'GAUSS-JORDAN MATRICE INVERSE
'=============================

'A original matrix, InvA inverse of A, B extended matrix, P = A * InvA

'let's have a square matrix A (3x3)         '2,        1,    -2
                                            '3         2      2
                                            '5         4      3

'--- creating the extended matrix B = A + unit matrix

DATA         2,     1,   -2,      1,     0,     0
DATA         3,     2,    2,      0,     1,     0
DATA         5,     4,    3,      0,     0,     1

n = 3

DIM A(n+1, n+1)
DIM B(n+1, n*2)

PRINT "EXTENDED MATRIX B :": PRINT
FOR i = 1 TO n
FOR j = 1 TO 2*n
    READ m
    B(i,j) = m
    PRINT "", B(i,j);    'and displaying it
NEXT j
    PRINT: PRINT
NEXT i
PRINT

'--- A is the left part of B, extraction

PRINT "ORIGINAL MATRIX A, extracted from B :"
PRINT
FOR i = 1 TO n
FOR j = 1 TO n
    A(i,j) = B(i,j)
    PRINT "", A(i,j);    'and displaying it
NEXT j
    PRINT: PRINT
NEXT i
PRINT

'--- calculation begins here 

'i = row index        'k = other row index (not-i)      'j = column index
'intermediary variables p, q, r to avoid divisions like A(i,j)/A(m,n) e.g., giving unpredictable results

FOR i = 1 TO n

'reducing each diagonal term B(i,i)  to 1 by dividing the hole row by B(i,i)
    p = B(i, i)
    FOR j = 1 TO 2*n
    q = B(i, j)
    B(i, j) = q / p
    NEXT j

'multiply row i by B(k,i)
'substract row i from row k

        FOR k = 1 TO n
        p = B(k, i)
        IF k <> i THEN
                  FOR j = 1 TO 2*n
                  q = B(k, j)
                   r = B(i, j)
                  B(k, j) = q - p * r
                  NEXT j
        END IF
    NEXT k

NEXT i

'--- end of calculation, let's see matrix B  now

PRINT "INTERMEDIAIRY MATRIX B :": PRINT
FOR i = 1 TO n
  FOR j= 1 TO 2*n
      m = B(i, j)
      PRINT "", USING("###.###", m);
  NEXT j
      PRINT: PRINT
NEXT i
PRINT

'---INVERSE of a, i.e. InvA, is the *right* part of B

PRINT "INVERSE of A (rounded) :": PRINT
DIM InvA(n+1,n+1)

FOR i = 1 TO n
FOR j = n+1 TO 2*n                            'copying the right part of B
          InvA(i, j-n) = B(i,j)
          m = InvA(i, j-n)
          PRINT "", USING("###.###", m);      'rounding
NEXT j
          PRINT: PRINT
NEXT i
PRINT

'--- Checking the result: A * InvA must guve P = the unit matrix

PRINT "CHECKING for InvA:  A * Inv must give the unit matrix P"
PRINT

DIM P(n+1, n+1)                  'empty

FOR i = 1 TO n                    'for each row i o f A and P
  FOR j = 1 TO n                  'for each horizontal term j of row i of P
    FOR k = 1 TO n                'for each vertical term k of column j of InvA
        x = A(k,i): y = InvA(j,k)
        P(i, j) = P(i, j) + x*y
     NEXT k
     m = P(i, j)
     PRINT "", USING("###.###", m);
  NEXT j
  PRINT: PRINT
NEXT i

END