unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls,unit1, Forms, Dialogs;
Procedure TRED2(var A:matrx2; N:integer;
                       var D:array of real;var E:array of real);
Procedure TQLI(var D:array of real; E:array of real; N:integer;var Z:matrx2);

implementation
Procedure TRED2(var A:matrx2; N:integer;
                       var D:array of real;var E:array of real);
var
     I,J,K,L:integer; H,F,G,SCALE1,ZZ,HH:real;
begin
     If N > 1 Then
     begin
        For I:=N DownTo 2 do
        begin
            L:=I - 1;
            H:=0;
            SCALE1:=0;
            If L > 1 Then 
            Begin 
                For K:=1 To L do
                    SCALE1:=SCALE1 + Abs(A[I, K]);
                If SCALE1 = 0  Then
                    E[I]:=A[I, L]
                else
                begin
                    For K:=1 To L do
                    begin
                        A[I, K]:=A[I, K] / SCALE1;
                        H:=H + Sqr(A[I, K]);
                    end;
                    F:=A[I, L];
                    If F >= 0 then
                        ZZ:=1
                    Else
                        ZZ:=-1;
                    G:=-Sqrt(H) * ZZ;
                    E[I]:=SCALE1 * G;
                    H:=H - F * G;
                    A[I, L]:=F - G;
                    F:=0; 
                    For J:=1 To L do
                    begin
                        A[J, I]:=A[I, J] / H;
                        G:=0;
                        For K:=1 To J do
                            G:=G + A[J, K] * A[I, K];
                        If L > J Then
                        begin
                            For K:=J + 1 To L do
                                G:=G + A[K, J] * A[I, K];
                        End;
                        E[J]:=G / H;
                        F:=F + E[J] * A[I, J];
                    end; 
                    HH:=F / (H + H);
                    For J:=1 To L do
                    begin
                        F:=A[I, J];
                        G:=E[J] - HH * F;
                        E[J]:=G;
                        For K:=1 To J do
                            A[J, K]:=A[J, K] - F * E[K] - G * A[I, K];
                    end; 
                end;
            end
            else
                E[I]:=A[I, L];
            D[I]:=H;
        end;
    end; 
    //Omit following line if finding only eigenvalues.
    D[1]:=0;
    E[1]:=0;
    For I:=1 To N do
    begin
    //Delete lines from here ...
        L:=I - 1;
        If D[I] <> 0  Then
        begin
            For J:=1 To L do
            begin
                G:=0;
                For K:=1 To L do
                    G:=G + A[I, K] * A[K, J];
                For K:=1 To L do
                    A[K, J]:=A[K, J] - G * A[K, I];
            end; 
        end; 
     //... to here when finding only eibenvalues.
        D[I]:=A[I, I];
    //Also delete lines from here ...
        A[I, I]:=1; 
        If L >= 1 Then
        begin
            For J:=1 To L do
            begin
                A[I, J]:=0;
                A[J, I]:=0;
            end;
        end;
    //... to here when finding only eigenvalues.
    end;
end;

Procedure TQLI(var D:array of real; E:array of real; N:integer;var Z:matrx2);
Label 1,2;
Var
    I,L,M,K,ITER:integer;  G,R,C,F,S,P,DD,ZZ,ZZZ,B:real;
begin
    If N > 1 Then
    begin
        For I:=2 To N do
            E[I - 1]:=E[I];
        E[N]:=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 showMessage(' too many iterations ');
                ITER:=ITER + 1;
                G:=(D[L + 1] - D[L]) / (2  * E[L]);
                R:=Sqrt(Sqr(G) + 1 );
                If G >= 0 then
                    ZZ:=1
                Else
                    ZZ:=-1;
                If G >= 0 then
                    ZZZ:=1
                Else
                    ZZZ:=-1;
                G:=D[M] - D[L] + E[L] / (G + ZZZ * ZZ);
                S:=1;
                C:=1;
                P:=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(Sqr(C) + 1 );
                        E[I + 1]:=F * R;
                        S:=1 / R;
                        C:=C * S;
                    end
                    Else
                    begin
                        S:=F / G;
                        R:=Sqrt(Sqr(S) + 1 );
                        E[I + 1]:=G * R;
                        C:=1 / R;
                        S:=S * C;
                    end;
                    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;
                    //Omit lines from here ...
                    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; 
                   //to here when finding only eigenvalues.
                end; 
                D[L]:=D[L] - P;
                E[L]:=G;
                E[M]:=0; 
                GoTo 1;
            end;
       end; 
    end;
end;
end.