'**************************************************** '* Calculate the eigenvalues and the eigenvectors * '* of a real symmetric tridiagonal matrix * '* ------------------------------------------------ * '* Ref.: "NUMERICAL RECIPES, Cambridge University * '* Press, 1986" [BIBLI 08]. * '* * '* SAMPLE RUN: * '* (find eigenvalues and eigenvectors of matrix: * '* 1 2 0 0 0 * '* 2 4 7 0 0 * '* A = 0 7 10 8 0 * '* 0 0 8 -0.75 -9 * '* 0 0 0 -9 10 ) * '* * '* Eigenvalues 1: -0.78071442 * '* Eigenvector: * '* 1.000000 -0.890357 0.322363 0.344649 0.287721 * '* * '* Eigenvalues 2: 2.53046878 * '* Eigenvector: * '* 1.000000 0.765234 -0.446362 -0.252815 -0.304616 * '* * '* Eigenvalues 3: -9.15659229 * '* Eigenvector: * '*-0.056408 0.286458 -0.522285 1.000000 0.469812 * '* * '* Eigenvalues 4: 12.53989066 * '* Eigenvector: * '* 0.097161 0.560616 0.656182 -0.282210 1.000000 * '* * '* Eigenvalues 5: 19.11694726 * '* Eigenvector: * '* 0.051876 0.469920 1.000000 0.728439 -0.719095 * '* * '* Basic version by J-P Moreau, Paris * '* (www.jpmoreau.fr) * '**************************************************** 'PROGRAM TEST_TQLI defdbl a-h,o-z defint i-n NMAX=20 dim D(NMAX), E(NMAX), Z(NMAX,NMAX) cls print input " Input size of matrix: ", n gosub 1000 'read data 'initialize Z to identity matrix for i=1 to n for j=1 to n if i<>j then Z(i,j) = 0# else Z(i,j) = 1# end if next j next i gosub 2000 'call TQLI(D,E,n,Z) 'print results for j=1 to n print print " Eigenvalue ";j;": "; print using "####.######"; D(j) print " Eigenvector:" 'normalize vector before printing xmax=Z(1,j) for i=2 to n if ABS(Z(i,j))>ABS(xmax) then xmax=Z(i,j) next i for i=1 to n print using "####.######"; Z(i,j)/xmax; next i print next j print print END 'of main program 'read vectors D and E from screen 'Subroutine Read_data(n,D,E) 1000 print print " Input "; n-1; " elements of subdiagonal:" E(1)=1# 'arbitrary value for i=1 to n-1 print " Element ";i; ": "; : input "", E(i+1) next i print " Input "; n; " elements of main diagonal:" for i=1 to n print " Element ";i; ": "; : input "", D(i) next i print return 1500 'y=Sign(r,g) if g>=0 then y=ABS(r) else y=-ABS(r) end if return '**************************************************** '* This subroutine implements the QL algorithm with * '* implicit shifts to determine the eigenvalues and * '* eigenvectors of a real symmetric tridiagonal * '* matrix. * '* ------------------------------------------------ * '* INPUTS: * '* D Elements of main diagonal * '* E Elements of subdiagonal (from E(2) * '* to E(n)) * '* n size of matrix * '* Z identity matrix * '* OUTPUTS * '* D vector storing the n eigenvalues * '* Z matrix storing the n eigenvectors * '* in columns. * '**************************************************** 'Subroutine TQLI(D,E,n,Z) 2000 if n>1 then for i=2 to n E(i-1)=E(i) next i E(n)=0# for l=1 to n iter=0 1 for m=l to n-1 dd=ABS(D(m))+ABS(D(m+1)) if ABS(E(m))+dd=dd then goto 2 next m m=n 2 if m<>l then if iter=30 then print " 30 iterations." return end if iter=iter+1 g=(D(l+1)-D(l))/(2#*E(l)) r=SQR(g*g+1#) gosub 1500 'y=Sign(r,g) g=D(m)-D(l)+E(l)/(g+y) s=1# : c=1# : p=0# for i=m-1 to l step -1 f=s*E(i) : b=c*E(i) if ABS(f)>=ABS(g) then c=g/f r=SQR(c*c+1#) E(i+1)=f*r s=1#/r c=c*s else s=f/g r=SQR(s*s+1#) E(i+1)=g*r c=1#/r s=s*c end if g=D(i+1)-p r=(D(i)-g)*s+2#*c*b p=s*r D(i+1)=g+p g=c*r-b for k=1 to n f=Z(k,i+1) Z(k,i+1)=s*Z(k,i)+c*f Z(k,i)=c*Z(k,i)-s*f next k next i D(l)=D(l)-p E(l)=g : E(m)=0# goto 1 end if next l end if return 'end of file elprotd.bas