home *** CD-ROM | disk | FTP | other *** search
/ ftp.barnyard.co.uk / 2015.02.ftp.barnyard.co.uk.tar / ftp.barnyard.co.uk / cpm / walnut-creek-CDROM / CPM / EDUCATIN / PR11.LBR / PR11.BQS / PR11.BAS
BASIC Source File  |  2000-06-30  |  6KB  |  221 lines
  1. 100 REM COPYRIGHT (C) 1986 ADAM FRITZ, 133 MAIN ST., AFTON, N.Y. 13730
  2. 110 REM
  3. 120 REM   PROGRAM: MAIN PROGRAM FOR PAN-REIF TEST
  4. 130 REM
  5. 140 REM   VERSION: MICROSOFT BASIC 1.1      DATE: 04/25/86
  6. 150 REM
  7. 160 REM   DESCRIPTION:
  8. 170 REM
  9. 180 REM     INITIALIZE
  10. 190 REM        GET CONTROL PARAMETERS
  11. 200 REM     GENERATE AND DISPLAY TEST MATRIX
  12. 210 REM     COMPUTE AND DISPLAY INVERSE MATRIX
  13. 220 REM     COMPUTE AND DISPLAY RESIDUALS
  14. 230 REM
  15. 240 REM   AUTHOR:
  16. 250 REM
  17. 260 REM      ADAM FRITZ
  18. 270 REM      133 MAIN STREET
  19. 280 REM      AFTON, NEW YORK 13730
  20. 290 REM
  21. 300 REM   REFERENCE:
  22. 310 REM
  23. 320 REM     THE INVERSION OF LARGE MATRICES                   
  24. 330 REM     THOMAS E. PHIPPS, JR.                             
  25. 340 REM     BYTE MAGAZINE                                     
  26. 350 REM     APRIL, 1986                                       
  27. 360 REM
  28. 370 REM   VARIABLES:
  29. 380 REM
  30. 390 REM      LDA    - LEADING DIMENSION OF A
  31. 400 REM      N      - ORDER OF A
  32. 410 REM
  33. 420 REM      A      - INPUT: REAL SINGLE PRECISION MATRIX
  34. 430 REM      B      - OUTPUT: REAL SINGLE PRECISION MATRIX INVERSE
  35. 440 REM         R      - RESIDUALS MATRIX
  36. 450 REM         X        - WORKING STORAGE MATRIX
  37. 460 REM
  38. 470 REM      ITERX   - NUMBER OF ITERATIONS
  39. 480 REM      RESIDX  - CONVERGENCE CRITERION
  40. 490 REM
  41. 500 REM      LIN    - FORTRAN LOGICAL CONSOLE INPUT
  42. 510 REM      LOUT   - FORTRAN LOGICAL CONSOLE OUTPUT
  43. 520 REM
  44. 530       LDA = 10
  45. 540       DIM A(10,10), B(10,10), R(10,10), X(10,10)
  46. 550       ITERXX = 100
  47. 560 REM
  48. 570       PRINT
  49. 580       PRINT "***** TEST PROGRAM FOR PAN-REIF *****"
  50. 590       PRINT
  51. 600 REM
  52. 610 REM   GET PARAMETERS
  53. 620 REM
  54. 630          INPUT "DIMENSION: ", N
  55. 640       IF (N < 1 OR N > LDA) GOTO 630
  56. 650 REM
  57. 660          INPUT "ITERATIONS: ", ITERX
  58. 670       IF (ITERX < 1 OR ITERX > ITERXX) GOTO 660
  59. 680 REM
  60. 690          INPUT "RESIDUAL: ", RESIDX
  61. 700       IF (RESIDX <= 0!) GOTO 690
  62. 710 REM
  63. 720 REM   GENERATE TEST MATRIX
  64. 730 REM
  65. 740       FOR I=1 TO N
  66. 750          FOR J=1 TO N
  67. 760             A(I,J) = 1!/(I+J-1)
  68. 770          NEXT
  69. 780       NEXT
  70. 790       PRINT
  71. 800       PRINT "ORIGINAL SYSTEM (BY COLUMN)"
  72. 810       PRINT
  73. 820       IC = (N+4)/5
  74. 830       ICB = 1
  75. 840       ICE = 5
  76. 850       FOR K=1 TO IC
  77. 860         IF (ICE > N) THEN ICE = N
  78. 870         FOR I=1 TO N
  79. 880           FOR J=ICB TO ICE
  80. 890              PRINT A(I,J);
  81. 900           NEXT
  82. 910           PRINT
  83. 920         NEXT
  84. 930         ICB = ICB + 5
  85. 940         ICE = ICE + 5
  86. 950         PRINT
  87. 960       NEXT
  88. 970 REM   INITIALIZE
  89. 980 REM
  90. 990       ITER = 0
  91. 1000 REM
  92. 1010 REM   INITIALIZE INVERSE  
  93. 1020 REM
  94. 1030       RX = 0!  
  95. 1040       CX = 0!  
  96. 1050       FOR I=1 TO N  
  97. 1060          RS = 0!
  98. 1070          CS = 0!
  99. 1080          FOR J=1 TO N
  100. 1090             RS = RS + ABS(A(I,J))
  101. 1100             CS = CS + ABS(A(J,I))
  102. 1110          NEXT
  103. 1120          IF (RS > RX) THEN RX = RS  
  104. 1130          IF (CS > CX) THEN CX = CS
  105. 1140       NEXT
  106. 1150       T = 1!/(RX*CX)  
  107. 1160       FOR I=1 TO N  
  108. 1170          FOR J=1 TO N
  109. 1180             B(I,J) = T*A(J,I)
  110. 1190          NEXT
  111. 1200       NEXT
  112. 1210 REM
  113. 1220 REM   ITERATION  
  114. 1230 REM
  115. 1240 REM      COMPUTE RESIDUAL MATRIX  
  116. 1250 REM
  117. 1260          RESID = 0!
  118. 1270          FOR I=1 TO N  
  119. 1280             FOR J=1 TO N  
  120. 1290                T = 0!
  121. 1300                FOR K=1 TO N
  122. 1310                   T = T + B(I,K)*A(K,J)
  123. 1320                NEXT
  124. 1330                R(I,J) = -T  
  125. 1340                T = ABS(T)
  126. 1350                IF (I <> J) THEN IF (T > RESID) THEN RESID = T
  127. 1360             NEXT
  128. 1370             T = 1! + R(I,I)  
  129. 1380             R(I,I) = T  
  130. 1390             T = ABS(T)
  131. 1400             IF (T > RESID) THEN RESID = T
  132. 1410          NEXT
  133. 1420 REM
  134. 1430 REM      DONE?
  135. 1440 REM
  136. 1450             IF (RESID <= RESIDX OR ITER >= ITERX) GOTO 1850
  137. 1460 REM
  138. 1470             ITER = ITER + 1  
  139. 1480 REM
  140. 1490 REM         UPDATE RESIDUAL MATRIX  
  141. 1500 REM
  142. 1510             FOR I=1 TO N  
  143. 1520                T = R(I,I)  
  144. 1530                T = 1! + T  
  145. 1540 REM            T = 1.0+(1.0+(1.0+(1.0+T)*T)*T)*T    
  146. 1550                FOR J=1 TO N
  147. 1560                   R(I,J) = T*R(I,J)
  148. 1570                NEXT
  149. 1580                R(I,I) = 1! + R(I,I)
  150. 1590             NEXT
  151. 1600 REM
  152. 1610 REM         COMPUTE NEW INVERSE  
  153. 1620 REM
  154. 1630             FOR I=1 TO N
  155. 1640                FOR J=1 TO N
  156. 1650                   T = 0!
  157. 1660                   FOR K=1 TO N
  158. 1670                      T = T + R(I,K)*B(K,J)
  159. 1680                   NEXT
  160. 1690                   X(I,J) = T
  161. 1700                NEXT
  162. 1710             NEXT
  163. 1720 REM
  164. 1730 REM         UPDATE INVERSE  
  165. 1740 REM
  166. 1750             FOR I=1 TO N
  167. 1760                FOR J=1 TO N
  168. 1770                   B(I,J) = X(I,J)
  169. 1780                NEXT
  170. 1790             NEXT
  171. 1800 REM
  172. 1810 REM      TEST FOR EXIT
  173. 1820 REM
  174. 1830          GOTO 1260
  175. 1840 REM
  176. 1850       PRINT "ITERATIONS: ", ITER
  177. 1860       PRINT
  178. 1870       PRINT "INVERSE (BY COLUMN)"
  179. 1880       PRINT
  180. 1890       IC = (N+4)/5
  181. 1900       ICB = 1
  182. 1910       ICE = 5
  183. 1920       FOR K=1 TO IC
  184. 1930         IF (ICE > N) THEN ICE = N
  185. 1940         FOR I=1 TO N
  186. 1950           FOR J=ICB TO ICE
  187. 1960              PRINT B(I,J);
  188. 1970           NEXT
  189. 1980           PRINT
  190. 1990         NEXT
  191. 2000         ICB = ICB + 5
  192. 2010         ICE = ICE + 5
  193. 2020         PRINT
  194. 2030       NEXT
  195. 2040 REM
  196. 2050 REM   DISPLAY RESIDUALS
  197. 2060 REM
  198. 2070       PRINT "RESIDUALS (BY COLUMN)"
  199. 2080       PRINT
  200. 2090       IC = (N+4)/5
  201. 2100       ICB = 1
  202. 2110       ICE = 5
  203. 2120       FOR K=1 TO IC
  204. 2130         IF (ICE > N) THEN ICE = N
  205. 2140         FOR I=1 TO N
  206. 2150           FOR J=ICB TO ICE
  207. 2160              PRINT R(I,J);
  208. 2170           NEXT
  209. 2180           PRINT
  210. 2190         NEXT
  211. 2200         ICB = ICB + 5
  212. 2210         ICE = ICE + 5
  213. 2220         PRINT
  214. 2230       NEXT
  215. 2240 REM
  216. 2250       PRINT "***** END OF TEST *****"
  217. 2260 REM
  218. 2270       END
  219. 2280 REM
  220. 2290 REM COPYRIGHT (C) 1986 ADAM FRITZ, 133 MAIN ST., AFTON, N.Y. 13730
  221.