-- Generic_Complex_Eigenvalues

-- The procedure EIGEN computes both the eigenvalues and eigenvectors
-- of arbitrary complex matricies. The generic formal parameter REAL
-- determines the type of each component of the type COMPLEX.
-- This is an Ada version of the NAG Fortran library subroutine
-- Some Fortran labels have been preserved for traceability

generic

  type REAL is digits <> ;
  type COMPLEX is private;
  type IMAGINARY is private;
  type COMPLEX_VECTOR is array(INTEGER range <>) of COMPLEX;
  type COMPLEX_MATRIX is array(INTEGER range <>,INTEGER range <>) of COMPLEX;
  with function  RE (X : COMPLEX) return REAL is <>;
  with function  IM (X : COMPLEX) return REAL is <>;
  with function  COMPOSE_FROM_CARTESIAN (RE, IM : REAL) return COMPLEX is <>;
  with function  "abs" (RIGHT : COMPLEX) return REAL is <>;
  with function  "-" (RIGHT : COMPLEX) return COMPLEX is <>;
  with function  "+" (LEFT, RIGHT : COMPLEX) return COMPLEX is <>;
  with function  "-" (LEFT, RIGHT : COMPLEX) return COMPLEX is <>;
  with function  "*" (LEFT, RIGHT : COMPLEX) return COMPLEX is <>;
  with function  "/" (LEFT, RIGHT : COMPLEX) return COMPLEX is <>;
  with function IDENTITY_MATRIX(I,J,K:INTEGER) return COMPLEX_MATRIX is <>;
  with function SQRT(X : COMPLEX) return COMPLEX is <>;
package GENERIC_COMPLEX_EIGENVALUES is

  procedure EIGEN ( A : in out COMPLEX_MATRIX ;
                    VALUES : in out COMPLEX_VECTOR ;
                    VECTORS : in out COMPLEX_MATRIX ;
                    FAIL : in out BOOLEAN ) ;

  function HESS ( A : COMPLEX_MATRIX ) return COMPLEX_MATRIX;

end GENERIC_COMPLEX_EIGENVALUES ;

with TEXT_IO; use TEXT_IO; --debug
with GENERIC_INTEGER_ARRAYS;
package body GENERIC_COMPLEX_EIGENVALUES is

  package INTEGER_ARRAYS is new GENERIC_INTEGER_ARRAYS(INTEGER);
  use INTEGER_ARRAYS;

  procedure CXHESS ( A : in out COMPLEX_MATRIX ;
                     INTERCHANGE : in out INTEGER_VECTOR ) ;

  procedure CXEIG2C ( A : in out COMPLEX_MATRIX ;
                      VALUES : in out COMPLEX_VECTOR ;
                      VECTORS : in out COMPLEX_MATRIX ;
                      INTERCHANGE : in out INTEGER_VECTOR ;
                      FAIL : in out BOOLEAN ) ;


  function SUMABS ( Z : COMPLEX ) return REAL is
  begin
    return abs(RE(Z))+abs(IM(Z));
  end SUMABS ;

  procedure EIGEN ( A : in out COMPLEX_MATRIX ;
                    VALUES : in out COMPLEX_VECTOR ;
                    VECTORS : in out COMPLEX_MATRIX ;
                    FAIL : in out BOOLEAN ) is

--  DRIVER PROCEDURE FOR EIGENVALUES AND EIGEN VECTORS OF COMPLEX MATRICIES
    INTERCHANGE : INTEGER_VECTOR ( A'FIRST .. A'LAST ) ;
  begin
    FAIL := FALSE ;
    CXHESS ( A , INTERCHANGE ) ;
    CXEIG2C ( A , VALUES , VECTORS , INTERCHANGE , FAIL ) ;
  end EIGEN ;

  procedure CXHESS ( A : in out COMPLEX_MATRIX ;
                     INTERCHANGE : in out INTEGER_VECTOR ) is
    K , I , N : INTEGER ;
    X : COMPLEX ;
    Y : COMPLEX ;
  begin
    K := A'FIRST ( 1 ) ;
    N := A'LAST ( 1 ) ;
    for M in K + 1 .. N - 1 loop
      I := M ;
      X := COMPOSE_FROM_CARTESIAN( 0.0 , 0.0 ) ;
      for J in M .. N loop
        if SUMABS ( A( J , M - 1 )) > SUMABS ( X ) then
          X := A ( J , M - 1 ) ;
          I := J ;
        end if ;
      end loop ;
      INTERCHANGE ( M ) := I ;
      if I /= M then

--              INTERCHANGE ROWS AND COLUMNS OF A
        for J in M - 1 .. N loop
          Y := A ( I , J ) ;
          A ( I , J ) := A ( M , J ) ;
          A ( M , J ) := Y ;
        end loop ;
        for J in 1 .. N loop
          Y := A ( J , I ) ;
          A ( J , I ) := A ( J , M ) ;
          A ( J , M ) := Y ;
        end loop ;
      end if ;
      if SUMABS ( X ) /= 0.0 then
        for I in M + 1 .. N loop
          Y := A ( I , M - 1 ) ;
          if SUMABS ( Y ) > 0.0 then
            Y := Y / X ;
            A ( I , M - 1 ) := Y ;
            for J in M .. N loop
              A ( I , J ) := A ( I , J ) - Y * A ( M , J ) ;
            end loop ;
            for J in 1 .. N loop
              A ( J , M ) := A ( J , M ) + Y * A ( J , I ) ;
            end loop ;
          end if ;
        end loop ;
      end if ;
    end loop ;
  end CXHESS ;

  procedure CXEIG2C ( A : in out COMPLEX_MATRIX ;
                      VALUES : in out COMPLEX_VECTOR ;
                      VECTORS : in out COMPLEX_MATRIX ;
                      INTERCHANGE : in out INTEGER_VECTOR ;
                      FAIL : in out BOOLEAN ) is
    N : INTEGER ;
    J : INTEGER ;
    M : INTEGER ;
    K : INTEGER ;
    MM : INTEGER ;
    LOW : INTEGER ;
    ITS : INTEGER ;
    ITN : INTEGER ;
    IEN : INTEGER ;
    ANORM : REAL := 0.0 ;
    AHR : REAL ;
    AAHR : REAL ;
    ACC : REAL ;
    XR : REAL ;
    XI : REAL ;
    YR : REAL ;
    YI : REAL ;
    ZR : REAL ;
    ACCNORM : COMPLEX ;
    X : COMPLEX ;
    Y : COMPLEX ;
    Z : COMPLEX ;
    YY : COMPLEX ;
    T : COMPLEX ;
    S : COMPLEX ;
  begin
    LOW := A'FIRST ;
    N := A'LAST ;
    ACC := 2.0 ** ( - 23 ) ;
    T := COMPOSE_FROM_CARTESIAN( 0.0 , 0.0 ) ;
    ITN := 30 * A'LENGTH ;
    VECTORS := IDENTITY_MATRIX ( VECTORS'LENGTH, 1, 1 ) ;

--          GET EIGENVALUE TRANSFORMATION FROM HESSENBERG REDUCTION  LEFT IN H
    for I in reverse A'FIRST + 1 .. A'LAST - 1 loop
      J := INTERCHANGE ( I ) ;
      for K in I + 1 .. A'LAST loop
        VECTORS ( K , I ) := A ( K , I - 1 ) ;
      end loop ;
      if I /= J then
        for K in I .. A'LAST loop
          VECTORS ( I , K ) := VECTORS ( J , K ) ;
          VECTORS ( J , K ) := COMPOSE_FROM_CARTESIAN( 0.0 , 0.0 ) ;
        end loop ;
        VECTORS ( J , I ) := COMPOSE_FROM_CARTESIAN( 1.0 , 0.0 ) ;
      end if ;
    end loop ;
    IEN := N ;

--         IEN IS DECREMENTED
    while LOW <= IEN loop  -- 260
      ITS := 0 ;

--     LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
      L280 : -- NO, THIS IS NOT FOR A GOTO, IT IS A LOOP LABEL FOR AN EXIT
      loop  --  280 
        K := LOW ;
        for KK in reverse LOW + 1 .. IEN loop  -- 300
          AHR := SUMABS ( A( KK , KK - 1 )) ;
          AAHR := ACC * ( SUMABS( A( KK - 1 , KK - 1 )) + SUMABS ( A( KK , KK
             ))) ;
          if AHR <= AAHR then
            K := KK ;
            exit ;
          end if ;
        end loop ;  --  300 
        exit L280 when K = IEN ;  -- 780
        if ITN <= 0 then
          FAIL := TRUE ;
          return ;
        end if ;

--     COMPUTE SHIFT
        if ITS = 10 or
           ITS = 20 then
          S := COMPOSE_FROM_CARTESIAN ( abs( RE( A( IEN , IEN - 1 ))) +
               abs ( RE( A( IEN - 1 , IEN - 2 ))) ,
               abs ( IM( A( IEN , IEN - 1 ))) +
               abs ( IM( A( IEN - 1 , IEN - 2 )))) ;
        else
          S := A ( IEN , IEN ) ;
          X := A ( IEN - 1 , IEN ) * A ( IEN , IEN - 1 ) ;
          if SUMABS ( X ) > 0.0 then
            Y := ( A( IEN - 1 , IEN - 1 ) - S) /COMPOSE_FROM_CARTESIAN(2.0,0.0);
            Z := SQRT ( Y * Y + X ) ;
            if RE ( Y ) * RE ( Z ) + IM ( Y ) *
               IM ( Z ) < 0.0 then
              Z := - Z ;
            end if ;
            YY := Y + Z ;
            S := S - X / YY ;
          end if ;
        end if ;  --  400 
        for I in LOW .. IEN loop  -- 420
          A ( I , I ) := A ( I , I ) - S ;
        end loop ;  --  420 
        T := T + S ;
        ITS := ITS + 1 ;
        ITN := ITN - 1 ;
        J := K + 1 ;

--     LOOK FOR TWO CONSECUTIVE SMALL SUB-DIAGONAL ELEMENTS
        XR := SUMABS ( A( IEN - 1 , IEN - 1 )) ;
        YR := SUMABS ( A( IEN , IEN - 1 )) ;
        ZR := SUMABS ( A( IEN , IEN )) ;
        M := K ;
        for MM in reverse J .. IEN - 1 loop  -- 460
          YI := YR ;
          YR := SUMABS ( A( MM , MM - 1 )) ;
          XI := ZR ;
          ZR := XR ;
          XR := SUMABS ( A( MM - 1 , MM - 1 )) ;
          if YR <= ( ACC * ZR / YI *( ZR + XR + XI )) then
            M := MM ;
            exit ;
          end if ;
        end loop ;  --  460 

--     TRIANGULAR DECOMPOSITION A := L*R
        for I in M + 1 .. IEN loop  -- 620
          X := A ( I - 1 , I - 1 ) ;
          Y := A ( I , I - 1 ) ;
          if SUMABS ( X ) >= SUMABS ( Y ) then
            Z := Y / X ;
            VALUES ( I ) := COMPOSE_FROM_CARTESIAN( - 1.0 , 0.0 ) ;
          else

--     INTERCHANGE ROWS OF H
            for J in I - 1 .. N loop  -- 540
              Z := A ( I - 1 , J ) ;
              A ( I - 1 , J ) := A ( I , J ) ;
              A ( I , J ) := Z ;
            end loop ;  --  540
            Z := X / Y ;
            VALUES ( I ) := COMPOSE_FROM_CARTESIAN( 1.0 , 0.0 ) ;
          end if ;
          A ( I , I - 1 ) := Z ;
          for J in I .. N loop  -- 600
            A ( I , J ) := A ( I , J ) - Z * A ( I - 1 , J ) ;
          end loop ;  --  600
        end loop ;  --  620 

--     COMPOSITION R*L := H
        for J in M + 1 .. IEN loop  -- 760
          X := A ( J , J - 1 ) ;
          A ( J , J - 1 ) := COMPOSE_FROM_CARTESIAN( 0.0 , 0.0 ) ;

--     INTERCHANGE COLUMNS OF A AND VECTORS , IF NECESSARY
          if RE ( VALUES( J )) > 0.0 then
            for I in LOW .. J loop  -- 660
              Z := A ( I , J - 1 ) ;
              A ( I , J - 1 ) := A ( I , J ) ;
              A ( I , J ) := Z ;
            end loop ;  --  660
            for I in LOW .. N loop  -- 680
              Z := VECTORS ( I , J - 1 ) ;
              VECTORS ( I , J - 1 ) := VECTORS ( I , J ) ;
              VECTORS ( I , J ) := Z ;
            end loop ;  --  680
          end if ;

--     END INTERCHANGE COLUMNS
          for I in LOW .. J loop  -- 720
            A ( I , J - 1 ) := A ( I , J - 1 ) + X * A ( I , J ) ;
          end loop ;  --  720
          for I in LOW .. N loop  -- 740
            VECTORS ( I , J - 1 ) := VECTORS ( I , J - 1 ) + X * VECTORS ( I ,
               J ) ;
          end loop ;  --  740

--     END ACCUMULATE TRANSFORMATIONS
        end loop ;  --  760 
      end loop L280 ;  -- 280

--     A ROOT FOUND
      VALUES ( IEN ) := A ( IEN , IEN ) + T ;
      IEN := IEN - 1 ;
    end loop ;  -- 260 WHILE

--     ALL ROOTS FOUND
    for I in A'RANGE loop
      ANORM := ANORM + SUMABS ( VALUES( I )) ;
      for J in I + 1 .. A'LAST loop
        ANORM := ANORM + SUMABS ( A( I , J )) ;
      end loop ;
    end loop ;
    ACCNORM := COMPOSE_FROM_CARTESIAN ( ANORM * 2.0 **( - 23 ) , 0.0) ;
    if ANORM = 0.0 or
       N < 2 then
      return ;  -- done
    end if ;

-- BACKSUBSTITUTE TO SET UP VECTORS OF UPPER TRIANGULAR FORM
    for IEN in reverse LOW + 1 .. N loop
      X := VALUES ( IEN ) ;
      for I in reverse LOW .. IEN - 1 loop
        Z := A ( I , IEN ) ;
        for J in I + 1 .. IEN - 1 loop
          Z := Z + A ( I , J ) * A ( J , IEN ) ;
        end loop ;
        Y := X - VALUES ( I ) ;
        if SUMABS ( Y ) = 0.0 then
          Y := ACCNORM ;
        end if ;
        A ( I , IEN ) := Z / Y ;
      end loop ;
    end loop ;

--  MULTIPLY BY TRANSFORMATION MATRIX TO GIVE VECTORS OF ORIGINAL FULL MATRIX
    for J in reverse A'RANGE loop
      for I in A'RANGE loop
        Z := VECTORS ( I , J ) ;
        for K in A'FIRST .. J - 1 loop
          Z := Z + VECTORS ( I , K ) * A ( K , J ) ;
        end loop ;
        VECTORS ( I , J ) := Z ;
      end loop ;
    end loop ;
  end CXEIG2C ;

  function HESS ( A : COMPLEX_MATRIX ) return COMPLEX_MATRIX is
    B : COMPLEX_MATRIX(1..A'LENGTH(1), 1..A'LENGTH(2));
    INTERCHANGE : INTEGER_VECTOR(1..A'LENGTH(1));
    N : INTEGER := A'LENGTH(1);
  begin
    if A'LENGTH(1) /= A'LENGTH(2) then
      raise ARRAY_INDEX_ERROR;
    end if;
    B := A;
    CXHESS ( B, INTERCHANGE ) ;
    put("CXHESS inter ");
    for I in 1..A'LENGTH(1) loop
      put(integer'image(INTERCHANGE(I)) & ' ');
    end loop;
    new_line;
    return B;
  end HESS;

end GENERIC_COMPLEX_EIGENVALUES ;