{*************************************************** * 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 * * * * Pascal version by J-P Moreau, Paris * * (www.jpmoreau.fr) * ***************************************************} PROGRAM TEST_TQLI; Uses WinCrt; Const NMAX = 30; Type pM = ^MAT; MAT = Array[1..NMAX,1..NMAX] of real; pV = ^VEC; VEC = Array[1..NMAX] of real; Var D, E : pV; Z : pM; i,j,n: integer; max : real; Function Sign(a,b:real):real; begin if b>=0 then Sign:=ABS(a) else Sign:=-ABS(a) end; {read vectors D and E from screen} Procedure Read_data(VAR n:integer; VAR D,E:pV); var i:integer; Begin writeln; write(' Input size of matrix: '); readln(n); writeln(' Input ',n-1,' elements of subdiagonal:'); E^[1]:=1.0; {arbitrary value} for i:=1 to n-1 do begin write(' Element ',i,': '); readln(E^[i+1]) end; writeln(' Input ',n,' elements of main diagonal:'); for i:=1 to n do begin write(' Element ',i,': '); readln(D^[i]) end End; {*************************************************** * 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. * ***************************************************} Procedure TQLI(Var D,E:pV;n:integer; VAR Z:pM); Label 1,2,fin; Var b,c,dd,f,g,p,r,s: real; i,iter,k,l,m:integer; Begin if n>1 then begin for i:=2 to n do E^[i-1]:=E^[i]; E^[n]:=0.0; for l:=1 to n do begin iter:=0; 1: for m:=l to n-1 do begin dd:=ABS(D^[m])+ABS(D^[m+1]); if ABS(E^[m])+dd=dd then goto 2 end; m:=n; 2: if m<>l then begin if iter=30 then goto fin; iter:=iter+1; g:=(D^[l+1]-D^[l])/(2.0*E^[l]); r:=SQRT(g*g+1.0); g:=D^[m]-D^[l]+E^[l]/(g+SIGN(r,g)); s:=1.0; c:=1.0; p:=0.0; for i:=m-1 Downto l do begin f:=s*E^[i]; b:=c*E^[i]; if ABS(f)>=ABS(g) then begin c:=g/f; r:=SQRT(c*c+1.0); E^[i+1]:=f*r; s:=1.0/r; c:=c*s end else begin s:=f/g; r:=SQRT(s*s+1.0); E^[i+1]:=g*r; c:=1.0/r; s:=s*c end; g:=D^[i+1]-p; r:=(D^[i]-g)*s+2.0*c*b; p:=s*r; D^[i+1]:=g+p; g:=c*r-b; for k:=1 to n do begin f:=Z^[k,i+1]; Z^[k,i+1]:=s*Z^[k,i]+c*f; Z^[k,i]:=c*Z^[k,i]-s*f end end; D^[l]:=D^[l]-p; E^[l]:=g; E^[m]:=0.0; goto 1 end end; {of l loop} fin:if iter>0 then writeln(' ',iter,' iterations.'); end end; {main program} BEGIN {allocate vectors and matrix} New(D); New(E); New(Z); read_data(n,D,E); {initialize Z to identity matrix} for i:=1 to n do for j:=1 to n do if i<>j then Z^[i,j] := 0.0 else Z^[i,j] := 1.0; {calculate eigenvalues and eigenvectors} TQLI(D,E,n,Z); Dispose(E); {print results} for j:=1 to n do begin writeln; writeln(' Eigenvalue ',j,': ', D^[j]:12:8); writeln(' Eigenvector:'); max:=ABS(Z^[1,j]); for i:=2 to n do if ABS(Z^[i,j])>ABS(max) then max:=Z^[i,j]; for i:=1 to n do write(' ',(Z^[i,j]/max):10:6); writeln end; writeln; readkey; Dispose(D); Dispose(Z); DoneWinCrt END. {end of file elprotd.pas}