begin
                    Y := Y / X;
                    A[I,Mm1] := Y;
                    for J := M to Ubound do
                      A[I,J] := A[I,J] - Y * A[M,J];
                    for J := Lbound to I_igh do
                      A[J,M] := A[J,M] + Y * A[J,I];
                  end;
              end;
          end;
      end;
  end;

  procedure Eltran(A : TMatrix; Lbound, Ubound, I_low, I_igh : Integer;
                   I_int : TIntVector; Z : TMatrix);
{ ----------------------------------------------------------------------
  This procedure is a translation of the EISPACK subroutine Eltran.

  This procedure accumulates the stabilized elementary similarity
  transformations used in the reduction of a real general matrix
  to upper Hessenberg form by Elmhes.

  On input:

    A contains the multipliers which were used in the reduction
    by Elmhes in its lower triangle below the subdiagonal.

    Lbound, Ubound are the lowest and highest indices
    of the elements of A

    I_low and I_igh are integers determined by the balancing procedure
    Balance. If Balance has not been used, set I_low=Lbound, I_igh=Ubound.

    I_int contains information on the rows and columns interchanged in
    the reduction by Elmhes. Only elements I_low through I_igh are used.

  On output:

    Z contains the transformation matrix produced in the reduction by
    Elmhes.
  ---------------------------------------------------------------------- }
  var
    I, J, Mp, Mp1 : Integer;
  begin
    { Initialize Z to identity matrix }
    for I := Lbound to Ubound do
      for J := Lbound to Ubound do
        if I = J then Z[I,J] := 1.0 else Z[I,J] := 0.0;

    if I_igh < I_low then Exit;

    for Mp := I_igh - 1 downto I_low + 1 do
      begin
        Mp1 := Mp + 1;
        for I := Mp1 to I_igh do
          Z[I,Mp] := A[I,Mp - 1];
        I := I_int[Mp];
        if I <> Mp then
          begin
            for J := Mp to I_igh do
              begin
                Z[Mp,J] := Z[I,J];
                Z[I,J] := 0.0;
              end;
            Z[I,Mp] := 1.0;
          end;
      end;
  end;

  function Hqr(H : TMatrix; Lbound, Ubound, I_low, I_igh : Integer;
               Wr, Wi : TVector) : Integer;
{ ----------------------------------------------------------------------
  This function is a translation of the EISPACK subroutine hqr.

  This function finds the eigenvalues of a real upper Hessenberg
  matrix by the QR method.

  On input:

    H contains the upper Hessenberg matrix.

    Lbound, Ubound are the lowest and highest indices
    of the elements of H

    I_low and I_igh are integers determined by the balancing subroutine
    Balance. If Balance has not been used, set I_low=Lbound, I_igh=Ubound

  On output:

    H has been destroyed.

    Wr and Wi contain the real and imaginary parts, respectively, of
    the eigenvalues. The eigenvalues are unordered except that complex
    conjugate pairs of values appear consecutively with the eigenvalue
    having the positive imaginary part first.

    The function returns an error code:
       zero       for normal return,
       -j         if the limit of 30*N iterations is exhausted
                  while the j-th eigenvalue is being sought.
                  (N being the size of the matrix). The eigenvalues
                  should be correct for indices j+1,...,Ubound.
  ----------------------------------------------------------------------
  Note: This is a crude translation. Many of the original goto's have
  been kept !
  ---------------------------------------------------------------------- }
  var
    I, J, K, L, M, N, En, Na, Itn, Its, Mp2, Enm2 : Integer;
    P, Q, R, S, T, W, X, Y, Zz, Norm, Tst1, Tst2 : Float;
    NotLas : Boolean;

  label
    60, 70, 100, 130, 150, 170, 225, 260, 270, 280, 320, 330;

  begin
    { Store roots isolated by Balance and compute matrix norm }
    K := Lbound;
    Norm := 0.0;
    for I := Lbound to Ubound do
      begin
        for J := K to Ubound do
          Norm := Norm + Abs(H[I,J]);
        K := I;
        if (I < I_low) or (I > I_igh) then
          begin
            Wr[I] := H[I,I];
            Wi[I] := 0.0;
          end;
      end;

    N := Ubound - Lbound + 1;
    Itn := 30 * N;
    En := I_igh;
    T := 0.0;

60: { Search for next eigenvalues }
    if En < I_low then
      begin
        Hqr := 0;
        Exit;
      end;

    Its := 0;
    Na := En - 1;
    Enm2 := Na - 1;

70: { Look for single small sub-diagonal element }
    for L := En downto I_low do
      begin
        if L = I_low then goto 100;
        S := Abs(H[L - 1,L - 1]) + Abs(H[L,L]);
        if S = 0.0 then S := Norm;
        Tst1 := S;
        Tst2 := Tst1 + Abs(H[L,L - 1]);
        if Tst2 = Tst1 then goto 100;
      end;

100: { Form shift }
    X := H[En,En];
    if L = En then goto 270;
    Y := H[Na,Na];
    W := H[En,Na] * H[Na,En];
    if L = Na then goto 280;

    if Itn = 0 then
      begin
        { Set error -- all eigenvalues have not
          converged after 30*N iterations }
        Hqr := - En;
        Exit;
      end;

    if (Its <> 10) and (Its <> 20) then goto 130;

    { Form exceptional shift }
    T := T + X;

    for I := I_low to En do
      H[I,I] := H[I,I] - X;

    S := Abs(H[En,Na]) + Abs(H[Na,Enm2]);
    X := 0.75 * S;
    Y := X;
    W := - 0.4375 * S * S;

130:
    Its := Its + 1;
    Itn := Itn - 1;

    { Look for two consecutive small sub-diagonal elements }
    for M := Enm2 downto L do
      begin
        Zz := H[M,M];
        R := X - Zz;
        S := Y - Zz;
        P := (R * S - W) / H[M + 1,M] + H[M,M + 1];
        Q := H[M + 1,M + 1] - Zz - R - S;
        R := H[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 150;
        Tst1 := Abs(P) * (Abs(H[M - 1,M - 1]) + Abs(Zz) + Abs(H[M + 1,M + 1]));
        Tst2 := Tst1 + Abs(H[M,M - 1]) * (Abs(Q) + Abs(R));
        if Tst2 = Tst1 then goto 150;
      end;

150:
    Mp2 := M + 2;

    for I := Mp2 to En do
      begin
        H[I,I - 2] := 0.0;
        if I <> Mp2 then H[I,I - 3] := 0.0;
      end;

    { Double QR step involving rows L to En and columns M to En }
    for K := M to Na do
      begin
        NotLas := (K <> Na);
        if (K = M) then goto 170;
        P := H[K,K - 1];
        Q := H[K + 1,K - 1];
        R := 0.0;
        if NotLas then R := H[K + 2,K - 1];
        X := Abs(P) + Abs(Q) + Abs(R);
        if X = 0.0 then goto 260;
        P := P / X;
        Q := Q / X;
        R := R / X;
170:    S := DSgn(Sqrt(P * P + Q * Q + R * R), P);
        if K <> M then
          H[K,K - 1] := - S * X
        else if L <> M then
          H[K,K - 1] := - H[K,K - 1];
        P := P + S;
        X := P / S;
        Y := Q / S;
        Zz := R / S;
        Q := Q / P;
        R := R / P;
        if NotLas then goto 225;

        { Row modification }
        for J := K to En do
          begin
            P := H[K,J] + Q * H[K + 1,J];
            H[K,J] := H[K,J] - P * X;
            H[K + 1,J] := H[K + 1,J] - P * Y;
          end;

        J := Min(En, K + 3);

        { Column modification }
        for I := L to J do
          begin
            P := X * H[I,K] + Y * H[I,K + 1];
            H[I,K] := H[I,K] - P;
            H[I,K + 1] := H[I,K + 1] - P * Q;
          end;
        goto 260;

225:
        { Row modification }
        for J := K to En do
          begin
            P := H[K,J] + Q * H[K + 1,J] + R * H[K + 2,J];
            H[K,J] := H[K,J] - P * X;
            H[K + 1,J] := H[K + 1,J] - P * Y;
            H[K + 2,J] := H[K + 2,J] - P * Zz;
          end;

        J := Min(En, K + 3);

        { Column modification }
        for I := L to J do
          begin
            P := X * H[I,K] + Y * H[I,K + 1] + Zz * H[I,K + 2];
            H[I,K] := H[I,K] - P;
            H[I,K + 1] := H[I,K + 1] - P * Q;
            H[I,K + 2] := H[I,K + 2] - P * R;
          end;

260:  end;

    goto 70;

270: { One root found }
    Wr[En] := X + T;
    Wi[En] := 0.0;
    En := Na;
    goto 60;

280: { Two roots found }
    P := 0.5 * (Y - X);
    Q := P * P + W;
    Zz := Sqrt(Abs(Q));
    X := X + T;
    if Q < 0.0 then goto 320;

    { Real pair }
    Zz := P + DSgn(Zz, P);
    Wr[Na] := X + Zz;
    Wr[En] := Wr[Na];
    if Zz <> 0.0 then Wr[En] := X - W / Zz;
    Wi[Na] := 0.0;
    Wi[En] := 0.0;
    goto 330;

320: { Complex pair }
    Wr[Na] := X + P;
    Wr[En] := X + P;
    Wi[Na] := Zz;
    Wi[En] := - Zz;

330:
    En := Enm2;
    goto 60;
  end;

  function Hqr2(H : TMatrix; Lbound, Ubound, I_low, I_igh : Integer;
                Wr, Wi : TVector; Z : TMatrix) : Integer;
{ ----------------------------------------------------------------------
  This function is a translation of the EISPACK subroutine hqr2

  This procedure finds the eigenvalues and eigenvectors of a real
  upper Hessenberg matrix by the QR method.

  On input:

    H contains the upper Hessenberg matrix.

    Lbound, Ubound are the lowest and highest indices
    of the elements of H

    I_low and I_igh are integers determined by the balancing subroutine
    Balance. If Balance has not been used, set I_low=Lbound, I_igh=Ubound

    Z contains the transformation matrix produced by Eltran after the
    reduction by Elmhes, or by Ortran after the reduction by Orthes, if
    performed. If the eigenvectors of the Hessenberg matrix are desired,
    Z must contain the identity matrix.

  On output:

    H has been destroyed.

    Wr and Wi contain the real and imaginary parts, respectively, of
    the eigenvalues. The eigenvalues are unordered except that complex
    conjugate pairs of values appear consecutively with the eigenvalue
    having the positive imaginary part first.

    Z contains the real and imaginary parts of the eigenvectors. If the
    i-th eigenvalue is real, the i-th column of Z contains its eigenvector.
    If the i-th eigenvalue is complex with positive imaginary part, the i-th
    and (i+1)-th columns of Z contain the real and imaginary parts of its
    eigenvector. The eigenvectors are unnormalized. If an error exit is made,
    none of the eigenvectors has been found.

    The function returns an error code:
       zero       for normal return,
       -j         if the limit of 30*N iterations is exhausted
                  while the j-th eigenvalue is being sought
                  (N being the size of the matrix). The eigenvalues
                  should be correct for indices j+1,...,Ubound.
  ----------------------------------------------------------------------
  Note: This is a crude translation. Many of the original goto's have
  been kept !
  ---------------------------------------------------------------------- }

    procedure Cdiv(Ar, Ai, Br, Bi : Float; var Cr, Ci : Float);
    { Complex division, (Cr,Ci) = (Ar,Ai)/(Br,Bi) }
    var
      S, Ars, Ais, Brs, Bis : Float;
    begin
      S := Abs(Br) + Abs(Bi);
      Ars := Ar / S;
      Ais := Ai / S;
      Brs := Br / S;
      Bis := Bi / S;
      S := Sqr(Brs) + Sqr(Bis);
      Cr := (Ars * Brs + Ais * Bis) / S;
      Ci := (Ais * Brs - Ars * Bis) / S;
    end;

  var
    I, J, K, L, M, N, En, Na, Itn, Its, Mp2, Enm2 : Integer;
    P, Q, R, S, T, W, X, Y, Ra, Sa, Vi, Vr, Zz, Norm, Tst1, Tst2 : Float;
    NotLas : Boolean;

  label
    60, 70, 100, 130, 150, 170, 225, 260, 270, 280, 320, 330, 340,
    600, 630, 635, 640, 680, 700, 710, 770, 780, 790, 795, 800;

  begin
    { Store roots isolated by Balance and compute matrix norm }
    K := Lbound;
    Norm := 0.0;
    for I := Lbound to Ubound do
      begin
        for J := K to Ubound do
          Norm := Norm + Abs(H[I,J]);
        K := I;
        if (I < I_low) or (I > I_igh) then
          begin
            Wr[I] := H[I,I];
            Wi[I] := 0.0;
          end;
      end;

    N := Ubound - Lbound + 1;
    Itn := 30 * N;
    En := I_igh;
    T := 0.0;

60: { Search for next eigenvalues }
    if En < I_low then goto 340;
    Its := 0;
    Na := En - 1;
    Enm2 := Na - 1;

70: { Look for single small sub-diagonal element }
    for L := En downto I_low do
      begin
        if L = I_low then goto 100;
        S := Abs(H[L - 1,L - 1]) + Abs(H[L,L]);
        if S = 0.0 then S := Norm;
        Tst1 := S;
        Tst2 := Tst1 + Abs(H[L,L - 1]);
        if Tst2 = Tst1 then goto 100;
      end;

100: { Form shift }
    X := H[En,En];
    if L = En then goto 270;
    Y := H[Na,Na];
    W := H[En,Na] * H[Na,En];
    if L = Na then goto 280;

    if Itn = 0 then
      begin
        { Set error -- all eigenvalues have not
          converged after 30*N iterations }
        Hqr2 := - En;
        Exit;
      end;

    if (Its <> 10) and (Its <> 20) then goto 130;

    { Form exceptional shift }
    T := T + X;

    for I := I_low to En do
      H[I,I] := H[I,I] - X;

    S := Abs(H[En,Na]) + Abs(H[Na,Enm2]);
    X := 0.75 * S;
    Y := X;
    W := - 0.4375 * S * S;

130:
    Its := Its + 1;
    Itn := Itn - 1;

    { Look for two consecutive small sub-diagonal elements }
    for M := Enm2 downto L do
      begin
        Zz := H[M,M];
        R := X - Zz;
        S := Y - Zz;
        P := (R * S - W) / H[M + 1,M] + H[M,M + 1];
        Q := H[M + 1,M + 1] - Zz - R - S;
        R := H[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 150;
        Tst1 := Abs(P) * (Abs(H[M - 1,M - 1]) + Abs(Zz) + Abs(H[M + 1,M + 1]));
        Tst2 := Tst1 + Abs(H[M,M - 1]) * (Abs(Q) + Abs(R));
        if Tst2 = Tst1 then goto 150;
      end;

150:
    Mp2 := M + 2;

    for I := Mp2 to En do
      begin
        H[I,I - 2] := 0.0;
        if I <> Mp2 then H[I,I - 3] := 0.0;