Procedure HQR(A:matrx2;N:integer;var WR,WI:array of real);
Label 1,2,3,4;
var
    I,J,NN,III,K,ITS,L,M:integer;
    ANORM,T,S,X,Y,W,R,Q,P,U,ZZ,Z,AAA,BBB,V:real;
begin
    ANORM:=Abs(A[1, 1]);
    For I:=2 To N do
        For J:=I - 1 To N do
            ANORM:=ANORM + Abs(A[I, J]);
    NN:=N;
    T:=0;
1:  If NN >= 1 Then
    begin
        ITS:=0;
2:      For L:=NN DownTo 2 do
        begin
            S:=Abs(A[L - 1, L - 1]) + Abs(A[L, L]);
            If S = 0  Then S:=ANORM;
            If Abs(A[L, L - 1]) + S = S Then GoTo 3;
        end;
        L:=1;
3:      X:=A[NN, NN];
        If L = NN Then
        begin
            WR[NN]:=X + T;
            WI[NN]:=0;
            NN:=NN - 1;
        end
        Else
        begin
            Y:=A[NN - 1, NN - 1];
            W:=A[NN, NN - 1] * A[NN - 1, NN];
            If L = NN - 1 Then
            begin
                P:=0.5 * (Y - X);
                Q:=P*P + W;
                Z:=Sqrt(Abs(Q));
                X:=X + T;
                If Q >= 0 Then
                begin
                    if P>=0 then
                      ZZ:=1
                    else
                      ZZ:=-1;
                    Z:=P + Abs(Z) * ZZ;
                    WR[NN]:=Z + X;
                    WR[NN - 1]:=WR[NN];
                    If Z <> 0 Then WR[NN]:=X - W / Z;
                    WI[NN]:=0;
                    WI[NN - 1]:=0;
                end
                Else
                begin
                    WR[NN]:=X + P;
                    WR[NN - 1]:=WR[NN];
                    WI[NN]:=Z;
                    WI[NN - 1]:=-Z;
                end;
                NN:=NN - 2;
            end
            Else
            begin
                If ITS = 30 Then ShowMessage(' too many iterations ');
                If (ITS = 10) Or (ITS = 20) Then
                begin
                    T:=T + X;
                    For I:=1 To NN do
                        A[I, I]:=A[I, I] - X;
                    S:=Abs(A[NN, NN - 1]) + Abs(A[NN - 1, NN - 2]);
                    X:=0.75 * S;
                    Y:=X;
                    W:=-0.4375 * S * S;
                end;
                ITS:=ITS + 1;
                For M:=NN - 2 DownTo L do
                begin
                    Z:=A[M, M];
                    R:=X - Z;
                    S:=Y - Z;
                    P:=(R * S - W) / A[M + 1, M] + A[M, M + 1];
                    Q:=A[M + 1, M + 1] - Z - R - S;
                    R:=A[M + 2, M + 1];
                    S:=Abs(P) + Abs(Q) + Abs(R);
                    P:=P / S;
                    Q:=Q / S;
                    R:=R / S;
                    If M = L Then GoTo 4;
                    U:=Abs(A[M, M - 1]) * (Abs(Q) + Abs(R));
                    BBB:=Abs(A[M + 1, M + 1]);
                    AAA:=Abs(A[M - 1, M - 1]) + Abs(Z) + BBB;
                    V:=Abs(P) * AAA;
                    If U + V = V Then GoTo 4;
                end;
4:              For I:=M + 2 To NN do
                begin
                    A[I, I - 2]:=0;
                    If I <> M + 2 Then A[I, I - 3]:=0;
                end;
                For K:=M To NN - 1 do
                begin
                    If K <> M Then
                    begin
                        P:=A[K, K - 1];
                        Q:=A[K + 1, K - 1];
                        R:=0;
                        If K <> NN - 1 Then R:=A[K + 2, K - 1];
                        X:=Abs(P) + Abs(Q) + Abs(R);
                        If X <> 0 Then
                        begin
                            P:=P / X;
                            Q:=Q / X;
                            R:=R / X;
                        end;
                    end;
                    if P>=0 then
                      ZZ:=1
                    else
                      ZZ:=-1;
                    S:=Sqrt(P*P + Q*Q + R*R) * ZZ;
                    If S <> 0 Then
                    begin
                        If K = M Then
                        begin
                            If L <> M Then
                              A[K, K - 1]:=-A[K, K - 1];
                        end
                        Else
                            A[K, K - 1]:=-S * X;
                        P:=P + S;
                        X:=P / S;
                        Y:=Q / S;
                        Z:=R / S;
                        Q:=Q / P;
                        R:=R / P;
                        For J:=K To NN do
                        begin
                            P:=A[K, J] + Q * A[K + 1, J];
                            If K <> NN - 1 Then
                            begin
                                P:=P + R * A[K + 2, J];
                                A[K + 2, J]:=A[K + 2, J] - P * Z;
                            end;
                            A[K + 1, J]:=A[K + 1, J] - P * Y;
                            A[K, J]:=A[K, J] - P * X;
                        end;
                        If NN > K + 3 Then
                           III:=K + 3
                        Else
                           III:=NN;
                        For I:=L To III do
                        begin
                            P:=X * A[I, K] + Y * A[I, K + 1];
                            If K <> NN - 1 Then
                            begin
                                P:=P + Z * A[I, K + 2];
                                A[I, K + 2]:=A[I, K + 2] - P * R;
                            end;
                            A[I, K + 1]:=A[I, K + 1] - P * Q;
                            A[I, K]:=A[I, K] - P;
                        end;
                    end;
                end;
                GoTo 2;
            end;
        end;
        GoTo 1;
    end;
end;