Date   : Wed, 02 Feb 1983 17:08:37 PST (Wed)
From   : David Allen Gewirtz <UCBARPA.dag@ucb-vax.arpa (David Allen Gewirtz)>
Subject: Re: MATRIX INVERTER NEEDED

Hope this helps...it's in BASIC, but can be converted with no
real problem.  It's out of "Some Common BASIC Programs" by Osborne,
McGraw Hill.  (you never can tell where you find cheap and dirty
hacks)..

This inverts a square matrix.  The inversion is performed by a modified
Gauss-Jordan elimination method.

10 print "matrix inversion"
20 print
25 REM A() and B() are the matrix dimensions
30 DIM A(10,10), B(10,10)
35 REM Matrix is square, so only one dimension is needed
40 print "Dimension of Matrix";
50 input r
60 print "enter matrix elements:"
70 for j = 1 to r
80 print "Row";j
90 for i=1 to r
100 print "value column";i;
110 input a(i,j)
120 next i
130 b(j,j) = 1
140 next j
145
146 REM *** Ignore 145 - statements 150 - 420 invert the matrix
147 REM
150
(let's try again)
150 For j = 1 to r
160 for i=j to r
170 if a(i,j) <> 0 then 210      ' Ain't BASIC 'onderful?
180 print i
190 Print "singular matrix"
200 goto 500
210 for k=1 to r
220 s=a(j,k)
230 a(j,k)=a(i,k)
240 a(i,k)=s
250 s=b(j,k)
260 b(j,k)=b(i,k)
270 b(i,k)=s                   ' a swap() function would be nice
280 next k
290 t=1/a(j,j)
300 for k = 1 to r
310 a(j,k)=t*a(j,k)
320 b(j,k)=t*b(j,k)
330 next k
340 for l = 1 to r
350 if l=j then 410
360 t = -a(l,j)
370 for k=1 to r
380 a(l,k)=a(l,k)+t*a(j,k)
390 b(l,k)=b(l,k)+t*b(j,k)
400 next k
410 next l
420 next j
430 print
440 rem - print the resultant matrix
(make the above line into 439)
440 for i=1 to r
450 for j = 1 to r
459 rem - round off and print
460 print int(b(i,j)*1000+.5)/1000;" "; 'Urgh...BASIC
470 next j
480 print
490 next i
500 end

Well, it certainly isn't glamorous, but it does work and can be
hacked into Pascal.

Good luck,

David