DIM A(120)
DIM B(120)
DIM A1(120)
DIM U(15)
DIM Z(15)

CLS

PRINT "INVERSE ITERATION WITH SHIFTS"
PRINT
PRINT
PRINT
PRINT "TYPE(1): DIAGONAL MATRIX"
PRINT "TYPE(2): FULL SYMMETRICAL MATRIX"
PRINT
PRINT

REM ==INPUT DATA==
PRINT
E5 = 1E-08
PRINT "ENTER TYPE OF PROBLEM "
PRINT
INPUT "TYPE= "; I9
PRINT
PRINT
PRINT "DIMENSION MATRIX A"
PRINT
INPUT "N= "; N
J0 = 0
PRINT
PRINT "ENTER MATRIX A UPPER TRIANGLE"
FOR J = 1 TO N
PRINT
PRINT "UPPER PART COLUMN"; J
PRINT
FOR I = 1 TO J
I1 = I + J0
PRINT "A("; I; ","; J; ") = ";
INPUT A1(I1)
NEXT I
PRINT
J0 = J0 + J
NEXT J
PRINT
PRINT
REM ==INPUT MATRIX OR VECTOR B==
IF I9 = 1 THEN GOSUB 2450
IF I9 = 2 THEN GOSUB 2610
PRINT
PRINT
PRINT "ENTER TOLERANCE "
PRINT
INPUT "TOL = "; T0
I8 = 0
540 :
PRINT
PRINT
PRINT "MAXIMUM NUMBER OF ITERATIONS"
PRINT
INPUT "NMAX = "; N8
IF I8 = 2 THEN 970
PRINT
600 :
PRINT "INPUT DESIRED SHIFT"
PRINT
INPUT "SHIFT = "; W0
PRINT
PRINT "-RUNNING-"
PRINT
REM ==INITIALIZE MATRIX A==
JJ = 0
FOR J = 1 TO N
FOR I = 1 TO J
IJ = JJ + I
A(IJ) = A1(IJ)
NEXT I
JJ = JJ + J
NEXT J
REM ==COMPUTE A - SHIFT B==
IF I9 = 1 THEN GOSUB 2550
IF I9 = 2 THEN GOSUB 2780
GOSUB 1740
REM ==CHECK POLARITY==
IF I8 = 1 THEN 600
M0 = 0
II = 0
REM ==CHECK STURM SEQUENCE==
FOR I = 1 TO N
II = II + I
IF A(II) > 0 THEN 820
M0 = M0 + 1
820 :
NEXT I
PRINT
PRINT "NUMBER OF EIGENVALUES"
PRINT
PRINT "SMALLER THAN SHIFT = "; M0
PRINT
PRINT "DETERMINANT = "; D8
900 :
PRINT
PRINT "WANT SHIFT CHANGED (Y/N)";
INPUT W$
PRINT
PRINT "-RUNNING-"
PRINT
IF W$ = "Y" THEN 600
IF I8 = 2 THEN 540
REM ==PERFORM INVERSE ITERATIONS==
970 :
GOSUB 1250
REM ==CHECK FOR CONVERGENCY==
IF I8 = 2 THEN 900
REM ==UNDO SHIFT==
E0 = E0 + W0
GOSUB 1060
PRINT
PRINT "WANT ANOTHER EIGENSOLUTION (Y/N) ";
INPUT W$
IF W$ = "Y" GOTO 540
END
1060 :
REM OUTPUT MODULE
PRINT
PRINT "SOLUTION"
PRINT "**********"
PRINT
PRINT "EIGENVALUE"
PRINT
PRINT " EIG = "; E0
PRINT
PRINT "EIGENVECTOR"
PRINT
FOR I = 1 TO N
PRINT "U("; I; ")= "; U(I)
NEXT I
PRINT
PRINT "NUM. ITERATIONS"
PRINT
PRINT " N = "; N0
PRINT
RETURN
1250 :
REM ==INVERSE ITERATION SUBROUTINE==
I8 = 0
REM ==INITIAL EIGENVECTOR==
II = 0
FOR I = 1 TO N
II = II + I
B1 = B(I)
IF I9 = 1 THEN 1320
B1 = B(II)
1320 :
U(I) = A(II) / B1
NEXT I
E1 = 1E+32
N0 = 0
R1 = 0
GOTO 1500
1380 :
REM ==SET ITERATION COUNTER==
N0 = N0 + 1
FOR I = 1 TO N
U(I) = Z(I)
NEXT I
REM ==FORWARD SOLUTION==
GOSUB 2180
REM ==COMPUTE RALEIGH QUOTIENT NUMERATOR==
R1 = 0
II = 0
FOR I = 1 TO N
II = II + I
R1 = R1 + U(I) * U(I) / A(II)
NEXT I
REM ==BACK SUBSTITUTION==
GOSUB 2280
REM ==COMPUTE B-NORM OF U==
1500 :
IF I9 = 1 THEN GOSUB 2870
IF I9 = 2 THEN GOSUB 2910
V2 = 0
FOR I = 1 TO N
V2 = V2 + Z(I) * U(I)
NEXT I
V3 = SQR(V2)
REM==SCALE Z BY NORM OF U==
FOR I = 1 TO N
Z(I) = Z(I) / V3
NEXT I
REM ==EIGENVALUE ESTIMATION==
E0 = R1 / V2
H1 = ABS(E1)
H0 = ABS(E0)
REM ==TEST OF CONVERGENCE==
IF ABS(H1 - H0) < (T0 * H0) THEN 1700
E1 = E0
IF N0 < N8 THEN 1380
REM ==ERROR MESSAGES==
I8 = 2
PRINT
PRINT " **WARNING** CONVERGENCE NOT"
PRINT "ATTAINED. CHANGE NMAX OR SHIFT"
RETURN
1700 :
REM ==FINAL EIGENVECTOR==
FOR I = 1 TO N
U(I) = U(I) / V3
NEXT I
RETURN
REM ==FACTORIZATION OF A ==
1740 :
I8 = 0
J0 = 1
D8 = A(J0)
REM: LOOP OVER COLUMNS
FOR J = 2 TO N
J2 = J - 1
IF J = 2 THEN 1910
I0 = 1
REM: LOOP OVER ROWS
FOR I = 2 TO J2
I2 = I - 1
I1 = J0 + I
REM REDUCE COLUMN 11
FOR K = 1 TO I2
K1 = I0 + K
K2 = J0 + K
A(I1) = A(I1) - A(K1) * A(K2)
NEXT K
I0 = I0 + I
NEXT I
REM ==CHECK FOR ZERO PIVOT==
1910 :
IF ABS(A(J0)) < E5 THEN 2120
REM ==REDUCE DIAGONAL ELEMENT AND COLUMNJ==
J3 = J0 + J
I3 = 0
D1 = A(J3)
FOR I = 1 TO J2
I3 = I3 + I
I1 = J0 + I
U = A(I1) / A(I3)
D1 = D1 - U * A(I1)
A(I1) = U
NEXT I
A(J3) = D1
D8 = D8 * D1
J0 = J0 + J
NEXT J
REM ==CHECK FOR ZERO PIVOT==
IF ABS(A(J0)) > E5 THEN 2160
REM ERROR MESSAGES
PRINT
PRINT "**WARNING** ZERO DIAGONAL"
PRINT "ELEMENT. CHANGE SHIFT"
PRINT
RETURN
2120 :
PRINT "UNSUCCESSFUL FACTORIZATION"
PRINT "CHOOSE NEW SHIFT"
I8 = 1
PRINT
2160 :
RETURN
2180 :
REM ==FORWARD SOLUTION==
J0 = 1
FOR J = 2 TO N
J2 = J - 1
FOR K = 1 TO J2
J4 = J0 + K
U(J) = U(J) - A(J4) * U(K)
NEXT K
J0 = J0 + J
NEXT J
RETURN
2280 :
REM ==BACK-SUBSTITUTION==
J0 = 0
FOR J = 1 TO N
J0 = J0 + J
U(J) = U(J) / A(J0)
NEXT J
J = N
J0 = N * (N + 1) / 2
2350 :
IF (J <= 1) THEN 2440
J2 = J - 1
J0 = J0 - J
FOR I = 1 TO J2
I1 = J0 + I
U(I) = U(I) - A(I1) * U(J)
NEXT I
J = J2
GOTO 2350
2440 :
RETURN
2450 :
REM ==VECTOR B INPUT==
PRINT "INPUT VECTOR B"
PRINT
FOR I = 1 TO N
PRINT "B("; I; ") =  ";
INPUT B(I)
NEXT I
PRINT
RETURN
2550 :
JJ = 0
FOR J = 1 TO N
JJ = JJ + J
A(JJ) = A(JJ) - W0 * B(J)
NEXT J
RETURN
2610 :
REM ==FULL MATRIX B INPUT==
PRINT "ENTER MATRIX B UPPER TRIANGLE"
J0 = 0
FOR J = 1 TO N
2640 :
PRINT
PRINT "UPPER PART COLUMN "; J
PRINT
FOR I = 1 TO J
I1 = I + J0
PRINT "B("; I; ","; J; ") =  ";
INPUT B(I1)
NEXT I
PRINT
INPUT "IS IT OK? (Y/N)"; W$
IF W$ = "N" THEN 2640
J0 = J0 + J
NEXT J
RETURN
2780 :
JJ = 0
FOR J = 1 TO N
FOR I = 1 TO J
IJ = JJ + I
A(IJ) = A(IJ) - W0 * B(IJ)
NEXT I
JJ = JJ + J
NEXT J
RETURN
2870 :
FOR I = 1 TO N
Z(I) = B(I) * U(I)
NEXT I
RETURN
2910 :
REM B IS A FULL MATRIX
J0 = 0
N1 = N - 1
FOR J = 1 TO N1
Z(J) = 0
FOR K = 1 TO J
J1 = J0 + K
Z(J) = Z(J) + B(J1) * U(K)
NEXT K
J0 = J0 + J
K0 = J + 1
J1 = J0
FOR K = K0 TO N
J1 = J1 + K - 1
Z(J) = Z(J) + B(J1) * U(K)
NEXT K
NEXT J
Z(N) = 0
FOR K = 1 TO N
J1 = J0 + K
Z(N) = Z(N) + B(J1) * U(K)
NEXT K
RETURN