-- inverse.adb  extracted from generic_real_linear_equations.adb
with Real_Arrays;
use  Real_Arrays;

function INVERSE(A : REAL_MATRIX)return REAL_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.
  --
  --      EXCEPTION : MATRIX_DATA_ERROR RAISED IF INVERSE CAN NOT BE COMPUTED
  --                  ARRAY_INDEX_ERROR IF 'A' IS NOT SQUARE
  --
  --      SAMPLE USE :  NEW_MATRIX := INVERSE(OLD_MATRIX);
  --                    MATRIX := INVERSE(MATRIX)* ANOTHER_MATRIX;
  --
  --      WRITTEN BY : JON SQUIRE,3 FEB 1983

  N : INTEGER := A'LENGTH;                 -- SIZE OF MATRIX
  AA : REAL_MATRIX(1..N,1..N);    -- WORKING MATRIX
  TEMP : REAL_VECTOR(1..N);           -- TEMPORARY FOR UNSCRAMBLING ROWS
  ROW,COL : array(1..N)of INTEGER; -- ROW,COLUMN INTERCHANGE INDICIES
  HOLD,I_PIVOT,J_PIVOT : INTEGER;      -- PIVOT INDICIES
  PIVOT : REAL;                            -- PIVOT ELEMENT VALUE
  ABS_PIVOT : REAL;                        -- ABS OF PIVOT ELEMENT
  NORM1 : REAL := 0.0;                     -- 1 NORM OF MATRIX
  MATRIX_DATA_ERROR : exception;
  ARRAY_INDEX_ERROR : exception;
begin
  if A'LENGTH(1)/= A'LENGTH(2)then
    raise ARRAY_INDEX_ERROR;
  end if;
  --                              BUILD WORKING DATA STRUCTURE
  AA := A;
  for I in 1..N loop
    for J in 1..N loop
      if abs AA(I,J)> NORM1 then
        NORM1 := abs AA(I,J);
      end if;
    end loop;
  end loop;
  --                              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;
    --                            TEST FOR SINGULAR
    if ABS_PIVOT < REAL'EPSILON * NORM1 then
      raise MATRIX_DATA_ERROR;
    end if;

    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)) := 1.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;
  --                              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;