CLS
DIM A(120)
DIM B(120)
DIM A1(120)
DIM U(15)
DIM Z(15)
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 = 1/10^8
PRINT "ENTER TYPE OF PROBLEM"
PRINT
INPUT "TYPE = ";I9
PRINT
PRINT
PRINT "DIMENSION MATRIX A"
PRINT
INPUT "N = ?";N
J0 = 0
PRINT
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 AAA
A1(I1)=AAA
NEXT I
PRINT
J0 = J0 + J
NEXT J
PRINT
PRINT
PRINT
REM ==INPUT MATRIX OR VECTOR B==
IF I9=1 THEN GOSUB [TWFOFIZE]
IF I9=2 THEN GOSUB [TWSIONZE]
PRINT
PRINT
PRINT ">ENTER TOLERANCE"
PRINT
INPUT "ORDER OF TOLERANCE X (AS IN 1e-X)";T0
T0=1/10^T0
I8 = 0
[FIFOZE]
PRINT
PRINT
PRINT ">MAXIMUM # OF ITERATIONS (PRESS ENTER FOR 50)"
PRINT
INPUT "NMAX = ";N8$
IF N8$="" THEN
 N8=50
ELSE
 N8=VAL(N8$)
END IF
IF I8 = 2 THEN [NISEZE]
[SIZEZE]
PRINT
PRINT ">INPUT DESIRED SHIFT"
PRINT
INPUT "SHIFT = ";W0
PRINT
PRINT "RUNNING"
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 [TWFIFIZE]
IF I9 = 2 THEN GOSUB [TWSEEIZE]
GOSUB [ONSEFOZE]
REM ==CHECK POLARITY==
IF I8 = 1 THEN [SIZEZE]
M0 = 0
II = 0
REM ==CHECK STURM SEQUENCE==
FOR I=1 TO N
II = II + I
IF A(II) > 0 THEN [EITWZE]
M0 = M0 + 1
[EITWZE]
NEXT I
PRINT
PRINT "NUMBER OF EIGENVALUES"
PRINT
PRINT "SMALLER THAN SHIFT = ";M0
PRINT
PRINT "DETERMINANT = ";D8
[NIZEZE]
PRINT
PRINT "WANT SHIFT CHANGED (Y/N)";
INPUT W$
PRINT
PRINT "RUNNING"
IF W$ ="Y" THEN [SIZEZE]
IF I8 = 2 THEN [FIFOZE]
REM ==PERFORM INVERSE ITERATIONS==
[NISEZE]
GOSUB [ONTWFIZE]
REM ==CHECK FOR CONVERGENCY==
IF I8 = 2 THEN [NIZEZE]
REM ==UNDO SHIFT==
E0 = E0 + W0
GOSUB [ONZESIZE]
PRINT
PRINT "WANT ANOTHER EIGENSOLUTION (Y/N) ";
INPUT W$
IF W$ = "Y" GOTO [FIFOZE]
END
REM OUTPUT M0DULE
[ONZESIZE]
PRINT
PRINT ">SOLUTION"
PRINT" ********"
PRINT
PRINT "EIGENVALUE"
PRINT "EIG = ";E0
PRINT
PRINT "EIGENVECTOR"
FOR I = 1 TO N
PRINT "U(";I; ")= ";U(I)
NEXT I
PRINT
PRINT "NUM. ITERATIONS"
PRINT "N = ";N0
PRINT
RETURN
REM ==INVERSE ITERATION SUBROUTINE==
[ONTWFIZE]
I8 = 0
REM ==INITIAL EIGENVECTOR==
II =0
FOR I = 1 TO N
II = II + I
B1 = B(I)
IF I9 = 1 THEN [ONTHTWZE]
B1 = B(II)
[ONTHTWZE]
U(I) = A(II) / B1
NEXT I
E1 = 1*10^32
N0 = 0
R1 = 0
GOTO [ONFIZEZE]
REM ==SET ITERATION COUNTER==
[ONTHEIZE]
N0 = N0+ 1
FOR I = 1 TO N
U(I) = Z(I)
NEXT I
REM ==FORWARD SOLUTION==
GOSUB [TWONEIZE]
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 [TWTWEIZE]
REM ==COMPUTE B-NORM OF U==
[ONFIZEZE]
IF I9 = 1 THEN GOSUB [TWEISEZE]
IF I9 = 2  THEN GOSUB [TWNIONZE]
V2 = 0
FOR I = 1 TO N
V2 = V2 + Z(I) * U(I)
NEXT I
V3 = V2^.5
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 [ONSEZEZE]
E1 = E0
IF N0 < N8 THEN [ONTHEIZE]
REM ==ERROR MESSAGES==
IS = 2
PRINT
PRINT "**WARNING** CONVERGENCE NOT"
PRINT "ATTAINED. CHANGE NMAX OR SHIFT"
RETURN
REM ==FINAL EIGENVECTOR==
[ONSEZEZE]
FOR I = 1 TO N
U(I) = U(I) / V3
NEXT I
RETURN
REM ==FACTORIZATION OF A ==
[ONSEFOZE]
I8 = 0
J0 = 1
D8 = A(J0)
REM LOOP OVER COLUMNS
FOR J = 2 TO N
J2 = J-1
IF J = 2 THEN [ONNIONZE]
I0 = 1
REM: LOOP OVER ROWS
FOR I = 2 TO J2
I2 = I -1
I1 = J0 + I
REM REDUCE COLUMN I1
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==
[ONNIONZE]
IF ABS(A(J0))< E5 THEN [TWONTWZE]
REM ==REDUCE DIAGONAL ELEMENT AND COLUMN J==
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 FORZERO PIVOT==
IF ABS(A(J0)) > E5 THEN [TWONSIZE]
REM ERROR MESSAGES
PRINT
PRINT "**WARNING** ZERO DIAGONAL"
PRINT "ELEMENT.CHANGESHIFT"
PRINT
RETURN
[TWONTWZE]
PRINT "**ERROR** UNSUCCESSFUL FACTORIZATION"
PRINT "CHOOSE NEW SHIFT"
I8 = 1
PRINT
[TWONSIZE]
RETURN
REM ==FORWARD SOLUTION==
[TWONEIZE]
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
REM ==BACK-SUBSTITUTION==
[TWTWEIZE]
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
[TWTHFIZE]
IF (J <=1) THEN [TWFOFOZE]
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 [TWTHFIZE]
[TWFOFOZE]
RETURN
REM ==VECTOR B INPUT==
[TWFOFIZE]
PRINT "INPUT VECTOR B"
PRINT
FOR I = 1 TO N
PRINT "B(";I;") =";
INPUT BBB
B(I)=BBB
NEXT I
PRINT
RETURN
[TWFIFIZE]
JJ = 0
FOR J = 1 TO N
JJ = JJ + J
A(JJ) = A(JJ)-W0*B(J)
NEXT J
RETURN
REM ==FULL MATRIX B INPUT==
[TWSIONZE]
PRINT "ENTER MATRIX B UPPER TRIANGLE"
J0 = 0
FOR J = 1 TO N
[TWSIFOZE]
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 [TWSIFOZE]
J0 = J0 + J
NEXT J
RETURN
[TWSEEIZE]
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
[TWEISEZE]
FOR I = 1 TO N
Z(I) = B(I) * U(I)
NEXT I
RETURN
REM ==B IS A FULL MATRIX==
[TWNIONZE]
J0 = 0
N1 = N-1
FOR J =1 TO N1
Z(J) = 0
FOR K = 1 TO J
J1 = J0 + K
Z(J) = I(J) + B(J1) * U(K)
NEXT K
J0 = J0 + J
K0 = J + 1
J1 = J0
FOR K = K0 TO N
J1 = J1 + K -1
I(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