-- This package is part of the test suite for Eigenvalue / Eigenvector
-- testing


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 COMPOSE_FROM_CARTESIAN(L,R:REAL) return COMPLEX is <>;
  with function RE(L:COMPLEX) return REAL is <>;
  with function IM(L:COMPLEX) return REAL is <>;
  with function "+"(L,R:COMPLEX) return COMPLEX is <>;  
  with function "-"(L,R:COMPLEX) return COMPLEX is <>;  
  with function "*"(L,R:COMPLEX) return COMPLEX is <>;  
  with function "/"(L,R:COMPLEX) return COMPLEX is <>;  
  with function "-"(L,R:COMPLEX_VECTOR) return COMPLEX_VECTOR is <>;  
  with function "*"(L:COMPLEX;
                    R:COMPLEX_VECTOR) return COMPLEX_VECTOR is <>;
  with function "*"(L:COMPLEX_MATRIX;
                    R:COMPLEX_VECTOR) return COMPLEX_VECTOR is <>;  
  with function "+" (LEFT : REAL ; RIGHT : IMAGINARY)
                     return COMPLEX is <>;
  with function "*" (LEFT : REAL; RIGHT : IMAGINARY)
                     return IMAGINARY is <>;

package GENERIC_COMPLEX_EIGEN_CHECK is


  procedure EIGDET ( B : COMPLEX_MATRIX ;
                     LAMBDA : COMPLEX_VECTOR ) ;

  procedure EVALVEC ( AA : COMPLEX_MATRIX ;
                      LAMBDA : COMPLEX_VECTOR ;
                      X : COMPLEX_MATRIX ) ;
end GENERIC_COMPLEX_EIGEN_CHECK ;

with TEXT_IO ; use TEXT_IO ;
with GENERIC_INTEGER_ARRAYS;
with COMPLEX_IO;
package body GENERIC_COMPLEX_EIGEN_CHECK is

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


  package INT_IO is new INTEGER_IO(INTEGER);
  use INT_IO;

  package FLT_IO is new FLOAT_IO(REAL);
  use FLT_IO;

  package CX_IO is new COMPLEX_IO( REAL , COMPLEX ) ;
  use CX_IO;

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



  procedure EIGDET ( B : COMPLEX_MATRIX ;
                     LAMBDA : COMPLEX_VECTOR ) is

--      PURPOSE : COMPUTE DETERMINANT OF A N BY N COMPLEX MATRIX 
--
--      INPUT :  B - THE COMPLEX MATRIX TO BE INVERTED
--               N - THE NUMBER OF COLUMNS ( NUMBER OF ROWS) USED
--               M - DIMENSION OF B
--
--      OUTPUT : DET - THE DETERMINANT OF THE ORIGIONAL MATRIX
--               SING - DEGREE OF SINGULARITY ( NOT IN DET )
--
--      METHOD : GAUSS / JORDAN  ELIMINATION USING MAXIMUM
--               PIVOTAL ELEMENT. USE OF DIVIDE IS MINIMIZED.
--              UNSCRAMBLING IS EFFICIENT AT THE EXPENSE OF
--               EXTRA STORAGE. DETERMINANT MAY HAVE WRONG SIGN.


    N : INTEGER := B'LENGTH ( 1 ) ;
    A : COMPLEX_MATRIX ( 1 .. N , 1 .. N ) ;
    IL : INTEGER ;
    DET : COMPLEX := COMPOSE_FROM_CARTESIAN( 1.0 , 0.0 ) ;
    LARGE : COMPLEX ;
    ROW : INTEGER_VECTOR ( 1 .. N ) ;
    TEMP : COMPLEX_VECTOR ( 1 .. N ) ;
    MTEMP : INTEGER ;
    SING : INTEGER ;
  begin

--            COPY B INTO A   SUBTRACTING LAMBDA FROM DIAGONAL
    for ILAMBDA in LAMBDA'RANGE loop
      for I in 1 .. N loop
        for J in 1 .. N loop
          A ( I , J ) := B ( I , J ) ;
          if I = J then
            A ( I , J ) := A ( I , J ) - LAMBDA ( ILAMBDA ) ;
          end if ;
        end loop ;
      end loop ;
      SING := 0 ;
      for K in 1 .. N loop
        ROW ( K ) := K ;
      end loop ;

--                BEGIN MAIN REDUCTION LOOP ON  K
      for K in 1 .. N loop

--          FIND LARGEST
        LARGE := A ( ROW( K ) , K) ;
        IL := K ;
        for I in K + 1 .. N loop
          if SUMABS ( A( ROW( I ) , K)) > SUMABS ( LARGE ) then
            IL := I ;
            LARGE := A ( ROW( I ) , K) ;
          end if ;
        end loop ;
        if IL /= K then
          MTEMP := ROW ( K ) ;
          ROW ( K ) := ROW ( IL ) ;
          ROW ( IL ) := MTEMP ;
        end if ;

--                        END MAX INTERCHANGE
--                     LARGE IS NOW THE  A(K,K) PIVOT
        if SUMABS ( LARGE ) < REAL'EPSILON then
          SING := SING + 1 ;
        else
          DET := DET * LARGE ;
          for J in K + 1 .. N loop
            A ( ROW( K ) , J) := A ( ROW( K ) , J) / A ( ROW( K ) , K) ;
          end loop ;
          for I in 1 .. N loop
            if I /= K then
              for J in K + 1 .. N loop
                A ( ROW( I ) , J) := A ( ROW( I ) , J) - A ( ROW( I ) , K) * A
                   ( ROW( K ) , J) ;
              end loop ;
            end if ;
          end loop ;
        end if ;
      end loop ;

--                    END OF MAIN REDUCTION LOOP ON K
      if SING > 0 then
        PUT ( SING ) ;
        PUT_LINE ( " = SINGULARITY" ) ;
      else
        if SUMABS ( DET ) < 2.0 * REAL'EPSILON then
          SING := 1 ;
          PUT ( SING ) ;
          PUT_LINE ( " = SINGULARITY" ) ;
        else
          PUT ( DET ) ;
          PUT_LINE ( " =  DETERMINANT " ) ;
        end if ;
      end if ;
    end loop ;
  exception
    when NUMERIC_ERROR | CONSTRAINT_ERROR =>
      PUT_LINE ( " NUMERIC ERROR IN EIGDET " ) ;
  end EIGDET ;

  procedure EVALVEC ( AA : COMPLEX_MATRIX ;
                      LAMBDA : COMPLEX_VECTOR ;
                      X : COMPLEX_MATRIX ) is

--             EIGENVALUE EIGENVECTOR ACCURACY CHECK ROUTINE
--             AA      - ORIGONAL COMPLEX MATRIX
--             LAMBDA - COMPLEX VECTOR OF EIGENVALUES
--             X      - COMPLEX MATRIX CONTAINING EIGENVECTORS
    N : INTEGER := AA'LENGTH ;
    T : COMPLEX_VECTOR ( 1 .. N ) ;
    T1 : COMPLEX_VECTOR ( 1 .. N ) ;
    T2 : COMPLEX_VECTOR ( 1 .. N ) ;
    D1 : COMPLEX_VECTOR ( 1 .. N ) ;
    D2 : COMPLEX_VECTOR ( 1 .. N ) ;
    XV : COMPLEX_VECTOR ( X'FIRST( 1 ) .. X'LAST ( 1 )) ;
    SUM : REAL ;
    TOTAL : REAL ;
  begin

--          CHECK ACCURACY OF  det | A - IDENT * LAMBDA(I) | = 0   I=1..N
    EIGDET ( AA , LAMBDA ) ;

--          CHECK ACCURACY OF  X A = LAMBDA X     1..Nth EIGENVECTOR 
    TOTAL := 0.0 ;
    for I in X'RANGE ( 2 ) loop
      for J in X'RANGE ( 1 ) loop
        XV ( J ) := X ( J , I ) ;
      end loop ;
      T1 := AA * XV ;
      T := LAMBDA ( I ) * XV ;
      D1 := T - T1 ;
      SUM := 0.0 ;
      for J in 1 .. N loop
        SUM := SUM + SUMABS ( D1( J )) ;
      end loop ;
      PUT ( SUM ) ;
      PUT_LINE ( " = SUM" ) ;
      TOTAL := TOTAL + SUM ;
    end loop ;
    PUT ( TOTAL ) ;
    PUT_LINE ( " = TOTAL EIGENVECTOR ERROR " ) ;
  end EVALVEC ;
end GENERIC_COMPLEX_EIGEN_CHECK ;