100 REM COPYRIGHT (C) 1986 ADAM FRITZ, 133 MAIN ST., AFTON, N.Y. 13730
110 REM
120 REM   PROGRAM: MAIN PROGRAM FOR PAN-REIF TEST
130 REM
140 REM   VERSION: MICROSOFT BASIC 1.1      DATE: 04/25/86
150 REM
160 REM   DESCRIPTION:
170 REM
180 REM     INITIALIZE
190 REM        GET CONTROL PARAMETERS
200 REM     GENERATE AND DISPLAY TEST MATRIX
210 REM     COMPUTE AND DISPLAY INVERSE MATRIX
220 REM     COMPUTE AND DISPLAY RESIDUALS
230 REM
240 REM   AUTHOR:
250 REM
260 REM      ADAM FRITZ
270 REM      133 MAIN STREET
280 REM      AFTON, NEW YORK 13730
290 REM
300 REM   REFERENCE:
310 REM
320 REM     THE INVERSION OF LARGE MATRICES                   
330 REM     THOMAS E. PHIPPS, JR.                             
340 REM     BYTE MAGAZINE                                     
350 REM     APRIL, 1986                                       
360 REM
370 REM   VARIABLES:
380 REM
390 REM      LDA    - LEADING DIMENSION OF A
400 REM      N      - ORDER OF A
410 REM
420 REM      A      - INPUT: REAL SINGLE PRECISION MATRIX
430 REM      B      - OUTPUT: REAL SINGLE PRECISION MATRIX INVERSE
440 REM	     R      - RESIDUALS MATRIX
450 REM	     X	    - WORKING STORAGE MATRIX
460 REM
470 REM      ITERX   - NUMBER OF ITERATIONS
480 REM      RESIDX  - CONVERGENCE CRITERION
490 REM
500 REM      LIN    - FORTRAN LOGICAL CONSOLE INPUT
510 REM      LOUT   - FORTRAN LOGICAL CONSOLE OUTPUT
520 REM
530       LDA = 10
540       DIM A(10,10), B(10,10), R(10,10), X(10,10)
550       ITERXX = 100
560 REM
570       PRINT
580       PRINT "***** TEST PROGRAM FOR PAN-REIF *****"
590       PRINT
600 REM
610 REM   GET PARAMETERS
620 REM
630          INPUT "DIMENSION: ", N
640       IF (N < 1 OR N > LDA) GOTO 630
650 REM
660          INPUT "ITERATIONS: ", ITERX
670       IF (ITERX < 1 OR ITERX > ITERXX) GOTO 660
680 REM
690          INPUT "RESIDUAL: ", RESIDX
700       IF (RESIDX <= 0!) GOTO 690
710 REM
720 REM   GENERATE TEST MATRIX
730 REM
740       FOR I=1 TO N
750          FOR J=1 TO N
760             A(I,J) = 1!/(I+J-1)
770          NEXT
780       NEXT
790       PRINT
800       PRINT "ORIGINAL SYSTEM (BY COLUMN)"
810       PRINT
820       IC = (N+4)/5
830       ICB = 1
840       ICE = 5
850       FOR K=1 TO IC
860         IF (ICE > N) THEN ICE = N
870         FOR I=1 TO N
880           FOR J=ICB TO ICE
890              PRINT A(I,J);
900           NEXT
910           PRINT
920         NEXT
930         ICB = ICB + 5
940         ICE = ICE + 5
950         PRINT
960       NEXT
970 REM   INITIALIZE
980 REM
990       ITER = 0
1000 REM
1010 REM   INITIALIZE INVERSE  
1020 REM
1030       RX = 0!  
1040       CX = 0!  
1050       FOR I=1 TO N  
1060          RS = 0!
1070          CS = 0!
1080          FOR J=1 TO N
1090             RS = RS + ABS(A(I,J))
1100             CS = CS + ABS(A(J,I))
1110          NEXT
1120          IF (RS > RX) THEN RX = RS  
1130          IF (CS > CX) THEN CX = CS
1140       NEXT
1150       T = 1!/(RX*CX)  
1160       FOR I=1 TO N  
1170          FOR J=1 TO N
1180             B(I,J) = T*A(J,I)
1190          NEXT
1200       NEXT
1210 REM
1220 REM   ITERATION  
1230 REM
1240 REM      COMPUTE RESIDUAL MATRIX  
1250 REM
1260          RESID = 0!
1270          FOR I=1 TO N  
1280             FOR J=1 TO N  
1290                T = 0!
1300                FOR K=1 TO N
1310                   T = T + B(I,K)*A(K,J)
1320                NEXT
1330                R(I,J) = -T  
1340                T = ABS(T)
1350                IF (I <> J) THEN IF (T > RESID) THEN RESID = T
1360             NEXT
1370             T = 1! + R(I,I)  
1380             R(I,I) = T  
1390             T = ABS(T)
1400             IF (T > RESID) THEN RESID = T
1410          NEXT
1420 REM
1430 REM      DONE?
1440 REM
1450             IF (RESID <= RESIDX OR ITER >= ITERX) GOTO 1850
1460 REM
1470             ITER = ITER + 1  
1480 REM
1490 REM         UPDATE RESIDUAL MATRIX  
1500 REM
1510             FOR I=1 TO N  
1520                T = R(I,I)  
1530                T = 1! + T  
1540 REM            T = 1.0+(1.0+(1.0+(1.0+T)*T)*T)*T    
1550                FOR J=1 TO N
1560                   R(I,J) = T*R(I,J)
1570                NEXT
1580                R(I,I) = 1! + R(I,I)
1590             NEXT
1600 REM
1610 REM         COMPUTE NEW INVERSE  
1620 REM
1630             FOR I=1 TO N
1640                FOR J=1 TO N
1650                   T = 0!
1660                   FOR K=1 TO N
1670                      T = T + R(I,K)*B(K,J)
1680                   NEXT
1690                   X(I,J) = T
1700                NEXT
1710             NEXT
1720 REM
1730 REM         UPDATE INVERSE  
1740 REM
1750             FOR I=1 TO N
1760                FOR J=1 TO N
1770                   B(I,J) = X(I,J)
1780                NEXT
1790             NEXT
1800 REM
1810 REM      TEST FOR EXIT
1820 REM
1830          GOTO 1260
1840 REM
1850       PRINT "ITERATIONS: ", ITER
1860       PRINT
1870       PRINT "INVERSE (BY COLUMN)"
1880       PRINT
1890       IC = (N+4)/5
1900       ICB = 1
1910       ICE = 5
1920       FOR K=1 TO IC
1930         IF (ICE > N) THEN ICE = N
1940         FOR I=1 TO N
1950           FOR J=ICB TO ICE
1960              PRINT B(I,J);
1970           NEXT
1980           PRINT
1990         NEXT
2000         ICB = ICB + 5
2010         ICE = ICE + 5
2020         PRINT
2030       NEXT
2040 REM
2050 REM   DISPLAY RESIDUALS
2060 REM
2070       PRINT "RESIDUALS (BY COLUMN)"
2080       PRINT
2090       IC = (N+4)/5
2100       ICB = 1
2110       ICE = 5
2120       FOR K=1 TO IC
2130         IF (ICE > N) THEN ICE = N
2140         FOR I=1 TO N
2150           FOR J=ICB TO ICE
2160              PRINT R(I,J);
2170           NEXT
2180           PRINT
2190         NEXT
2200         ICB = ICB + 5
2210         ICE = ICE + 5
2220         PRINT
2230       NEXT
2240 REM
2250       PRINT "***** END OF TEST *****"
2260 REM
2270       END
2280 REM
2290 REM COPYRIGHT (C) 1986 ADAM FRITZ, 133 MAIN ST., AFTON, N.Y. 13730