10 REM Real-entry matrix inversion (Pan-Reif method)
20 REM Program by Thomas E. Phipps, Jr.
30 REM Ref= Proc. 17th Annual ACM Sympos. on Theory
40 REM of Computing, Providence, RI, May, 1985
50 REM Random-entry matrices, single precision
60 REM Quartic modification of Newton's iteration
70 INPUT "Run number = ";RN
80 REM N data is in line 640, D1 data in line 660
90 READ N, D1:RANDOMIZE RN:OPTION BASE 1
100 REM Can optionally insert DEFDBL A,B,E,X here
110 REM A is the matrix to invert, B is the
	approximate inverse if A, E is the error
	matix for the inverse, and X is a temporary
	storage matrix as explained in the article.
120 DIM A(N,N),B(N,N),E(N,N),X(N,N)
130 REM Generate random-entry A matrix. For real
	data, replace 140 by an INPUT statement
	and input normalized data (each element
	of A is divided by the largest element, L).
	Delete line 150. B and D1 must be
	multiplied by L to find the inverse and error
140 FOR I=1 TO N:FOR J=1 TO N:A(I,J)=RND(1):NEXT:NEXT
150 FOR I=1 TO N:FOR J=1 TO N:IF RND(1)<.5 THEN 
                                   A(I,J)=-A(I,J)
160 NEXT :NEXT
170 CLS:TIME$="00"   'Start clock
180 REM Eval. t by Pan-Reif eq 8 in BYTE article
190 FOR I=1 TO N:R0=0:S0=0:FOR J=1 TO N
200 R0=R0+ABS(A(I,J)):S0=S0+ABS(A(J,I)):NEXT J
210 X(1,I)=R0:E(1,I)=S0:NEXT I:T1=0:T2=0
220 FOR I=1 TO N:IF X(1,I)>T1 THEN T1=X(1,I)
230 IF E(1,I)>T2 THEN T2=E(1,I)
240 NEXT I:T=1/(T1*T2)
250 REM Eval. initial B-matrix
260 FOR I=1 TO N:FOR J=1 TO N:B(I,J)=T*A(J,I):NEXT:
                                      NEXT :H=0
270 PRINT "ITER."TAB(10);"E(1,1)"TAB(28);"E(1,N)"
271 TAB(46);"B(1,1)"TAB(64);"B(1,N)"
280 REM Eval. error matrix E
290 FOR I=1 TO N:FOR K=1 TO N:Z=0:FOR J=1 TO N
300 Z=Z+B(I,J)*A(J,K):NEXT J:E(I,K)=-Z:NEXT K:
                                       NEXT I
310 FOR I=1 TO N:E(I,I)=1+E(I,I):NEXT I
320 BEEP:PRINT H;TAB(6);E(1,1);TAB(24);E(1,N);
              TAB(42);B(1,1);TAB(60);B(1,N)
330 IF H>50 THEN PRINT "STUCK!":BEEP:BEEP:
                                  GOTO 620
340 REM Test for escape from loop
350 FOR I=1 TO N:FOR J=1 TO N:IF ABS(E(I,J))>D1
                                 THEN 370
360 NEXT :NEXT :GOTO 440
370 H=H+1  'Newton's iteration
                    (quartic modification)
380 FOR I=1 TO N:
    X(1,I)=1+(1+(1+(1+E(I,I))*E(I,I))*E(I,I))*E(I,I):
    NEXT
390 FOR I=1 TO N:FOR J=1 TO N:E(I,J)=X(1,I)*E(I,J):
                     NEXT :E(I,I)=1+E(I,I):NEXT
400 FOR I=1 TO N:FOR K=1 TO N:W=0:FOR J=1 TO N
410 W=W+E(I,J)*B(J,K):NEXT J:X(I,K)=W:NEXT K:NEXT I
420 FOR I=1 TO N:FOR J=1 TO N:B(I,J)=X(I,J):NEXT:
                                            NEXT
430 GOTO 280
440 REM output follows
450 T$=TIME$:BEEP:BEEP:BEEP:PRINT "Run number";RN,
                                 "Duration = ";T$
460 Print "No. iter.=";H,"Matrix dim.="N,
                            "Max. error=";D1
470 INPUT "Another iteration (y,n)";C$:
     IF C$="Y" OR C$="y" THEN TIME$="00":GOTO 370
480 PRINT "The maximum error was calculated by
            multiplying A^-1*A."
490 PRINT "As a check on the true error in the
           approximated inverse matrix"
500 PRINT "the maximum error in A*A^-1 will now
           be determined."
510 ERASE X      'recycle the temporary storage
                    matrix for errors in A*A^-1
520 DIM X(N,N)
530 FOR I=1 TO N:FOR J=1 TO N:FOR K=1 TO N
540 X(I,K)=X(I,K)+A(I,J)*B(J,K)
550 NEXT :NEXT :NEXT
560 REM You may print out the matrix and its
        inverse by placing a PRINT command for
        A and B in the following loop.
569 MAX.ERROR=0
570 FOR I=1 TO N:FOR J=1 TO N
580 TEST.VAL=ABS(X(I,J)):IF I=J THEN
                     TEST.VAL=ABS(1-TEST.VAL)
590 IF TEST.VAL>MAX.ERROR THEN MAX.ERROR=TEST.VAL
600 NEXT :NEXT
610 PRINT "The max. error found this way is ";
                                     MAX.ERROR
620 END
630 REM Matrix dim. N =
640 DATA 10
650 REM Error criterion D1 =
660 DATA 3E-6