{ **********************************************************************
  *                           Unit EIGEN.PAS                           *
  *                            Version 2.2d                            *
  *                      (c) J. Debord, July 2004                      *
  **********************************************************************
          Procedures for computing eigenvalues and eigenvectors
  **********************************************************************
  References:
  1) 'Mathematiques et Statistiques' by H. Haut (PSI ed.) : Jacobi
  2) EISPACK (http://www.netlib.org/eispack) : EigenVals, EigenVect
  ********************************************************************** }

unit eigen;

interface

uses
  fmath, matrices, stat;

function Jacobi(A : TMatrix; Lbound, Ubound, MaxIter : Integer;
                Tol : Float; V : TMatrix; Lambda : TVector) : Integer;
{ ----------------------------------------------------------------------
  Eigenvalues and eigenvectors of a symmetric matrix by the iterative
  method of Jacobi
  ----------------------------------------------------------------------
  Input parameters  : A       = matrix
                      Lbound  = index of first matrix element
                      Ubound  = index of last matrix element
                      MaxIter = maximum number of iterations
                      Tol     = required precision
  ----------------------------------------------------------------------
  Output parameters : V      = matrix of eigenvectors (stored by columns)
                      Lambda = eigenvalues in decreasing order
  ----------------------------------------------------------------------
  Possible results  : MAT_OK
                      MAT_NON_CONV
  ----------------------------------------------------------------------
  NB : 1. The eigenvectors are normalized, with their first component > 0
       2. This procedure destroys the original matrix A
  ---------------------------------------------------------------------- }

function EigenVals(A : TMatrix; Lbound, Ubound : Integer;
                   Lambda_Re, Lambda_Im : TVector) : Integer;
{ ----------------------------------------------------------------------
  Eigenvalues of a general square matrix
  ----------------------------------------------------------------------
  Input parameters  : A      = matrix
                      Lbound = index of first matrix element
                      Ubound = index of last matrix element
  ----------------------------------------------------------------------
  Output parameters : Lambda_Re = real part of eigenvalues
                      Lambda_Im = imaginary part of eigenvalues

                      The eigenvalues are unordered, except that complex
                      conjugate pairs appear consecutively with the
                      value having the positive imaginary part first.
  ----------------------------------------------------------------------
  Possible results  :  0 : No error
                      -i : Non-convergence has been encountered during
                           the search for the i-th eigenvalue. The
                           eigenvalues should be correct for indices
                           (i+1)..Ubound.
  ----------------------------------------------------------------------
  NB : This procedure destroys the original matrix A
  ---------------------------------------------------------------------- }

function EigenVect(A : TMatrix; Lbound, Ubound : Integer;
                   Lambda_Re, Lambda_Im : TVector; V : TMatrix) : Integer;
{ ----------------------------------------------------------------------
  Eigenvalues and eigenvectors of a general square matrix
  ----------------------------------------------------------------------
  Input parameters  : A      = matrix
                      Lbound = index of first matrix element
                      Ubound = index of last matrix element
  ----------------------------------------------------------------------
  Output parameters : Lambda_Re = real part of eigenvalues
                      Lambda_Im = imaginary part of eigenvalues
                      V         = matrix of eigenvectors

  If the i-th eigenvalue is real, the i-th column of V contains its
  eigenvector. If the i-th eigenvalue is complex with positive imaginary
  part, the i-th and (i+1)-th columns of V contain the real and imaginary
  parts of its eigenvector. The eigenvectors are unnormalized.
  ----------------------------------------------------------------------
  Possible results  :  0 : No error
                      -i : Non-convergence has been encountered during
                           the search for the i-th eigenvalue. None of
                           the eigenvectors has been found. The
                           eigenvalues should be correct for indices
                           (i+1)..Ubound.
  ----------------------------------------------------------------------
  NB : This procedure destroys the original matrix A
  ---------------------------------------------------------------------- }

procedure DivLargest(V : TVector; Lbound, Ubound : Integer;
                     var Largest : Float);
{ ----------------------------------------------------------------------
  Normalizes an eigenvector V by dividing by the element with the
  largest absolute value
  ---------------------------------------------------------------------- }

function RootPol(Coef : TVector; Deg : Integer;
                 Xr, Xi : TVector) : Integer;
{ ----------------------------------------------------------------------
  Real and complex roots of a real polynomial by the method of the
  companion matrix
  ----------------------------------------------------------------------
  Input parameters  : Coef = coefficients of polynomial
                      Deg  = degree of polynomial
  ----------------------------------------------------------------------
  Output parameters : Xr = real parts of root
                      Xi = imaginary parts of root
  ----------------------------------------------------------------------
  If no error occurred, the function returns the number of real roots.
  The N real roots are returned in Xr[1..N] in increasing order. The
  complex roots are returned in Xr[N+1..Deg], Xi[N+1..Deg] and are
  unordered.

  If an error occurred during the search for the i-th root, the function
  returns (-i). The roots should be correct for indices (i+1)..Deg. The
  roots are unordered.
  ---------------------------------------------------------------------- }

implementation

  function Jacobi(A : TMatrix; Lbound, Ubound, MaxIter : Integer;
                  Tol : Float; V : TMatrix; Lambda : TVector) : Integer;
  var
    I, J, K, Im, Jm, Iter : Integer;
    B, C, C2, Na, Nd, P, Q, S, S2, R, T : Float;
  begin
    Iter := 0;
    Na := 0.0;
    Nd := 0.0;
    R := 0.0;

    for I := Lbound to Ubound do
      begin
        V[I,I] := 1.0;
        Nd := Nd + Sqr(A[I,I]);
        if I <> Ubound then
          for J := Succ(I) to Ubound do
            begin
              R := R + Sqr(A[I,J]);
              V[I,J] := 0.0;
              V[J,I] := 0.0;
            end;
      end;

    Na := Nd + 2.0 * R;

    repeat
      R := 0.0;
      for I := Lbound to Pred(Ubound) do
        for J := Succ(I) to Ubound do
          begin
            T := Abs(A[I,J]);
            if T > R then
              begin
                R := T;
                Im := I;
                Jm := J;
              end;
          end;

      B := A[Im,Im] - A[Jm,Jm];

      if B = 0 then
        begin
          C := SQRT2DIV2;
          S := C * Sgn(A[Im,Jm]);
        end
      else
        begin
          P := 2.0 * A[Im,Jm] * Sgn(B);
          Q := Abs(B);
          R := Pythag(P, Q);
          C := Sqrt(0.5 * (1.0 + Q / R));
          S := 0.5 * P / (R * C);
        end;

      for K := Lbound to Ubound do
        begin
          R := V[K,Im];
          V[K,Im] := C * R + S * V[K,Jm];
          V[K,Jm] := C * V[K,Jm] - S * R;
        end;

      if Im <> Lbound then
        for K := Lbound to Pred(Im) do
          begin
            R := A[K,Im];
            A[K,Im] := C * R + S * A[K,Jm];
            A[K,Jm] := C * A[K,Jm] - S * R;
          end;

      if Jm <> Succ(Im) then
        for K := Succ(Im) to Pred(Jm) do
          begin
            R := A[Im,K];
            A[Im,K] := C * R + S * A[K,Jm];
            A[K,Jm] := C * A[K,Jm] - S * R;
          end;

      if Jm <> Ubound then
        for K := Succ(Jm) to Ubound do
          begin
            R := A[Im,K];
            A[Im,K] := C * R + S * A[Jm,K];
            A[Jm,K] := C * A[Jm,K] - S * R;
          end;

      Nd := Nd + 2.0 * Sqr(A[Im,Jm]);

      C2 := Sqr(C);
      S2 := Sqr(S);
      P := 2.0 * S * C * A[Im,Jm];
      R := A[Im,Im];
      A[Im,Im] := C2 * R + S2 * A[Jm,Jm] + P;
      A[Jm,Jm] := S2 * R + C2 * A[Jm,Jm] - P;
      A[Im,Jm] := 0.0;

      Inc(Iter);
      if Iter > MaxIter then
        begin
          Jacobi := MAT_NON_CONV;
          Exit;
        end;
    until Abs(1.0 - Na / Nd) < Tol;

    { The diagonal terms of the transformed matrix are the eigenvalues }
    for I := Lbound to Ubound do
      Lambda[I] := A[I,I];

    { Sort eigenvalues and eigenvectors }
    for I := Lbound to Pred(Ubound) do
      begin
        K := I;
        R := Lambda[I];
        for J := Succ(I) to Ubound do
          if Lambda[J] > R then
            begin
              K := J;
              R := Lambda[J];
            end;
        Swap(Lambda[I], Lambda[K]);
        SwapCols(I, K, V, Lbound, Ubound);
      end;

    { Make sure that the first component of each eigenvector is > 0 }
    for J := Lbound to Ubound do
      if V[Lbound, J] < 0.0 then
        for I := Lbound to Ubound do
          V[I,J] := - V[I,J];

    Jacobi := MAT_OK;
  end;

  procedure Balance(A : TMatrix; Lbound, Ubound : Integer;
                    var I_low, I_igh : Integer; Scale : TVector);
{ ----------------------------------------------------------------------
  This procedure is a translation of the EISPACK procedure Balanc.

  This procedure balances a real matrix and isolates eigenvalues
  whenever possible.

  On input:

    A contains the input matrix to be balanced.

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

  On output:

    A contains the balanced matrix.

    I_low and I_igh are two integers such that A[i,j]
    is equal to zero if
      (1) i is greater than j and
      (2) j=Lbound,...,I_low-1 or i=I_igh+1,...,Ubound.

    Scale contains information determining the permutations
    and scaling factors used.

    Suppose that the principal submatrix in rows I_low through I_igh
    has been balanced, that P[j] denotes the index interchanged
    with j during the permutation step, and that the elements
    of the diagonal matrix used are denoted by D[i,j].  then
        Scale[j] = P[j],    for j = Lbound,...,I_low-1
                 = D[j,j],      j = I_low,...,I_igh
                 = P[j]         j = I_igh+1,...,Ubound.
    the order in which the interchanges are made is
    Ubound to I_igh+1, then Lbound to I_low-1.

    Note that Lbound is returned for I_igh if I_igh is < Lbound formally
  ---------------------------------------------------------------------- }
  const
    RADIX = 2;  { Base used in floating number representation }

  var
    I, J, M : Integer;
    C, F, G, R, S, B2 : Float;
    Flag, Found, Conv : Boolean;

    procedure Exchange;
    { Row and column exchange }
    var
      I : Integer;
    begin
      Scale[M] := J;
      if J = M then Exit;

      for I := Lbound to I_igh do
        Swap(A[I,J], A[I,M]);

      for I := I_low to Ubound do
        Swap(A[J,I], A[M,I]);
    end;

  begin
    B2 := RADIX * RADIX;
    I_low := Lbound;
    I_igh := Ubound;

    { Search for rows isolating an eigenvalue and push them down }
    repeat
      J := I_igh;
      repeat
        I := Lbound;
        repeat
          Flag := (I <> J) and (A[J,I] <> 0.0);
          I := I + 1;
        until Flag or (I > I_igh);
        Found := not Flag;
        if Found then
          begin
            M := I_igh;
            Exchange;
            I_igh := I_igh - 1;
          end;
        J := J - 1;
      until Found or (J < Lbound);
    until (not Found) or (I_igh < Lbound);

    if I_igh < Lbound then I_igh := Lbound;
    if I_igh = Lbound then Exit;

    { Search for columns isolating an eigenvalue and push them left }
    repeat
      J := I_low;
      repeat
        I := I_low;
        repeat
          Flag := (I <> J) and (A[I,J] <> 0.0);
          I := I + 1;
        until Flag or (I > I_igh);
        Found := not Flag;
        if Found then
          begin
            M := I_low;
            Exchange;
            I_low := I_low + 1;
          end;
        J := J + 1;
      until Found or (J > I_igh);
    until (not Found);

    { Now balance the submatrix in rows I_low to I_igh }
    for I := I_low to I_igh do
      Scale[I] := 1.0;

    { Iterative loop for norm reduction }
    repeat
      Conv := True;

      for I := I_low to I_igh do
        begin
          C := 0.0;
          R := 0.0;

          for J := I_low to I_igh do
            if J <> I then
              begin
                C := C + Abs(A[J,I]);
                R := R + Abs(A[I,J]);
              end;

          { Guard against zero C or R due to underflow }
          if (C <> 0.0) and (R <> 0.0) then
            begin
              G := R / RADIX;
              F := 1.0;
              S := C + R;

              while C < G do
                begin
                  F := F * RADIX;
                  C := C * B2;
                end;

              G := R * RADIX;

              while C >= G do
                begin
                  F := F / RADIX;
                  C := C / B2;
                end;

              { Now balance }
              if (C + R) / F < 0.95 * S then
                begin
                  G := 1.0 / F;
                  Scale[I] := Scale[I] * F;
                  Conv := False;
                  for J := I_low to Ubound do A[I,J] := A[I,J] * G;
                  for J := Lbound to I_igh do A[J,I] := A[J,I] * F;
                end;
            end;
        end;
    until Conv;
  end;

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

  Given a real general matrix, this procedure reduces a submatrix
  situated in rows and columns I_low through I_igh to upper Hessenberg
  form by stabilized elementary similarity transformations.

  On input:

    A contains the input matrix.

    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.

  On output:

    A contains the Hessenberg matrix. The multipliers which were used
    in the reduction are stored in the remaining triangle under the
    Hessenberg matrix.

    I_int contains information on the rows and columns interchanged in
    the reduction. Only elements I_low through I_igh are used.
  ---------------------------------------------------------------------- }
  var
    I, J, M, La, Kp1, Mm1, Mp1 : Integer;
    X, Y : Float;
  begin
    La := I_igh - 1;
    Kp1 := I_low + 1;
    if La < Kp1 then Exit;

    for M := Kp1 to La do
      begin
        Mm1 := M - 1;
        X := 0.0;
        I := M;

        for J := M to I_igh do
          if Abs(A[J,Mm1]) > Abs(X) then
            begin
              X := A[J,Mm1];
              I := J;
            end;

        I_int[M] := I;

        { Interchange rows and columns of A }
        if I <> M then
          begin
            for J := Mm1 to Ubound do
              Swap(A[I,J], A[M,J]);

            for J := Lbound to I_igh do
              Swap(A[J,I], A[J,M]);
          end;

        if X <> 0.0 then
          begin
            Mp1 := M + 1;
            for I := Mp1 to I_igh do
              begin
                Y := A[I,Mm1];
                if Y <> 0.0 then