-- Generic_Complex_Linear_Equations

-- this should be instantiated where an instantiation of the
-- COMPLEX_TYPES and COMPLEX_ARRAYS is visible
--
--   package YOUR_NAME is new GENERIC_COMPLEX_LINEAR_EQUATIONS
--     ( REAL, COMPLEX, COMPLEX_VECTOR, COMPLEX_MATRIX ) ;
--
--   all the rest are defaults

with ARRAY_EXCEPTIONS;

generic

  type REAL is digits <>;
  type COMPLEX 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 "-"(L: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) return COMPLEX is <>;
  with function "abs"(L:COMPLEX) return REAL is <>;

package GENERIC_COMPLEX_LINEAR_EQUATIONS is


  function LINEAR_EQUATIONS ( A : COMPLEX_MATRIX ;
                              Y : COMPLEX_VECTOR ) return COMPLEX_VECTOR ;

  function CROUT_SOLUTION ( A : COMPLEX_MATRIX ;
                            Y : COMPLEX_VECTOR ) return COMPLEX_VECTOR ;

  function DETERMINANT ( A : COMPLEX_MATRIX ) return COMPLEX ;

  function INVERSE ( A : COMPLEX_MATRIX ) return COMPLEX_MATRIX ;

  ARRAY_INDEX_ERROR : exception renames ARRAY_EXCEPTIONS.ARRAY_INDEX_ERROR;

  --
  --  function LINEAR_EQUATIONS ( A : COMPLEX_MATRIX ;
  --                              Y : COMPLEX_VECTOR ) return COMPLEX_VECTOR ;
  --
  --  function CROUT_SOLUTION ( A : COMPLEX_MATRIX ;
  --                            Y : COMPLEX_VECTOR ) return COMPLEX_VECTOR ;
  --
  --      PURPOSE : SOLVE THE LINEAR SYSTEM OF EQUATIONS WITH COMPLEX
  --                COEFFICIENTS   [A] * |X| = |Y|
  --
  --      INPUT  : THE COMPLEX MATRIX  A
  --               THE COMPLEX VECTOR  Y
  --
  --      OUTPUT : THE COMPLEX VECTOR  X
  --
  --      METHOD1 : GUASS-JORDAN ELIMINATION USING MAXIMUM ELEMENT
  --                FOR PIVOT.
  --      METHOD2 : CROUT REDUCTION USING MAXIMIM ELEMENT FOR
  --                DIAGONAL
  --
  --      USAGE  :     X := LINEAR_EQUATIONS ( A , Y ) ;
  --                   X := CROUT_SOLUTION ( A , Y ) ;
  --
  --      EXCEPTION : ARRAY_INDEX_ERROR  RASIED IF RESULT CAN NOT BE COMPUTED
  --
  --  function DETERMINANT ( A : COMPLEX_MATRIX ) return COMPLEX ;
  --
  --      PURPOSE : COMPUTE THE DETERMINANT OF A COMPLEX MATRIX
  --                Z := |A|
  --
  --      INPUT  : THE COMPLEX MATRIX  A   ( UNCHANGED )
  --
  --      OUTPUT : THE COMPLEX NUMBER Z 
  --
  --      METHOD : GUASS-JORDAN ELIMINATION USING MAXIMUM ELEMENT
  --                FOR PIVOT.
  --
  --      USAGE  :     Z := DETERMINANT ( A ) ;
  --
  --      EXCEPTION : ARRAY_INDEX_ERROR  RASIED IF RESULT CAN NOT BE COMPUTED
  --
  --  function INVERSE ( A : COMPLEX_MATRIX ) return COMPLEX_MATRIX ;
  --
  --      PURPOSE : INVERT AN N BY N  COMPLEX MATRIX
  --
  --      INPUT : THE COMPLEX MATRIX  A  ( NOT CHANGED )
  --
  --      OUTPUT : THE INVERSE OF MATRIX  A
  --
  --      METHOD : GUASS-JORDAN ELIMINATION USING MAXIMUM ELEMENT
  --               FOR PIVOT.
  --
  --      EXCEPTION : ARRAY_INDEX_ERROR RAISED IF INVERSE CAN NOT BE COMPUTED
  --
  --      SAMPLE USE :  NEW_MATRIX := INVERSE ( OLD_MATRIX ) ;
  --                    MATRIX := INVERSE ( MATRIX ) * ANOTHER_MATRIX ;
  --
end GENERIC_COMPLEX_LINEAR_EQUATIONS ;

package body GENERIC_COMPLEX_LINEAR_EQUATIONS is

  function LINEAR_EQUATIONS ( A : COMPLEX_MATRIX ;
                              Y : COMPLEX_VECTOR ) return COMPLEX_VECTOR is

  --      PURPOSE : SOLVE THE LINEAR SYSTEM OF EQUATIONS WITH COMPLEX
  --                COEFFICIENTS   [A] * |X| = |Y|
  --
  --      INPUT  : THE COMPLEX MATRIX  A
  --               THE COMPLEX VECTOR  Y
  --
  --      OUTPUT : THE COMPLEX VECTOR  X
  --
  --      METHOD : GUASS-JORDAN ELIMINATION USING MAXIMUM ELEMENT
  --               FOR PIVOT.
  --
  --      USAGE  :     X := LINEAR_EQUATIONS ( A , Y ) ;
  --
  --
  --      EXCEPTION : ARRAY_INDEX_ERROR  RASIED IF RESULT CAN NOT BE COMPUTED
  --
  --    WRITTEN BY : JON SQUIRE , 28 MAY 1983
  --                 updated to new complex spec 11/30/88
  --
  --
  --                          LOCAL SPECIFICATIONS

    N : INTEGER := A'LENGTH ;  -- NUMBER OF EQUATIONS
    X : COMPLEX_VECTOR ( 1 .. N ) ;  -- RESULT BEING COMPUTED
    B : COMPLEX_MATRIX ( 1 .. N , 1 .. N + 1 ) ;  -- WORKING MATRIX
    ROW : array ( 1 .. N ) of INTEGER ;  -- ROW INTERCHANGE INDICIES
    HOLD , I_PIVOT : INTEGER ;  -- PIVOT INDICIES
    PIVOT : COMPLEX ;  -- PIVOT ELEMENT VALUE
    ABS_PIVOT : REAL ;
  begin
    if A'LENGTH ( 1 ) /= A'LENGTH ( 2 ) then
      raise ARRAY_INDEX_ERROR ;
    end if ;

    --                        BUILD WORKING DATA STRUCTURE
    for I in 1 .. N loop
      for J in 1 .. N loop
        B ( I , J ) := A ( I - 1 + A'FIRST( 1 ) , J - 1 + A'FIRST ( 2 )) ;
      end loop ;
      B ( I , N + 1 ) := Y ( I - 1 + Y'FIRST ) ;
    end loop ;

    --                        SET UP ROW  INTERCHANGE VECTORS
    for K in 1 .. N loop
      ROW ( K ) := K ;
    end loop ;

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

      --                             FIND LARGEST ELEMENT FOR PIVOT
      PIVOT := B ( ROW( K ) , K) ;
      ABS_PIVOT := abs ( PIVOT ) ;
      I_PIVOT := K ;
      for I in K .. N loop
        if abs ( B( ROW( I ) , K)) > ABS_PIVOT then
          I_PIVOT := I ;
          PIVOT := B ( ROW( I ) , K) ;
          ABS_PIVOT := abs ( PIVOT ) ;
        end if ;
      end loop ;

      --                                HAVE PIVOT, INTERCHANGE ROW POINTERS
      HOLD := ROW ( K ) ;
      ROW ( K ) := ROW ( I_PIVOT ) ;
      ROW ( I_PIVOT ) := HOLD ;

      --                                CHECK FOR NEAR SINGULAR
      if ABS_PIVOT < 1.0E-5 then
        for J in K + 1 .. N + 1 loop
          B ( ROW( K ) , J) := COMPOSE_FROM_CARTESIAN( 0.0 , 0.0 ) ;
        end loop ;
      else

        --                                REDUCE ABOUT PIVOT
        for J in K + 1 .. N + 1 loop
          B ( ROW( K ) , J) := B ( ROW( K ) , J) / B ( ROW( K ) , K) ;
        end loop ;

        --                               INNER REDUCTION LOOP
        for I in 1 .. N loop
          if I /= K then
            for J in K + 1 .. N + 1 loop
              B ( ROW( I ) , J) := B ( ROW( I ) , J) - B ( ROW( I ) , K) * B (
                 ROW( K ) , J) ;
            end loop ;
          end if ;
        end loop ;
      end if ;

    --                               FINISHED INNER REDUCTION
    end loop ;

    --                               END OF MAIN REDUCTION LOOP
    --                               BUILD  X  FOR RETURN, UNSCRAMBLING ROWS
    for I in 1 .. N loop
      X ( I ) := B ( ROW( I ) , N + 1) ;
    end loop ;
    return X ;
  end LINEAR_EQUATIONS ;



  function CROUT_SOLUTION ( A : COMPLEX_MATRIX ;
                            Y : COMPLEX_VECTOR ) return COMPLEX_VECTOR is

  --      PURPOSE : SOLVE THE LINEAR SYSTEM OF EQUATIONS WITH COMPLEX
  --                COEFFICIENTS   [A] * |X| = |Y|
  --
  --      INPUT  : THE COMPLEX MATRIX  A
  --               THE COMPLEX VECTOR  Y
  --
  --      OUTPUT : THE COMPLEX VECTOR  X
  --
  --      METHOD : CROUT REDUCTION AND BACK SUBSTITUTION USING
  --               MAXIMUM ELEMENT FOR DIAGONAL
  --
  --      USAGE  :     X := CROUT_SOLUTION ( A , Y ) ;
  --
  --
  --      EXCEPTION : ARRAY_INDEX_ERROR  RASIED IF RESULT CAN NOT BE COMPUTED
  --
  --    WRITTEN BY : JON SQUIRE , 28 MAY 1983
  --                 updated to new complex spec 11/30/88
  --
  --
  --                          LOCAL SPECIFICATIONS

    N : INTEGER := A'LENGTH ;  -- NUMBER OF EQUATIONS
    X : COMPLEX_VECTOR ( 1 .. N ) ;  -- RESULT BEING COMPUTED
    B : COMPLEX_MATRIX ( 1 .. N , 1 .. N + 1 ) ;  -- WORKING MATRIX
    ROW : array ( 1 .. N ) of INTEGER ;  -- ROW INTERCHANGE INDICIES
    HOLD , I_PIVOT : INTEGER ;  -- PIVOT INDICIES
    PIVOT : COMPLEX ;  -- PIVOT ELEMENT VALUE
    ABS_PIVOT : REAL ;
    SUM:COMPLEX;
    I:INTEGER;
    J:INTEGER;
  begin
    if A'LENGTH ( 1 ) /= A'LENGTH ( 2 ) then
      raise ARRAY_INDEX_ERROR ;
    end if ;

    --                        BUILD WORKING DATA STRUCTURE
    for I in 1 .. N loop
      for J in 1 .. N loop
        B ( I , J ) := A ( I - 1 + A'FIRST( 1 ) , J - 1 + A'FIRST ( 2 )) ;
      end loop ;
      B ( I , N + 1 ) := Y ( I - 1 + Y'FIRST ) ;
    end loop ;

    --                        SET UP ROW  INTERCHANGE VECTORS
    for K in 1 .. N loop
      ROW ( K ) := K ;
    end loop ;

    --                         BEGIN MAIN REDUCTION LOOP
    for M in 1..N loop
      --                             FIND LARGEST ELEMENT FOR PIVOT
      PIVOT := B ( ROW( M ) , M) ;
      ABS_PIVOT := abs ( PIVOT ) ;
      I_PIVOT := M ;
      for I in M .. N loop
        if abs ( B( ROW( I ) , M)) > ABS_PIVOT then
          I_PIVOT := I ;
          PIVOT := B ( ROW( I ) , M) ;
          ABS_PIVOT := abs ( PIVOT ) ;
        end if ;
      end loop ;

      --                                HAVE PIVOT, INTERCHANGE ROW POINTERS
      HOLD := ROW ( M ) ;
      ROW ( M ) := ROW ( I_PIVOT ) ;
      ROW ( I_PIVOT ) := HOLD ;

      --                                CHECK FOR NEAR SINGULAR
      if ABS_PIVOT < 1.0E-5 then -- degeneration
        for J in M + 1 .. N + 1 loop
          B ( ROW( M ) , J) := COMPOSE_FROM_CARTESIAN( 0.0 , 0.0 ) ;
        end loop ;
      else -- no degeneration
        -- below diagonal computations
        J := M;
        for I in M..N loop
          SUM := COMPOSE_FROM_CARTESIAN(0.0,0.0);
          for K in 1..J-1 loop
             SUM := SUM + B(ROW(I),K)*B(ROW(K),J);
          end loop; -- k
          B(ROW(I),J) := B(ROW(I),J) - SUM;
        end loop; -- i

        -- above diagonal computations
        I := M;
        for J in M+1..N+1 loop
          SUM := COMPOSE_FROM_CARTESIAN(0.0,0.0);
          for K in 1..I-1 loop
            SUM := SUM + B(ROW(I),K)*B(ROW(K),J);
          end loop; -- k
          B(ROW(I),J) := (B(ROW(I),J)-SUM)/B(ROW(I),I);
        end loop; -- j
      end if; -- degeneration
    end loop; -- m

    -- back substitute
    for I in reverse 1..N loop
      SUM := COMPOSE_FROM_CARTESIAN(0.0,0.0);
      for K in I+1..N loop
        SUM := SUM + B(ROW(I),K)*X(K);
      end loop; -- k
      X(I) := B(ROW(I),N+1) - SUM;
    end loop; -- i
    return X;
  end CROUT_SOLUTION;


  function DETERMINANT ( A : COMPLEX_MATRIX ) return COMPLEX is

    --                          LOCAL SPECIFICATIONS
    N : INTEGER := A'LENGTH ;  -- NUMBER OF EQUATIONS
    D : COMPLEX := COMPOSE_FROM_CARTESIAN( 1.0 , 0.0 ) ;  -- DETERMINANT
    SINGULAR : COMPLEX := COMPOSE_FROM_CARTESIAN( 0.0 , 0.0 ) ;
    B : COMPLEX_MATRIX ( 1 .. N , 1 .. N ) ;  -- WORKING MATRIX
    ROW : array ( 1 .. N ) of INTEGER ;  -- ROW INTERCHANGE INDICIES
    HOLD , I_PIVOT : INTEGER ;  -- PIVOT INDICIES
    PIVOT : COMPLEX ;  -- PIVOT ELEMENT VALUE
    ABS_PIVOT : REAL ;
  begin
    if A'LENGTH ( 1 ) /= A'LENGTH ( 2 ) then
      raise ARRAY_INDEX_ERROR ;
    end if ;
    B := A ;
    D := COMPOSE_FROM_CARTESIAN( 1.0 , 0.0 ) ;

    --                        SET UP ROW  INTERCHANGE VECTORS
    for K in 1 .. N loop
      ROW ( K ) := K ;
    end loop ;

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

      --                             FIND LARGEST ELEMENT FOR PIVOT
      PIVOT := B ( ROW( K ) , K) ;
      ABS_PIVOT := abs ( PIVOT ) ;
      I_PIVOT := K ;
      for I in K .. N loop
        if abs ( B( ROW( I ) , K)) > ABS_PIVOT then
          I_PIVOT := I ;
          PIVOT := B ( ROW( I ) , K) ;
          ABS_PIVOT := abs ( PIVOT ) ;
        end if ;
      end loop ;

      --                                HAVE PIVOT, INTERCHANGE ROW POINTERS
      if I_PIVOT /= K then
        HOLD := ROW ( K ) ;
        ROW ( K ) := ROW ( I_PIVOT ) ;
        ROW ( I_PIVOT ) := HOLD ;
        D := - D ;
      end if ;

      --                                CHECK FOR NEAR SINGULAR
      if ABS_PIVOT < 1.0E-7 then
        return SINGULAR ;
      else
        D := D * PIVOT ;

        --                                REDUCE ABOUT PIVOT
        for J in K + 1 .. N loop
          B ( ROW( K ) , J) := B ( ROW( K ) , J) / B ( ROW( K ) , K) ;
        end loop ;

        --                               INNER REDUCTION LOOP
        for I in 1 .. N loop
          if I /= K then
            for J in K + 1 .. N loop
              B ( ROW( I ) , J) := B ( ROW( I ) , J) - B ( ROW( I ) , K) * B (
                 ROW( K ) , J) ;
            end loop ;
          end if ;
        end loop ;
      end if ;

    --                               FINISHED INNER REDUCTION
    end loop ;

    --                               END OF MAIN REDUCTION LOOP
    return D * B ( ROW( N ) , N) ;
  end DETERMINANT ;


  function INVERSE ( A : COMPLEX_MATRIX ) return COMPLEX_MATRIX is

  --      PURPOSE : INVERT AN N BY N MATRIX
  --
  --      INPUT : THE MATRIX  A
  --
  --      OUTPUT : THE INVERSE OF MATRIX  A
  --
  --      METHOD : GUASS-JORDAN ELIMINATION USING MAXIMUM ELEMENT
  --               FOR PIVOT.
  --
  --    WRITTEN BY : JON SQUIRE , 3 FEB 1983
  --                              30 Nov 1988, revised for new COMPLEX spec
  --
  --
  --                                              LOCAL SPECIFICATIONS

    N : INTEGER := A'LENGTH ;  -- SIZE OF MATRIX
    AA : COMPLEX_MATRIX ( 1 .. N , 1 .. N ) ;  -- WORKING MATRIX
    -- ROW,COLUMN INTERCHANGE INDICIES
    ROW , COL : array ( 1 .. N ) of INTEGER ;
    TEMP : COMPLEX_VECTOR ( 1 .. N ) ;  -- TEMP ARRAY FOR UNSCRAMBLING
    HOLD , I_PIVOT , J_PIVOT : INTEGER ;  -- PIVOT INDICIES
    PIVOT : COMPLEX ;  -- PIVOT ELEMENT VALUE
    ABS_PIVOT : REAL ;
  begin
    if A'LENGTH ( 1 ) /= A'LENGTH ( 2 ) then
      raise ARRAY_INDEX_ERROR ;
    end if ;
    AA := A ;

    --                  SET UP ROW AND COLUMN INTERCHANGE VECTORS
    for K in 1 .. N loop
      ROW ( K ) := K ;
      COL ( K ) := K ;
    end loop ;

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

      --                             FIND LARGEST ELEMENT FOR PIVOT
      PIVOT := AA ( ROW( K ) , COL ( K )) ;
      I_PIVOT := K ;
      J_PIVOT := K ;
      for I in K .. N loop
        for J in K .. N loop
          ABS_PIVOT := abs ( PIVOT ) ;
          if abs ( AA( ROW( I ) , COL ( J ))) > ABS_PIVOT then
            I_PIVOT := I ;
            J_PIVOT := J ;
            PIVOT := AA ( ROW( I ) , COL ( J )) ;
          end if ;
        end loop ;
      end loop ;
      HOLD := ROW ( K ) ;
      ROW ( K ) := ROW ( I_PIVOT ) ;
      ROW ( I_PIVOT ) := HOLD ;
      HOLD := COL ( K ) ;
      COL ( K ) := COL ( J_PIVOT ) ;
      COL ( J_PIVOT ) := HOLD ;

      --                                REDUCE ABOUT PIVOT
      AA ( ROW( K ) , COL ( K )) := COMPOSE_FROM_CARTESIAN(1.0,0.0) / PIVOT ;
      for J in 1 .. N loop
        if J /= K then
          AA ( ROW( K ) , COL ( J )) := AA ( ROW( K ) , COL ( J )) * AA ( ROW
            ( K ) , COL ( K )) ;
        end if ;
      end loop ;

      --                       INNER REDUCTION LOOP
      for I in 1 .. N loop
        if K /= I then
          for J in 1 .. N loop
            if K /= J then
              AA ( ROW( I ) , COL ( J )) := AA ( ROW( I ) , COL ( J )) - AA (
                 ROW( I ) , COL ( K )) * AA ( ROW( K ) , COL ( J )) ;
            end if ;
          end loop ;
          AA ( ROW( I ) , COL ( K )) := - AA ( ROW( I ) , COL ( K )) * AA (
             ROW( K ) , COL ( K )) ;
        end if ;
      end loop ;

    --                               FINISHED INNER REDUCTION
    end loop ;

    --                   END OF MAIN REDUCTION LOOP
    --                              UNSCRAMBLE ROWS
    for J in 1 .. N loop
      for I in 1 .. N loop
        TEMP ( COL( I )) := AA ( ROW( I ) , J) ;
      end loop ;
      for I in 1 .. N loop
        AA ( I , J ) := TEMP ( I ) ;
      end loop ;
    end loop ;

    --                              UNSCRAMBLE COLUMNS
    for I in 1 .. N loop
      for J in 1 .. N loop
        TEMP ( ROW( J )) := AA ( I , COL( J )) ;
      end loop ;
      for J in 1 .. N loop
        AA ( I , J ) := TEMP ( J ) ;
      end loop ;
    end loop ;
    return AA ;
  end INVERSE ;

end GENERIC_COMPLEX_LINEAR_EQUATIONS ;