home *** CD-ROM | disk | FTP | other *** search
/ Best Objectech Shareware Selections / UNTITLED.iso / boss / educ / math / 023 / matop.bas (.txt) < prev    next >
Encoding:
GW-BASIC  |  1990-02-27  |  23.9 KB  |  405 lines
  1. 1  '        MATRIX OPERATIONS  ---  MATOP.BAS  ---  by  Dr Russell Langley
  2. 2  GOTO 400
  3. 4  '<UNK! {000A}>--- Press Enter ---
  4. 5  IF PR THEN RETURN ELSE PRINT TAB(40);:PRINT "Press <Enter> to continue.";:IN$=INKEY$:WHILE INKEY$<>CHR$(13):WEND:LOCATE,40:PRINT SPACE$(26):RETURN
  5. 7  '<UNK! {000A}>*** Redirect to Block ***
  6. 9  ON QB GOTO 10,177,429,409,488:STOP  '=start,printout,etc - CLOSE (exc 177)<UNK! {000A}><UNK! {000A}>--- Another go? ---
  7. 10  CLOSE:IF HEAD=1 THEN LPRINT STRING$(79,61)STRING$(4,10)
  8. 11  DO$="run this program again now":GOSUB 20:IF Z$="Y" THEN 2 ELSE 30
  9. 19  '<UNK! {000A}>--- Yes/No? ---
  10. 20  PRINT:PRINT"Do you want to "+DO$;
  11. 21  INPUT" (Y/N)";Z$:IF Z$="" THEN Z$="N":RETURN ELSE Z$=CHR$(ASC(Z$) AND 95):IF Z$="Y" OR Z$="N" THEN RETURN ELSE PRINT"WHAT?  ";:GOTO 21
  12. 29  '<UNK! {000A}>--- Errors & End ---
  13. 30  IF ERR THEN BEEP ELSE RUN "MENU"
  14. 31  IF ERR=70 THEN LINE INPUT"Can't write to that disk.  Remove its Write-Protect tab, then press <Enter>.";Z$:RESUME
  15. 32  IF ERR=71 THEN LINE INPUT"That drive is empty or its gate is open.  Fix, then press <Enter>.";Z$:RESUME
  16. 33  IF ERR=210 THEN RESUME 9  'from #86
  17. 39  ON ERROR GOTO 0:END'<UNK! {000A}><UNK! {000A}>--- Messages ---
  18. 40  BEEP:PRINT "---> Sorry, that entry is illegal.":RETURN
  19. 41  BEEP:PRINT "---> Sorry, double quotes are not allowed here.":RETURN
  20. 42  BEEP:PRINT"* * * Can't Do That.":QB=4:GOTO 9
  21. 43  COLOR 23,0:PRINT:PRINT"Working";:COLOR 7,0:LOCATE ,1:RETURN
  22. 44  LOCATE,1:PRINT"Ok, done.";:GOTO 5
  23. 46  ZZ$=STRING$(37-LEN(Z$)\2,177):LOCATE 1,1:PRINT ZZ$"  ";:COLOR 15,0:PRINT Z$;:COLOR 7,0:PRINT"  "ZZ$:RETURN 'Display brightened Z$ at top of screen
  24. 49  '<UNK! {000A}>*** Vetted Decoding of FF ..C(I,J).. from X$ ***<UNK! {000A}>    Needs I, M, Q(0) from #96 or #256, & UT>0 if UT matrix.
  25. 50  K=1:L=M:IF UT>0 THEN L=M-I+1
  26. 51  KX=0:FOR J=K TO L:Y$=SPACE$(10):KY=0
  27. 52  KX=KX+1:IF KX>LEN(X$) THEN IF J=L THEN 54 ELSE 57
  28. 53  Z$=MID$(X$,KX,1):IF INSTR("-.0123456789",Z$) THEN KY=KY+1:MID$(Y$,KY,1)=Z$:GOTO 52 ELSE IF Z$<>" " THEN 58 ELSE IF KY=0 THEN 52
  29. 54  IF Q(0) THEN Q(J)=VAL(Y$) ELSE C(I,J)=VAL(Y$)
  30. 55  NEXT J:IF KX>=LEN(X$) THEN 60
  31. 56  PLAY"L8O3CO2C":PRINT"Only the first"L"values have been read in that line.  Re-do it";:GOSUB 21:IF Z$="Y" THEN 59 ELSE 60
  32. 57  PLAY"L32O4CEG>C":PRINT"Not enough values in the line above.  Please re-do whole line.":GOTO 59
  33. 58  PLAY"L16O3CEL4>B":PRINT"That line contains a `non-numeric' entry.  Please re-do whole line."
  34. 59  INPUT X$:IF RIGHT$(X$,1)<>"/" THEN 51 ELSE X$=LEFT$(X$,LEN(X$)-1):N=I+(X$=""):IF X$>"" THEN 51
  35. 60  RETURN
  36. 69  '<UNK! {000A}>--- K/b Input of all X(I,J) in FF ---<UNK! {000A}>    Needs first I, M, & if UT matrix UT>0.  Returns N.
  37. 70  PRINT"Enter data from keyboard, ";:IF M=1 THEN PRINT "pressing <Enter> after each number.":GOTO 72
  38. 71  PRINT"in Free Format, pressing <Enter> at end of each row.":IF UT>0 THEN 73
  39. 72  PRINT"Null entry duplicates previous row.  Signal `end-of-data' by entering a `/'"
  40. 73  PRINT "Row"STR$(I);:INPUT X$:IF UT>0 AND I=M THEN N=M:GOTO 75 ELSE IF X$=""THEN IF I>1 THEN FOR J=1 TO M:X(I,J)=X(I-1,J):NEXT J:I=I+1:LOCATE CSRLIN-1,POS(0)+9:PRINT"Ditto":GOTO 73
  41. 74  IF RIGHT$(X$,1)="/" THEN X$=LEFT$(X$,LEN(X$)-1):N=I+(X$=""):IF X$="" THEN 76
  42. 75  GOSUB 50:IF N=0 THEN IF I<MXR THEN I=I+1:GOTO 73 ELSE N=MXR
  43. 76  RETURN
  44. 79  '<UNK! {000A}>*** Disk Input of C(I,J), N, M, etc ***
  45. 80  QD=1:IO$="I":GOSUB 110:CLS:PRINT "Loading data from disk.":INPUT #1,FL$,VR$
  46. 81  IF LEFT$(VR$,4)<>"(RL," THEN BEEP:PRINT "Unreadable file --- not made by our Data Filer/Editor program.":GOTO 86
  47. 82  INPUT #1,DT$,ID$,N,M,UT,VN$:PRINT "Filename: "FL$,"Made: "DT$,"Version: "VR$:PRINT"ID: "ID$:PRINT:IF MT=1 THEN NA=N:MA=M ELSE NB=N:MB=M
  48. 83  PRINT"File has"N"rows of data.":IF N>MXR THEN PRINT"--- Too many!":GOTO 85
  49. 84  PRINT:PRINT"File has"M"column variables.":IF M>MXC THEN PRINT"--- Too many!" ELSE IF ABS(UT)<=1 THEN 86 ELSE PRINT"This matrix based on"ABS(UT)"persons.":GOTO 86
  50. 85  CLOSE:BEEP:GOSUB 5:ERROR 210
  51. 86  IF MT=2 THEN 476 'test legality of B before proceeding.
  52. 87  '   Select variables
  53. 88  PRINT:IF VN$="N" THEN PRINT "The"M"variables are not named.":GOTO 90
  54. 89  PRINT "Variables in file are:":FOR J=1 TO M:INPUT #1,VN$(J):PRINT J;VN$(J),:NEXT J:PRINT
  55. 90  IF M=MNC OR UT>0 THEN 100 ELSE PRINT
  56. 91  PRINT"How many filed column variables are to be IGNORED (0-"MID$(STR$(M-MNC),2)")";:INPUT ND:IF ND<0 OR ND>M-MNC THEN 91 ELSE IF ND=0 THEN 100
  57. 92  '
  58. 93  '
  59. 94  IF ND=1 THEN PRINT "Number of the variable to be IGNORED (1-"MID$(STR$(M),2)")";:INPUT X$:Q(1)=VAL(X$):IF Q(1)<1 OR Q(1)>M THEN GOSUB 40:GOTO 94 ELSE 97
  60. 95  PRINT MID$(STR$(ND),2)" numbers of variables to be IGNORED (in ascending order & Free Format):":INPUT X$:IF VAL(X$)<1 THEN GOSUB 40:GOTO 95
  61. 96  Q(0)=1:MM=M:M=ND:GOSUB 50:Q(0)=0:M=MM
  62. 97  KK=1:L=1:FOR J=1 TO M:IF J=Q(KK) THEN KK=KK+1 ELSE VN$(L)=VN$(J):L=L+1
  63. 98  NEXT J
  64. 99  '   Now read numerical data from disk
  65. 100  COLOR 23,0:PRINT"Loading numbers";:COLOR 7,0:FOR I=1 TO N:KK=1:LL=1:IF UT<=0 THEN L=M ELSE L=M-I+1
  66. 101  FOR J=1 TO L:INPUT #1,Z
  67. 102  IF J=Q(KK) THEN KK=KK+1 ELSE C(I,LL)=Z:LL=LL+1
  68. 103  NEXT J:NEXT I:CLOSE:LOCATE,1:M=M-ND:RETURN
  69. 109  '<UNK! {000A}>--- Get Filespec ---
  70. 110  IF IO$="O" AND FL$>"" THEN PRINT "Will you file this data under the name "FL$;:GOSUB 21:IF Z$="Y" THEN 115
  71. 111  LINE INPUT "Filename (I will add .DAT extension)? ";FL$:IF FL$="" THEN 111 ELSE IF MID$(FL$,2,1)=":" THEN DR$=LEFT$(FL$,1):FL$=MID$(FL$,3)
  72. 112  ER=0:FOR I=1 TO LEN(FL$):Z$=MID$(FL$,I,1):IF INSTR(" .,/\|?*:;[]+="+CHR$(34),Z$) THEN ER=1
  73. 113  NEXT I
  74. 114  IF ER=0 AND FL$>"" AND LEN(FL$)<9 THEN FL$=FL$+".DAT" ELSE BEEP:PRINT "Invalid filename.  Will you try again";:GOSUB 21:IF Z$="Y" THEN 111 ELSE 2
  75. 115  INPUT "Which drive (A,B,C,D)";DR$:IF DR$="" THEN 115
  76. 116  DR$=CHR$(ASC(DR$) AND 95):IF INSTR("ABCD",DR$)=0 THEN 115
  77. 117  INPUT "Which directory (e.g. WORK, MYDATA, or Null Entry if root) ";DR2$:IF DR2$="" THEN DR$=DR$+":" ELSE DR$=DR$+":\"+DR2$+"\"
  78. 129  '<UNK! {000A}>--- Open File, IO$= "I" or "O" ---
  79. 130  IF IO$="I" THEN 134 ELSE ON ERROR GOTO 132:OPEN DR$+FL$ FOR INPUT AS #1
  80. 131  CLOSE:DO$="<OVERWRITE> existing "+FL$:GOSUB 20:IF Z$="N" THEN 110 ELSE 133
  81. 132  RESUME 133  'OK to start new file, since FL$ not present.
  82. 133  ON ERROR GOTO 30:OPEN DR$+FL$ FOR OUTPUT AS #1:RETURN  'print #1,Q$A$Q$Q$B$Q$:close
  83. 134  ON ERROR GOTO 136:OPEN DR$+FL$ FOR INPUT AS #1  'for input
  84. 135  ON ERROR GOTO 30:RETURN   'input #1,A$,B$:close
  85. 136  PRINT FL$" not found on Drive "DR$:RESUME 137
  86. 137  GOSUB 5:ON ERROR GOTO 30:CLS:ERASE X:SHELL "DIR "+DR$+"/W":GOSUB 5
  87. 138  DIM X(MXR,MXC):GOTO 110
  88. 139  '<UNK! {000A}>*** Save data A(I,J) ***
  89. 140  IF ID$="" THEN LINE INPUT"Problem ID? ";ID$:GOTO 142
  90. 141  LINE INPUT"Problem ID (null = unchanged)? ";Z$:IF Z$>"" THEN ID$=Z$
  91. 142  Q$=CHR$(34):IO$="O":GOSUB 110
  92. 143  PRINT "Saving data on disk ..... ";:PRINT #1,Q$FL$Q$Q$VER$Q$:PRINT #1,Q$DAT$Q$Q$ID$Q$;N;M;UT;Q$VN$Q$
  93. 144  IF VN$="Y" THEN FOR J=1 TO M:PRINT #1,Q$VN$(J)Q$;:NEXT J:PRINT #1,""
  94. 145  FOR I=1 TO N  'always file full matrix in this program.
  95. 146  FOR J=1 TO M:PRINT #1,A(I,J);:NEXT J:NEXT I:CLOSE:PRINT "done."
  96. 147  DO$="backup this file under the name "+FL$:GOSUB 20:IF Z$="Y" THEN GOSUB 115:GOTO 143 ELSE RETURN
  97. 159  '<UNK! {000A}>*** Show/Print Answers ***
  98. 160  PR=0:CLOSE:OPEN "SCRN:" FOR OUTPUT AS #2:QB=-(QB<>2)*QB-(QB=2)*QBB:RETURN
  99. 161  DO$="print these results":GOSUB 20:IF Z$="N" THEN PR=0:RETURN
  100. 162  PR=1:IF ID$="" THEN LINE INPUT "Problem ID? ";ID$
  101. 163  '
  102. 164  RETURN
  103. 165  QBB=QB:QB=2:CLS:LOCATE 8,1
  104. 166  PRINT TAB(12)"KEYTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENCLOSE"
  105. 167  PRINT TAB(12)"OPEN                                                       OPEN"
  106. 168  PRINT TAB(12)"OPEN                TURN  PRINTER  ON.                     OPEN"
  107. 169  PRINT TAB(12)"OPEN                                                       OPEN"
  108. 170  PRINT TAB(12)"OPEN    Then PRESS <ENTER> to start printing ..... or ..   OPEN"
  109. 171  PRINT TAB(12)"OPEN                                                       OPEN"
  110. 172  PRINT TAB(12)"OPEN    To send Printer Codes in Basic before printing,    OPEN"
  111. 173  PRINT TAB(12)"OPEN    press <Ctrl-Break>, & start printing by GOTO 9.    OPEN"
  112. 174  PRINT TAB(12)"OPEN                                                       OPEN"
  113. 175  PRINT TAB(12)"SCREENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENLOAD":IN$=INKEY$
  114. 176  IN$=INKEY$:IF IN$<>CHR$(13) THEN 176
  115. 177  CLS:PRINT"Printing .....":LOCATE 4,1:CLOSE:OPEN "LPT1:" FOR OUTPUT AS #2
  116. 178  IF HEAD=0 THEN PRINT #2,STRING$(79,61):PRINT #2,DAT$;TAB(42-LEN(HD$)\2);HD$;TAB(73)VER$:PRINT #2,STRING$(79,61):HEAD=1
  117. 179  IF ID$<>CHR$(168) THEN PRINT #2,CHR$(10)"Problem: "ID$:PRINT #2,"":ID$=CHR$(168)
  118. 180  RETURN
  119. 199  '<UNK! {000A}>--- Show a Row of Data ---
  120. 200  PRINT"Row"STR$(I)": ";:FOR J=1 TO M:PRINT X(I,J);:NEXT J:PRINT:RETURN
  121. 209  '<UNK! {000A}>--- Varnames ---
  122. 210  IF VN$="Y" THEN RETURN ELSE FOR L=1 TO M:IF L<10 THEN VN$(L)="Var #"+STR$(L) ELSE VN$(L)="Var #"+MID$(STR$(L),2)
  123. 211  NEXT L:RETURN
  124. 339  '<UNK! {000A}>--- Date ---
  125. 340  DAT$=MID$(DATE$,4,2)+" "+MID$("JanFebMarAprMayJunJulAugSepOctNovDec",-2+3*VAL(LEFT$(DATE$,2)),3)+" "+RIGHT$(DATE$,4):RETURN
  126. 349  '<UNK! {000A}>--- Menu ---
  127. 350  Z$=HD$:CLS:GOSUB 46:PRINT
  128. 351  I=23:J=42:PRINT"1  A'"TAB(I)"7  A * k"TAB(J)"13  A' * B(Inv) * A<UNK! {000A}>2  A * A'"TAB(I)"8  A / k"TAB(J)"14  Normalize<UNK! {000A}>3  A(Inv)"TAB(I)"9  A * B"TAB(J)"15  Sqrt Diag Any Square A"
  129. 352  I=22:PRINT"4  A + B"TAB(I)"10  B * A"TAB(J)"16  Inv Sqrt Diag Any Square A<UNK! {000A}>5  A - B"TAB(I)"11  A * B(Inv)"TAB(J)"17  Eigenstructure Symmetric A<UNK! {000A}>6  B - A"TAB(I)"12  A' * B * A"TAB(J)"18  Eigenstructure B(Inv) * A":PRINT:RETURN
  130. 399  '<UNK! {000A}>--- Start ---
  131. 400  KEY OFF:CLEAR:SCREEN 0,0,0,0:CLS:ON ERROR GOTO 30:DEFINT I-N,Q:I=22
  132. 401  DIM J,K,S1,N,L,X,S2,M,IJ,IC,XX,J1,Y,IR,IC,M1,Z,NM,MA,JJ,IV,AM,LL,F$,J2,NP,MJ,S,BD,LN,D,KK,P,M2,NV,SU,IR,X$,Y$,Z$,TR,KL,ME,M4,EM,Q,QP,R$,EI,MT
  133. 402  DIM CC(I),BB(I),C(I,I),V(I,I),DD(I),E(I),AA(I),VV(I*I),R(I*I),D(I,I),A(I,I),B(I,I),IP(I),P(I),KR(I),KC(I),Y(I),M$(2),VN$(I)
  134. 403  I=I-2:MD=I:MXR=I:MXC=I:MNR=1:MNC=1:NOP=18
  135. 404  F4$="####.##### ":F5$="#####.#### ":F8$="########.# "
  136. 405  D$="Determinant of ":E$="****>  ERROR.   ":L$=", with each row on a new line:":M$(1)="MATRIX A":M$(2)="MATRIX B"
  137. 406  NC$="These matrices are NOT CONFORMABLE.":NS$=" is not square.":P$="PRODUCT, ":EQA$=" = New A =":R$="Row##:   ":U$="Upper Triangular "
  138. 407  '--------
  139. 408  HD$="  M A T R I X    O P E R A T I O N S  ":VER$="(RL,6)"
  140. 409  QB=1:CLOSE:CLS:GOSUB 340:PRINT DAT$;TAB(40-LEN(HD$)\2);:COLOR 0,7:PRINT HD$;:COLOR 7,0:PRINT TAB(73)VER$:LOCATE 4,1,0:K=12
  141. 410  PRINT TAB(K)"KEYTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENCLOSE"
  142. 411  PRINT TAB(K)"OPEN  A menu with 18 OPTIONS will appear whenever needed.  OPEN"
  143. 412  PRINT TAB(K)"OPEN           MATRICES can be up to 20 x 20.              OPEN"
  144. 413  PRINT TAB(K)"OPEN VECTORS are handled as matrices with 1 row or column. OPEN"
  145. 414  PRINT TAB(K)"OPEN   Square symmetric matrices can be UPPER TRIANGULAR.  OPEN"
  146. 415  PRINT TAB(K)"OPEN                                                       OPEN"
  147. 416  PRINT TAB(K)"OPEN  Matrices are to be entered BY ROWS, in Free Format,  OPEN"
  148. 417  PRINT TAB(K)"OPEN          from keyboard, or from datafile.             OPEN"
  149. 418  PRINT TAB(K)"OPEN   Simple data CORRECTIONS can be made in all cases.   OPEN"
  150. 419  PRINT TAB(K)"OPEN                                                       OPEN"
  151. 420  PRINT TAB(K)"OPEN       The 1st matrix entered will be CALLED A.        OPEN"
  152. 421  PRINT TAB(K)"OPEN    If a 2nd matrix is needed, it will be CALLED B.    OPEN"
  153. 422  PRINT TAB(K)"OPEN  Outcome matrices are automatically RE-NAMED as A to  OPEN"
  154. 423  PRINT TAB(K)"OPEN  allow CONTINUED processing, like A x B(Inv) + C - D. OPEN"
  155. 424  PRINT TAB(K)"OPEN                                                       OPEN"
  156. 425  PRINT TAB(K)"OPEN    Printouts &/or filing can be done after displays.  OPEN"
  157. 426  PRINT TAB(K)"OPEN         If trouble, try <Ctrl-Break> & GOTO 9.        OPEN"
  158. 427  PRINT TAB(K)"SCREENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENTHENLOAD":PRINT:LOCATE ,,1:GOSUB 5
  159. 428  '<UNK! {000A}>--- Get Data ---
  160. 429  CLOSE:GOSUB 350:MT=1:ID$=""
  161. 430  DO$="enter "+CHR$(34)+M$(MT)+CHR$(34)+" from Disk":GOSUB 20:IF Z$="N" THEN 434
  162. 431  '<UNK! {000A}>--- Disk Entry ---
  163. 432  QB=3:GOSUB 80:GOSUB 5:GOTO 441
  164. 433  '<UNK! {000A}>--- Keyboard Entry ---
  165. 434  QB=3:QD=0:I=0:M=2:UT=0:PRINT"# of Rows & # of Cols in "M$(MT)" (Free Format, both max "MID$(STR$(MXC),2)")";:INPUT X$:GOSUB 50:N=C(0,1):M=C(0,2):IF N<1 OR N>MXR OR M<1 OR M>MXC THEN BEEP:PRINT"Size out of bounds!":GOTO 434
  166. 435  IF MT=1 THEN NA=N:MA=M ELSE NB=N:MB=M:GOTO 476
  167. 436  I=1:IF N=M THEN PRINT U$"Matrix";:GOSUB 21:IF Z$="Y" THEN 438
  168. 437  UT=0:PRINT"Enter "M$(MT)L$:GOTO 439
  169. 438  UT=1:PRINT"Enter "U$M$(MT)L$
  170. 439  PRINT"Row"I;:INPUT X$:GOSUB 50:IF I<N THEN I=I+1:GOTO 439
  171. 440  '<UNK! {000A}>--- Show Matrix ---
  172. 441  QB=4:GOSUB 160:PRINT #2,
  173. 442  IF QR=0 THEN PRINT #2,M$(MT);" was read as:" ELSE PRINT #2,"Revised ";M$(MT);" is:"
  174. 443  MM=M:L=1:FOR I=1 TO N:PRINT #2,USING R$;I,:LN=0:JJ=9:IF UT>0 THEN MM=M-I+1:IF QD=1 THEN JJ=4
  175. 444  J1=1+LN*(JJ+1):J2=-(MM<J1+JJ)*MM-(MM>=J1+JJ)*(J1+JJ):FOR J=J1 TO J2:PRINT #2,C(I,J);:NEXT J:PRINT #2,:L=L+1:IF L>19 THEN GOSUB 5:L=0
  176. 445  IF MM>J2 THEN LN=LN+1:PRINT #2,SPACE$(9);:GOTO 444
  177. 446  NEXT I:PRINT #2,:IF PR=1 THEN RETURN
  178. 447  '<UNK! {000A}>--- Edit Matrix ---
  179. 448  PRINT"Any data changes";:GOSUB 21:IF Z$="N" THEN 458
  180. 449  IF N=1 THEN I=1:GOTO 451
  181. 450  INPUT"Which row";I:IF I<1 OR I>N THEN PRINT"Silly.":GOTO 450
  182. 451  IF UT>0 THEN L=M-I+1 ELSE L=M
  183. 452  PRINT USING R$;I;:FOR J=1 TO L:PRINT C(I,J);:NEXT J:PRINT
  184. 453  QR=1:IF L=1 THEN INPUT"Right value";C(I,L):GOTO 441
  185. 454  PRINT"Change which variable # (1-"MID$(STR$(L),2)") ";:INPUT Z$:J=VAL(Z$):IF J<1 OR J>L THEN BEEP:PRINT"What ?":GOTO 454
  186. 455  PRINT"Old value =";C(I,J);TAB(28);"New value";:INPUT Z$:FOR K=1 TO LEN(Z$):IF INSTR("-.0123456789",MID$(Z$,K,1)) THEN NEXT K:C(I,J)=VAL(Z$) ELSE PLAY"L16O3CEL4>B":PRINT"That contains a `non-numeric' entry.  Please re-do.":GOTO 455
  187. 456  GOTO 441
  188. 457  '<UNK! {000A}>--- Print Data? ---
  189. 458  DO$="print this data":GOSUB 20:IF Z$="N" THEN 461
  190. 459  GOSUB 162:GOSUB 165:GOSUB 442:GOSUB 160
  191. 460  '<UNK! {000A}>--- Copy C into A or B ---
  192. 461  CLS:QR=0:IF UT>0 THEN 465 ELSE ON MT GOTO 462,463
  193. 462  FOR I=1 TO N:FOR J=1 TO M:A(I,J)=C(I,J):NEXT J:NEXT I:GOTO 470
  194. 463  FOR I=1 TO N:FOR J=1 TO M:B(I,J)=C(I,J):NEXT J:NEXT I:GOTO 481
  195. 464  '     First assign UTM values to their proper subscripts.
  196. 465  FOR I=2 TO M:FOR J=M TO I STEP -1:C(I,J)=C(I,J-I+1):NEXT J:NEXT I:ON MT GOTO 467,468
  197. 466  '     Now expand UTM to full matrices, C and A or B.
  198. 467  FOR I=1 TO M:FOR J=I TO M:C(J,I)=C(I,J):A(I,J)=C(I,J):A(J,I)=A(I,J):NEXT J:NEXT I:GOTO 470
  199. 468  FOR I=1 TO M:FOR J=I TO M:C(J,I)=C(I,J):B(I,J)=C(I,J):B(J,I)=B(I,J):NEXT J:NEXT I:GOTO 481
  200. 469  '<UNK! {000A}>--- Opcodes ---
  201. 470  QB=5:GOSUB 350:GOSUB 160
  202. 471  PRINT USING"Which ! (1-##)";"#";NOP;:INPUT OP:PRINT
  203. 472  ON OP GOTO 490,492,475,473,473,473,504,504,473,473,473,473,473,529,475,475,475,473:GOTO 471
  204. 473  MT=2:GOTO 430
  205. 474  '<UNK! {000A}>--- Test Legality ---
  206. 475  IF NA<>MA THEN BEEP:PRINT E$M$(1)NS$:GOTO 471 ELSE IF OP=3 THEN 494 ELSE ON (OP-14) GOTO 533,538,543
  207. 476  ON OP GOTO 2,2,2,477,477,477,2,2,478,479,478,479,479,2,2,2,2,479
  208. 477  IF NA<>NB OR MA<>MB THEN PRINT E$NC$:GOTO 430 ELSE 480
  209. 478  IF MA<>NB THEN PRINT E$NC$:GOTO 430 ELSE IF OP=9 THEN 480 ELSE IF NB<>MB THEN PRINT E$M$(2)NS$:GOTO 430 ELSE 480
  210. 479  IF NA<>MB THEN PRINT E$NC$:GOTO 430 ELSE IF OP=10 THEN 480 ELSE IF NB<>MB THEN PRINT E$M$(2)NS$:GOTO 430
  211. 480  IF QD=0 THEN 436 ELSE ON MT GOTO 441,88 'Continue data entry from disk if B.
  212. 481  GOSUB 43:ON OP GOTO 2,2,2,498,500,502,2,2,510,514,516,520,525,2,2,2,2,552
  213. 482  '<UNK! {000A}>--- Finale ---
  214. 483  DO$="file this matrix":GOSUB 20:IF Z$="N" THEN 487
  215. 484  INPUT "How many vectors of measurements contributed to this matrix";UT:IF UT<1 THEN BEEP:GOTO 484
  216. 485  IF VN$="" THEN VN$="N"
  217. 486  UT=-UT:GOSUB 140  'hold # persons as -ut, then file full matrix (not UTM).
  218. 487  PRINT:PRINT"Another operation on ";:COLOR 23,0:PRINT"this";:COLOR 7,0:PRINT" matrix";:GOSUB 21:IF Z$="Y" THEN 470 ELSE IF HEAD=1 THEN HEAD=0:LPRINT STRING$(79,61)STRING$(4,10)
  219. 488  CLOSE:BEEP:PRINT:PRINT"===> ANOTHER MATRIX PROBLEM";:GOSUB 21:IF Z$="Y" THEN CLS:GOTO 429 ELSE 30
  220. 489  ' ------------------- Operations -------------------<UNK! {000A}><UNK! {000A}>--1  A'
  221. 490  FOR I=1 TO N:FOR J=1 TO M:C(J,I)=A(I,J):NEXT J:NEXT I:SWAP NA,MA:N=NA:M=MA:FOR I=1 TO N:FOR J=1 TO M:A(I,J)=C(I,J):NEXT J:NEXT I
  222. 491  PRINT #2,"TRANSPOSE, A'"EQA$:GOSUB 564:IF PR=1 THEN 491 ELSE 483'<UNK! {000A}><UNK! {000A}>--2  A * A'
  223. 492  GOSUB 43:FOR I=1 TO N:FOR J=1 TO N:C(I,J)=0:FOR K=1 TO M:C(I,J)=C(I,J)+A(I,K)*A(J,K):NEXT K:NEXT J:NEXT I:MA=NA:M=MA:FOR I=1 TO N:FOR J=1 TO M:A(I,J)=C(I,J):NEXT J:NEXT I
  224. 493  PRINT #2,P$"A * A'"EQA$:GOSUB 564:IF PR=1 THEN 493 ELSE 483'<UNK! {000A}><UNK! {000A}>--3  A(Inv)
  225. 494  GOSUB 43:FOR I=1 TO M:FOR J=1 TO M:D(I,J)=A(I,J):NEXT J:NEXT I:GOSUB 574
  226. 495  PRINT #2,D$M$(MT)" = ";D:PRINT #2,:IF D=0 THEN GOSUB 5:GOTO 483
  227. 496  IF PR=0 THEN FOR I=1 TO M:FOR J=1 TO M:A(I,J)=C(I,J):NEXT J:NEXT I
  228. 497  PRINT #2,"INVERSE MATRIX, A(Inv)"EQA$:GOSUB 564:GOTO 594'<UNK! {000A}><UNK! {000A}>--4  A+B
  229. 498  FOR I=1 TO N:FOR J=1 TO M:A(I,J)=A(I,J)+B(I,J):NEXT J:NEXT I
  230. 499  PRINT #2,"SUM, A + B"EQA$:GOSUB 564:IF PR=1 THEN 499 ELSE 483'<UNK! {000A}><UNK! {000A}>--5  A-B
  231. 500  FOR I=1 TO N:FOR J=1 TO M:A(I,J)=A(I,J)-B(I,J):NEXT J:NEXT I
  232. 501  PRINT #2,"DIFFERENCE, A - B"EQA$:GOSUB 564:IF PR=1 THEN 501 ELSE 483'<UNK! {000A}><UNK! {000A}>--6  B-A
  233. 502  FOR I=1 TO N:FOR J=1 TO M:A(I,J)=B(I,J)-A(I,J):NEXT J:NEXT I
  234. 503  PRINT #2,"REVERSE DIFFERENCE, B - A"EQA$:GOSUB 564:IF PR=1 THEN 503 ELSE 483'<UNK! {000A}><UNK! {000A}>--7  A * k
  235. 504  INPUT"Constant (any number, whole or decimal)";C:IF OP=8 THEN 508
  236. 505  FOR I=1 TO N:FOR J=1 TO M:A(I,J)=A(I,J)*C:NEXT J:NEXT I
  237. 506  GOSUB 507:PRINT #2,"* k"EQA$:GOSUB 564:IF PR=1 THEN 506 ELSE 483
  238. 507  PRINT #2,:PRINT #2,"With k ="STR$(C)", SCALED MATRIX A ";:RETURN'<UNK! {000A}><UNK! {000A}>--8  A / k
  239. 508  FOR I=1 TO N:FOR J=1 TO M:A(I,J)=A(I,J)/C:NEXT J:NEXT I
  240. 509  GOSUB 507:PRINT #2,"/ k"EQA$:GOSUB 564:IF PR=1 THEN 509 ELSE 483'<UNK! {000A}><UNK! {000A}>--9  A * B
  241. 510  FOR I=1 TO NA:FOR J=1 TO MB:C(I,J)=0:FOR K=1 TO MA:C(I,J)=C(I,J)+A(I,K)*B(K,J):NEXT K:NEXT J:NEXT I:N=NA:MA=MB:M=MA:FOR I=1 TO N:FOR J=1 TO M:A(I,J)=C(I,J):NEXT J:NEXT I
  242. 511  IF OP=11 THEN 513
  243. 512  PRINT #2,P$"A * B"EQA$:GOSUB 564:IF PR=1 THEN 512 ELSE 483
  244. 513  PRINT #2,P$"A * B(Inv)"EQA$:GOSUB 564:IF PR=1 THEN 513 ELSE 483'<UNK! {000A}><UNK! {000A}>--10  B * A
  245. 514  FOR I=1 TO NB:FOR J=1 TO MA:C(I,J)=0:FOR K=1 TO MB:C(I,J)=C(I,J)+B(I,K)*A(K,J):NEXT K:NEXT J:NEXT I:NA=NB:N=NA:M=MA:FOR I=1 TO N:FOR J=1 TO M:A(I,J)=C(I,J):NEXT J:NEXT I
  246. 515  PRINT"REVERSE "P$"B * A"EQA$:GOSUB 564:IF PR=1 THEN 513 ELSE 483'<UNK! {000A}><UNK! {000A}>--11  A * B(Inv)
  247. 516  M=NB:FOR I=1 TO M:FOR J=1 TO M:D(I,J)=B(I,J):NEXT J:NEXT I:GOSUB 574:GOSUB 160
  248. 517  PRINT #2,D$M$(MT)" = "D:PRINT #2,:IF D=0 THEN GOSUB 5:GOTO 483
  249. 518  IF PR=0 THEN FOR I=1 TO M:FOR J=1 TO M:A(I,J)=C(I,J):NEXT J:NEXT I
  250. 519  PRINT #2,P$"A * B(Inv)"EQA$:GOSUB 564:GOTO 594'<UNK! {000A}><UNK! {000A}>--12  A' * B * A
  251. 520  FOR I=1 TO MA:FOR J=1 TO MB:C(I,J)=0:FOR K=1 TO NA:C(I,J)=C(I,J)+A(K,I)*B(K,J):NEXT K:NEXT J:NEXT I
  252. 521  FOR I=1 TO MA:FOR J=1 TO MA:B(I,J)=0:FOR K=1 TO NA:B(I,J)=B(I,J)+C(I,K)*A(K,J):NEXT K:NEXT J:NEXT I
  253. 522  NA=MA:N=NA:M=MA:FOR I=1 TO N:FOR J=1 TO M:A(I,J)=B(I,J):NEXT J:NEXT I:GOSUB 160
  254. 523  IF OP=13 THEN 528
  255. 524  PRINT #2,"QUADRATIC "P$"A' * B * A"EQA$:GOSUB 564:IF PR=1 THEN 524 ELSE 483'<UNK! {000A}><UNK! {000A}>--13  A' * B(Inv) * A
  256. 525  M=NB:FOR I=1 TO M:FOR J=1 TO M:C(I,J)=B(I,J):D(I,J)=B(I,J):NEXT J:NEXT I:GOSUB 574:GOSUB 160
  257. 526  PRINT #2,D$M$(MT)" = "D:PRINT #2,:IF D=0 THEN GOSUB 5:GOTO 483
  258. 527  IF PR=0 THEN FOR I=1 TO M:FOR J=1 TO M:B(I,J)=C(I,J):NEXT J:NEXT I:GOTO 520
  259. 528  PRINT #2,"QUADRATIC "P$"A' * B(Inv) * A"EQA$:GOSUB 564:GOTO 594'<UNK! {000A}><UNK! {000A}>--14  Normalize
  260. 529  NA=N:MA=M:GOSUB 160:IF N>1 THEN 531 ELSE S=0:FOR J=1 TO M:S=S+A(1,J)*A(1,J):NEXT J:S=SQR(S):FOR J=1 TO M:A(1,J)=A(1,J)/S:NEXT J
  261. 530  PRINT #2,"NORMALIZED VECTOR"EQA$:GOSUB 564:IF PR=1 THEN 530 ELSE 483
  262. 531  FOR J=1 TO M:E(J)=0:FOR I=1 TO N:E(J)=E(J)+A(I,J)*A(I,J):NEXT I:NEXT J:FOR J=1 TO M:E(J)=SQR(E(J)):FOR I=1 TO N:A(I,J)=A(I,J)/E(J):NEXT I:NEXT J:IF M=1 THEN 530
  263. 532  PRINT #2,"NORMALIZED MATRIX"EQA$:GOSUB 564:IF PR=1 THEN 532 ELSE 483'<UNK! {000A}><UNK! {000A}>--15  Sqrt Diag A
  264. 533  FOR I=1 TO N:FOR J=1 TO M:IF I<>J THEN A(I,J)=0:GOTO 535 'diagonalize the matrix.
  265. 534  IF A(I,I)>=0 THEN A(I,I)=SQR(A(I,I)) ELSE FLAG=I:I=N:J=M
  266. 535  NEXT J:NEXT I:IF FLAG>0 THEN 537 ELSE GOSUB 160
  267. 536  PRINT #2,"SQRT OF DIAGONAL A"EQA$:GOSUB 564:IF PR=1 THEN 536 ELSE 483
  268. 537  BEEP:PRINT"Can't SQRT, because diagonal element"FLAG"is negative.":FLAG=0:GOTO 488'<UNK! {000A}><UNK! {000A}>--16  Inv Sqrt Diag A
  269. 538  FOR I=1 TO N:FOR J=1 TO M:IF I<>J THEN A(I,J)=0:GOTO 541 'diagonalize the matrix.
  270. 539  IF A(I,I)<0 THEN FLAG=I:I=N:J=M:GOTO 541 ELSE IF A(I,I)<1E-10 THEN A(I,I)=0
  271. 540  A(I,I)=1/SQR(A(I,I))
  272. 541  NEXT J:NEXT I:IF FLAG>0 THEN 537 ELSE GOSUB 160
  273. 542  PRINT #2,"INVERSE SQRT OF DIAGONAL A"EQA$:GOSUB 564:IF PR=1 THEN 542 ELSE 483'<UNK! {000A}><UNK! {000A}>--17  Eigen Symmetric A
  274. 543  IF M=2 THEN NV=2 ELSE EMIN=M:PRINT"How many Eigenvalues required (2-"MID$(STR$(M),2)"), or Min Value (less than 2)";:INPUT EMIN:IF EMIN>M THEN 543 ELSE IF EMIN<2 THEN NV=M ELSE NV=EMIN
  275. 544  GOSUB 43:TR=0:FOR I=1 TO M:TR=TR+A(I,I):NEXT I:GOSUB 599:CLS:GOSUB 160
  276. 545  PRINT #2,D$M$(MT)" = ";:IF NV<M THEN PRINT #2,"?" ELSE D=1:FOR I=1 TO M:D=D*E(I):NEXT I:PRINT #2,D
  277. 546  PRINT #2,"EIGENVALUES ";:S=0:Q=0:IF OP=18 THEN PRINT #2,"=":T$=LEFT$(P$,7):GOTO 548
  278. 547  T$="A":IF M=2 THEN PRINT #2,"=" ELSE IF EMIN<2 THEN PRINT #2,"(No. determined by MIN VALUE ="STR$(EMIN)"):" ELSE PRINT #2,"("MID$(STR$(NV),2)" values, as requested) ="
  279. 548  IF PR=0 THEN GOSUB 160
  280. 549  GOSUB 559:PRINT #2,"SUM OF EIGENVALUES = "S"    TRACE OF "T$" = "TR:PRINT #2,:GOSUB 5
  281. 550  PRINT #2,"NORMALIZED EIGENVECTORS (columnwise)"EQA$:N=NA:M=NV:MA=M:Q=1:GOSUB 559:FOR I=1 TO N:FOR J=1 TO M:A(I,J)=V(I,J):NEXT J:NEXT I
  282. 551  GOSUB 570:IF PR=1 THEN 545 ELSE 483'<UNK! {000A}><UNK! {000A}>--18  Eigen B(Inv) * A
  283. 552  M=MA:FOR I=1 TO NA:FOR J=1 TO M:C(I,J)=A(I,J):NEXT J:NEXT I:GOSUB 574:DA=D:FOR I=1 TO NB:FOR J=1 TO MB:C(I,J)=B(I,J):NEXT J:NEXT I:GOSUB 574:DB=D
  284. 553  FOR I=1 TO NA:FOR J=1 TO M:V(I,J)=0:FOR K=1 TO M:V(I,J)=V(I,J)+C(I,K)*A(K,J):NEXT K:NEXT J:NEXT I
  285. 554  TR=0:FOR I=1 TO NA:TR=TR+V(I,I):NEXT I
  286. 555  GOSUB 662:CLS:GOSUB 160
  287. 556  PRINT #2,"EIGENSTRUCTURE OF B(Inv) * A":PRINT #2,:PRINT #2,D$M$(1)" = "DA:PRINT #2,D$M$(2)" = "DB:GOTO 546
  288. 557  '<UNK! {000A}>          -------------- More Subroutines ---------------
  289. 558  '<UNK! {000A}>--- Print Eigenstructure ---
  290. 559  FOR K=1 TO NV STEP 5:IF NV<K+5 THEN L=NV ELSE L=K+4
  291. 560  FOR I=K TO L:PRINT #2,USING"     ##     ";I;:NEXT I:PRINT #2,:IF Q=0 THEN FOR J=K TO L:PRINT #2,USING" "+F4$;E(J);:S=S+E(J):NEXT J:PRINT #2,:GOTO 562
  292. 561  FOR I=1 TO N:FOR J=K TO L:PRINT #2,USING" "+F4$;V(I,J);:NEXT J:PRINT #2,:NEXT I:PRINT #2,
  293. 562  NEXT K:RETURN
  294. 563  '<UNK! {000A}>--- Show/Print Outcome Matrix,  A(I,J) ---
  295. 564  IF ABS(A(1,1))<10 THEN F$=F4$ ELSE IF ABS(A(1,1))<1000 THEN F$=F5$ ELSE F$=F8$
  296. 565  L=1:FOR I=1 TO NA:PRINT #2,USING R$;I;:LN=0
  297. 566  J1=1+LN*5:J2=-(MA<J1+4)*MA-(MA>=J1+4)*(J1+4):IF J1>5 THEN PRINT #2,SPACE$(9);
  298. 567  FOR J=J1 TO J2:PRINT #2,USING F$;A(I,J);:NEXT J:PRINT #2,
  299. 568  GOSUB 572:IF MA>J2 THEN LN=LN+1:GOTO 566
  300. 569  NEXT I:PRINT #2,:IF OP=3 OR OP=11 OR OP=13 THEN 571
  301. 570  IF PR=1 THEN GOSUB 160 ELSE GOSUB 161:IF PR=1 THEN GOSUB 165
  302. 571  RETURN'<UNK! {000A}>--- Line Count ---
  303. 572  IF PR=1 THEN RETURN ELSE IF L<20 THEN L=L+1:RETURN ELSE GOSUB 5:L=1:RETURN
  304. 573  '<UNK! {000A}>--- Inverse of C(M,M)   (Gauss-Jordan Method) ---<UNK! {000A}>    Uses I-M,D,Q,Z,FLAG,IR,IC,IP(M),KR(M),KC(M),P(M).
  305. 574  D=1:FLAG=0:FOR I=1 TO M:IP(I)=0:NEXT I:FOR I=1 TO M:Z=0:FOR J=1 TO M:IF IP(J)=1 THEN 579
  306. 575  FOR K=1 TO M:ON SGN(IP(K)-1)+2 GOTO 577,578,576
  307. 576  I=M:J=M:K=M:FLAG=1:GOTO 578
  308. 577  IF Z < ABS(C(J,K)) THEN IR=J:IC=K:Z=C(J,K)
  309. 578  NEXT K
  310. 579  NEXT J:IF FLAG=1 THEN 585
  311. 580  IP(IC)=IP(IC)+1:IF IR<>IC THEN D=-D:FOR L=1 TO M:Z=C(IR,L):C(IR,L)=C(IC,L):C(IC,L)=Z:NEXT L
  312. 581  KR(I)=IR:KC(I)=IC:P(I)=C(IC,IC):D=D*P(I):IF ABS(P(I))<=9.99E-07 THEN D=0:I=M:GOTO 585
  313. 582  C(IC,IC)=1:FOR L=1 TO M:C(IC,L)=C(IC,L)/P(I):NEXT L
  314. 583  FOR K=1 TO M:IF K<>IC THEN Z=C(K,IC):C(K,IC)=0:FOR L=1 TO M:C(K,L)=C(K,L)-C(IC,L)*Z:NEXT L
  315. 584  NEXT K
  316. 585  NEXT I:IF FLAG=1 OR D=0 THEN 589
  317. 586  FOR I=1 TO M:Q=M-I+1:IF KR(Q)=KC(Q) THEN 588 ELSE K=KR(Q):L=KC(Q)
  318. 587  FOR J=1 TO M:Z=C(J,K):C(J,K)=C(J,L):C(J,L)=Z:NEXT J
  319. 588  NEXT I
  320. 589  '      Verify inverse
  321. 590  FOR I=1 TO M:FOR J=1 TO M:V(I,J)=0:FOR K=1 TO M:V(I,J)=V(I,J)+D(I,K)*C(K,J):NEXT K:NEXT J:NEXT I:FOR I=1 TO M:FOR J=1 TO M:V(I,J)=ABS(V(I,J)):NEXT J:NEXT I
  322. 591  IF ABS(D) < 9.99E-07 THEN D=0
  323. 592  RETURN
  324. 593  '<UNK! {000A}>---  Print C*C(Inv)  ---
  325. 594  IF OP=3 THEN X$="A*A(Inv)" ELSE X$="B*B(Inv)":IF OP=13 THEN M=NB
  326. 595  GOSUB 572:PRINT #2,"Above results are Ok if the following "X$" product is an Identity Matrix:"
  327. 596  FOR I=1 TO M:FOR J=1 TO M:PRINT #2,USING"###.##";V(I,J);:NEXT J:PRINT #2,:GOSUB 572:NEXT I:IF PR=1 THEN PRINT #2,
  328. 597  GOSUB 570:IF PR=1 THEN 495 ELSE 483
  329. 598  '<UNK! {000A}>---  Sub HOW - (Householder, Ortega & Wilkinson, after W W Cooley & P R Lohnes)<UNK! {000A}>        Input = C(I,J), M, MD, NV, EMIN.   Output = V(I,J), E(I), NV.
  330. 599  FOR I=1 TO M:FOR J=1 TO M:R(I+(J-1)*MD)=C(I,J):NEXT J:NEXT I  'Equates uni-dimensional R(I) with C(I,J).
  331. 600  IF EMIN<2 THEN EIG=EMIN ELSE EIG= -1E+20
  332. 601  '<UNK! {000A}>      Tri-diagonalize.
  333. 602  M1=M-1:M2=M1*MD+M:M3=M2-MD:M4=MD+1:L=0:FOR I=1 TO M2 STEP M4:L=L+1:AA(L)=R(I):NEXT I:BB(1)=0:IF M=2 THEN 614 ELSE KK=0
  334. 603  FOR K=2 TO M1:KL=KK+K:KU=KK+M:KJ=K+1:SUM=0
  335. 604  FOR J=KL TO KU:SUM=SUM+R(J)*R(J):NEXT J:S=SQR(SUM):Z=R(KL):IF Z>0 THEN BB(K)=-S ELSE BB(K)=S
  336. 605  S=1/S:CC(K)=SQR(ABS(Z)*S+1):R(KL)=CC(K):X=ABS(S/CC(K)):IF Z<0 THEN X=-X
  337. 606  FOR I=KJ TO M:JJ=I+KK:CC(I)=X*R(JJ):R(JJ)=CC(I):NEXT I
  338. 607  FOR J=K TO M:JJ=J+1:DD(J)=0:L=KK+J
  339. 608  FOR I=K TO J:L=L+MD:DD(J)=DD(J)+R(L)*CC(I):NEXT I:IF J=M THEN 610
  340. 609  FOR I=JJ TO M:L=L+1:DD(J)=DD(J)+R(L)*CC(I):NEXT I
  341. 610  NEXT J:X=0:FOR J=K TO M:X=X+CC(J)*DD(J):NEXT J:X=0.5*X:FOR I=K TO M:DD(I)=X*CC(I)-DD(I):NEXT I:LL=KK:KK=KK+MD
  342. 611  FOR I=K TO M:LL=LL+MD:FOR J=I TO M:L=LL+J:R(L)=R(L)+DD(I)*CC(J)+DD(J)*CC(I):NEXT J:NEXT I
  343. 612  NEXT K:L=1:FOR I=1 TO M:X=AA(I):AA(I)=R(L):R(L)=X:L=L+M4:NEXT I
  344. 613  '<UNK! {000A}>       Eigenvalues.
  345. 614  EM6=9.99E-07:EM9=0:EM15=0:E15=1E+15:BB(M)=R(M3):BD=ABS(AA(1))
  346. 615  FOR I=2 TO M:XX=ABS(AA(I))+BB(I)*BB(I):IF BD<XX THEN BD=XX
  347. 616  NEXT I:BD=BD+1:FOR I=1 TO M:AA(I)=AA(I)/BD:BB(I)=BB(I)/BD:DD(I)=1:E(I)=-1:NEXT I
  348. 617  FOR K=1 TO M
  349. 618  XX=ABS(DD(K)):IF ABS(E(K))>XX THEN XX=ABS(E(K))
  350. 619  IF EM9 > XX THEN XX=EM9
  351. 620  IF((DD(K)-E(K))/XX) < EM6 THEN 633
  352. 621  X=(DD(K)+E(K))*0.5:S2=1:CC(1)=AA(1)-X:IF CC(1)<0 THEN S1=-1:N=0 ELSE S1=1:N=1
  353. 622  FOR I=2 TO M
  354. 623  IF BB(I)=0 THEN CC(I)=(AA(I)-X)*SGN(S1):GOTO 627
  355. 624  IF BB(I-1)=0 THEN CC(I)=(AA(I)-X)*CC(I-1)-BB(I)*BB(I)*SGN(S2):GOTO 627
  356. 625  IF ABS(CC(I-1))+ABS(CC(I-2)) < EM15 THEN CC(I-1)=CC(I-1)*E15:CC(I-2)=CC(I-2)*E15
  357. 626  CC(I)=(AA(I)-X)*CC(I-1)-BB(I)*BB(I)*CC(I-2)
  358. 627  S2=S1:IF CC(I) THEN S1=ABS(S1)*SGN(CC(I)):IF S1=-S2 THEN 629
  359. 628  N=N+1
  360. 629  NEXT I:N=M-N:IF N>=K THEN FOR J=K TO N:DD(J)=X:NEXT J
  361. 630  N=N+1:IF M<N THEN 618
  362. 631  FOR J=N TO M:IF X<=E(J) THEN J=M ELSE E(J)=X
  363. 632  NEXT J:GOTO 618
  364. 633  NEXT K
  365. 634  FOR I=1 TO M:AA(I)=AA(I)*BD:BB(I)=BB(I)*BD:CC(I)=(DD(I)+E(I))*BD*0.5:NEXT I:M1=M:K=1
  366. 635  I=1
  367. 636  FLAG=0:FOR J=1 TO M1:IF I-J THEN IF CC(I)-CC(J)>0 THEN I=J:J=M1:FLAG=1
  368. 637  NEXT J:IF FLAG=1 THEN 636
  369. 638  E(K)=CC(I):K=K+1:M1=M1-1:IF I-M1<1 THEN FOR M2=1 TO M1:CC(M2)=CC(M2+1):NEXT M2
  370. 639  IF M1>1 THEN 635 ELSE E(K)=CC(1):IF NV>=0 THEN FOR I=1 TO M:CC(I)=E(I):NEXT I:J=M:FOR I=1 TO M:E(I)=CC(J):IF E(I)>=EIG THEN J=J-1 ELSE NV=I-1:I=M
  371. 640  NEXT I
  372. 641  '<UNK! {000A}>       Eigenvectors.
  373. 642  EM10=0:IF NV=0 THEN 658 ELSE KX=ABS(NV):J=1
  374. 643  FOR IV=1 TO KX:X=AA(1)-E(IV):Y=BB(2):M1=M-1
  375. 644  FOR I=1 TO M1:IJ=I+J-1:ON SGN(ABS(X)-ABS(BB(I+1)))+2 GOTO 645,646,647
  376. 645  CC(I)=BB(I+1):DD(I)=AA(I+1)-E(IV):VV(IJ)=BB(I+2):Z=-X/CC(I):X=Z*DD(I)+Y:IF M1=I THEN 648 ELSE Y=Z*VV(IJ):GOTO 648
  377. 646  IF X=0 THEN X=EM10
  378. 647  CC(I)=X:DD(I)=Y:VV(IJ)=0:X=AA(I+1)-(BB(I+1)/X*Y+E(IV)):Y=BB(I+2)
  379. 648  NEXT I:MJ=M+J-1:IF X=0 THEN 651 ELSE VV(MJ)=1/X
  380. 649  I=M1:IJ=I+J-1:VV(IJ)=(1-DD(I)*VV(MJ))/CC(I):X=VV(MJ)*VV(MJ)+VV(IJ)*VV(IJ)
  381. 650  I=I-1:IJ=I+J-1:IF I=0 THEN 652 ELSE VV(IJ)=(1-(DD(I)*VV(IJ+1)+VV(IJ)*VV(IJ+2)))/CC(I):X=X+VV(IJ)*VV(IJ):GOTO 650
  382. 651  VV(MJ)=1E+10:GOTO 649
  383. 652  X=SQR(X):FOR I=1 TO M:IJ=I+J-1:VV(IJ)=VV(IJ)/X:NEXT I:J1=M1*MD-MD:K=M:GOTO 654
  384. 653  K=K-1:J1=J1-MD:Y=0:FOR I=K TO M:IJ=I+J-1:L=J1+I:Y=Y+VV(IJ)*R(L):NEXT I:FOR I=K TO M:IJ=I+J-1:L=J1+I:VV(IJ)=VV(IJ)-Y*R(L):NEXT I
  385. 654  IF J1 THEN 653 ELSE NPLUS=0:NMIN=0:FOR I=1 TO M:IJ=I+J-1:IF VV(IJ)<0 THEN NMIN=NMIN+1 ELSE NPLUS=NPLUS+1
  386. 655  NEXT I:IF NPLUS<NMIN THEN FOR I=1 TO M:IJ=I+J-1:VV(IJ)=-VV(IJ):NEXT I
  387. 656  J=J+MD:NEXT IV
  388. 657  '<UNK! {000A}>      Restore Matrices.
  389. 658  ME=MD+1:JJ=ME:M1=M*MD:FOR I=2 TO M1 STEP ME:K=I:FOR J=JJ TO M1 STEP MD:R(K)=R(J):K=K+1:NEXT J:JJ=JJ+ME:NEXT I
  390. 659  FOR I=1 TO M:FOR J=1 TO NV:V(I,J)=VV(I+(J-1)*MD):NEXT J:NEXT I:FOR J=1 TO NV:IF V(1,J)<0 THEN FOR I=1 TO M:V(I,J)=-V(I,J):NEXT I
  391. 660  NEXT J:RETURN
  392. 661  '<UNK! {000A}>---   Sub DIRNM    (After P R Lohnes) ---
  393. 662  FOR I=1 TO M:FOR J=1 TO M:C(I,J)=B(I,J):NEXT J:NEXT I:NV=M:EMIN=M
  394. 663  GOSUB 599:IF NV<M THEN 672
  395. 664  FOR I=1 TO M:D=SQR(ABS(E(I))):FOR J=1 TO M:B(J,I)=V(J,I)*D:NEXT J:NEXT I
  396. 665  FOR I=1 TO M:E(I)=1/SQR(ABS(E(I))):NEXT I
  397. 666  FOR I=1 TO M:FOR J=1 TO M:B(I,J)=V(I,J)*E(J):NEXT J:NEXT I
  398. 667  FOR I=1 TO M:FOR J=1 TO M:V(I,J)=0:FOR K=1 TO M:V(I,J)=V(I,J)+B(K,I)*A(K,J):NEXT K:NEXT J:NEXT I
  399. 668  FOR I=1 TO M:FOR J=1 TO M:C(I,J)=0:FOR K=1 TO M:C(I,J)=C(I,J)+V(I,K)*B(K,J):NEXT K:NEXT J:NEXT I
  400. 669  GOSUB 599:IF NV<M THEN 672
  401. 670  FOR I=1 TO M:FOR J=1 TO M:A(I,J)=0:FOR K=1 TO M:A(I,J)=A(I,J)+B(I,K)*V(K,J):NEXT K:NEXT J:NEXT I
  402. 671  FOR I=1 TO M:S=0:FOR J=1 TO M:S=S+A(J,I)*A(J,I):NEXT J:FOR J=1 TO M:V(J,I)=A(J,I)/SQR(S):NEXT J:NEXT I
  403. 672  RETURN
  404. 673  'end
  405.