-- cxinverse.adb
with Ada.Numerics.Long_Complex_Types;
use  Ada.Numerics.Long_Complex_Types;
with Complex_Arrays; -- mine
use  Complex_Arrays; -- mine

function CXINVERSE ( 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 : GAUSS-JORDAN ELIMINATION USING MAXIMUM ELEMENT
  --               FOR PIVOT.

  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 : Long_float;
begin
  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 CXINVERSE ;