1 REM  Eigenvalues and eigenvectors for a symmetric matrix. Eispack algorithm, taken after IMSL's EIGRS (1979)
2 N=5
3 N1=N*(N+1)/2
4 REPS=1
5 IF REPS+1=1 THEN GOTO 8
6 REPS=REPS/2
7 GOTO 5
8 REPS=REPS*2
30 DIM A(121),A1(121),D(11),E(11),Z(11,11)
31 REM  The Matrix is for Hueckel's method for the Butadiene molecule Hamiltonian
32 A(1)=1.4 :: A(2)=.6 :: A(5)=.8 :: A(9)=.6 :: A(11)=.6 :: A(14)=.8
100 FOR I=1 TO N1
110 A1(I)=A(I)
120 NEXT I
130 REM
140 REM
150 N9=N1-1
160 N2=N1-N
170 FOR I1=1 TO N
180 I=N+1-I1
190 L=I-1
200 H=0
210 S=0
220 IF L<1 THEN 290
230 N8=N9
240 FOR K=1 TO L
250 S=S+ABS(A(N8))
260 N8=N8+1
270 NEXT K
280 IF S<>0 THEN 310
290 E(I)=0
300 GOTO 820
310 N8=N9
320 S1=1/S
330 FOR K=1 TO L
340 A(N8)=A(N8)*S1
350 H=H+A(N8)*A(N8)
360 N8=N8-1
370 NEXT K
380 F=A(N9)
390 G=-SQR(H)
400 IF F<0 THEN G=-G
410 E(I)=S*G
420 H=H-F*G
430 A(N9)=F-G
440 IF L=1 THEN 790
450 F=0
460 J2=1
470 FOR J=1 TO L
480 G=0
490 I2=N2+1
500 J1=J2
510 FOR K=1 TO J
520 G=G+A(J1)*A(I2)
530 J1=J1+1
540 I2=I2+1
550 NEXT K
560 J3=J+1
570 IF L<J3 THEN 640
580 J1=J1+J-1
590 FOR K=J3 TO L
600 G=G+A(J1)*A(I2)
610 J1=J1+K
620 I2=I2+1
630 NEXT K
640 E(J)=G/H
650 F=F+E(J)*A(N2+J)
660 J2=J2+J
670 NEXT J
680 H1=F/(H+H)
690 J1=1
700 FOR J=1 TO L
710 F=A(N2+J)
720 G=E(J)-H1*F
730 E(J)=G
740 FOR K=1 TO J
750 A(J1)=A(J1)-F*E(K)-G*A(N2+K)
760 J1=J1+1
770 NEXT K
780 NEXT J
790 FOR K=1 TO L
800 A(N2+K)=S*A(N2+K)
810 NEXT K
820 D(I)=A(N2+I)
830 A(N2+I)=H*S*S
840 N2=N2-I+1
850 N9=N9-I
860 NEXT I1
870 REM
880 FOR I=1 TO N
890 FOR J=1 TO N
900 Z(I,J)=0
910 NEXT J
920 Z(I,I)=1
930 NEXT I
960 R4=REPS
970 IF N=1 THEN 1850
980 FOR I=2 TO N
990 E(I-1)=E(I)
1000 NEXT I
1010 E(N)=0
1020 B=0
1030 F=0
1040 FOR L=1 TO N
1050 J=0
1060 H=R4*(ABS(D(L))+ABS(E(L)))
1070 IF B<H THEN B=H
1080 FOR M=L TO N
1090 K=M
1100 IF ABS(E(K))<=B THEN 1120
1110 NEXT M
1120 M=K
1130 IF M=L THEN 1600
1140 IF J=30 THEN 1820
1150 J=J+1
1160 L1=L+1
1170 G=D(L)
1180 P=(D(L1)-G)/(E(L)+E(L))
1190 R=SQR(P*P+1)
1200 IF P>=0 THEN D(L)=E(L)/(P+R)
1210 IF P<0 THEN D(L)=E(L)/(P-R)
1220 H=G-D(L)
1230 FOR I=L1 TO N
1240 D(I)=D(I)-H
1250 NEXT I
1260 F=F+H
1270 P=D(M)
1280 C=1
1290 S=0
1300 M1=M-1
1310 M2=M1+L
1320 IF L>M1 THEN 1570
1330 FOR I1=L TO M1
1340 I=M2-I1
1350 G=C*E(I)
1360 H=C*P
1370 IF ABS(P)<ABS(E(I))THEN 1440
1380 C=E(I)/P
1390 R=SQR(C*C+1)
1400 E(I+1)=S*P*R
1410 S=C/R
1420 C=1/R
1430 GOTO 1490
1440 C=P/E(I)
1450 R=SQR(C*C+1)
1460 E(I+1)=S*E(I)*R
1470 S=1/R
1480 C=C*S
1490 P=C*D(I)-S*G
1500 D(I+1)=H+S*(C*G+S*D(I))
1510 FOR K=1 TO N
1520 H=Z(K,I+1)
1530 Z(K,I+1)=S*Z(K,I)+C*H
1540 Z(K,I)=C*Z(K,I)-S*H
1550 NEXT K
1560 NEXT I1
1570 E(L)=S*P
1580 D(L)=C*P
1590 IF ABS(E(L))>B THEN 1140
1600 D(L)=D(L)+F
1610 NEXT L
1620 FOR I=1 TO N
1630 K=I
1640 P=D(I)
1650 I3=I+1
1660 IF I3>N THEN 1720
1670 FOR J=I3 TO N
1680 IF D(J)>=P THEN 1710
1690 K=J
1700 P=D(J)
1710 NEXT J
1720 IF K=I THEN 1800
1730 D(K)=D(I)
1740 D(I)=P
1750 FOR J=1 TO N
1760 P=Z(J,I)
1770 Z(J,I)=Z(J,K)
1780 Z(J,K)=P
1790 NEXT J
1800 NEXT I
1810 GOTO 1850
1820 PRINT "error";L
1830 REM
1840 REM
1850 IF N=1 THEN 2020
1860 FOR I=2 TO N
1870 L=I-1
1880 I1=I*L/2
1890 H=A(I1+I)
1900 IF H=0 THEN 2010
1910 FOR J=1 TO N
1920 S=0
1930 FOR K=1 TO L
1940 S=S+A(I1+K)*Z(K,J)
1950 NEXT K
1960 S=S/H
1970 FOR K=1 TO L
1980 Z(K,J)=Z(K,J)-S*A(I1+K)
1990 NEXT K
2000 NEXT J
2010 NEXT I
2020 REM
2050 REM
2080 FOR J=1 TO N
2100 CALL CLEAR :: PRINT "Eigenvalue and eigenvector";J :: PRINT "Eig.val = ",D(J)
2110 !PRINT "Eigenvector"
2120 FOR I=1 TO N :: PRINT Z(I,J)
2130 NEXT I
2140 R=0
2150 FOR I=1 TO N
2160 S=0
2170 FOR K=1 TO N
2180 IF I<K THEN F=A1(K*(K-1)/2+I)
2190 IF I>=K THEN F=A1(I*(I-1)/2+K)
2200 IF I=K THEN F=F-D(J)
2210 S=S+F*Z(K,J)
2220 NEXT K
2230 IF R<ABS(S)THEN R=ABS(S)
2240 NEXT I
2250 PRINT "rezidue";R :: PRINT :: PRINT "Press any key to continue or <Ctrl><BREAK> to finish"
2255 INPUT A$
2260 NEXT J
2265 GOTO 2080
2270 END