'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
DATA 0, 1, -2, 1, 0, 0
DATA 3, 2, 2, 0, 1, 0
DATA 5, 4, 3, 0, 0, 1