-- hilbert_inverse.adb   watch accuracy decrease as n increases
-- A(i,j) := 1.0/(real(i)+real(j))  1/2 1/3 1/4 ...
--                                  1/3 1/4 1/5 ...
with Ada.Text_Io;
use Ada.Text_Io;
with Real_Arrays;
use Real_Arrays;
-- with Inverse; tailored version physically included for data collection

procedure Hilbert_Inverse is
  N : Integer;
  Abs_Pivot : Real;               -- ABS OF PIVOT ELEMENT
  Norm1: Real := 0.0;             -- 1 NORM OF MATRIX
  Avgerr : Real := 0.0;           -- abs sum AII-I /n^2
  Maxerr : Real := 0.0;           -- max error in AII-I
  Mindiag : Real := 0.0;          -- smallest Abs_Pivot
  State : Real := 30000000.0;     -- for random number generator

-- inverse.adb  extracted from generic_real_linear_equations.adb

function INVERSE(A : REAL_MATRIX)return REAL_MATRIX is
  --      PURPOSE : INVERT AN N BY N MATRIX
  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
  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;
  Norm1 := 0.0;
  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 NEAR SINGULAR
    if ABS_PIVOT < mindiag then
      Mindiag := Abs_Pivot;
    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;

-- udrnrt.adb uniformly disributed random numbers 0.0 .. 1.0
function Udrnrt return Real is
   Modulus : Real := 2147483647.0; -- 2**31-1
   Multiplier : Real := 16807.0; -- prime to prime power
   T : Long_Float;
   I : Integer;
begin
   T := State * Multiplier;
   I := Integer(T / Modulus);
   State := T - (Long_Float(I) * Modulus);
   if State < 0.0 then
      State := State + Modulus;
   end if;
   return State / Modulus;
end Udrnrt;


begin
   Put_Line("hilbert_inverse.adb accuracy test");
   Put_Line("IEEE 64-bit floating point, Hilbert vs random matrix");
   Put_Line("for N by N matrices, N=2, 4, 8, 16, 32, 64, 126, 256");
   New_Line;
   for NI in 1..8 loop -- will be powers of 2
      N := 1;
      for NJ in 1..NI loop
         N := N*2;
      end loop;
      declare
         A : Real_Matrix(1..N,1..N);
         AI : Real_Matrix(1..N,1..N);
         AII : Real_Matrix(1..N,1..N);
         Err : Real;
      begin
         -- Hilbert (then later, random)
         for I in 1..N loop
            for J in 1..N loop
               A(I,J) := 1.0/(Real(I)+Real(J));
            end loop; -- J
         end loop; -- I
         Avgerr := 0.0;
         Maxerr := 0.0;
         Mindiag := 1.0e20;
         AI := Inverse(A);
         AII := A*AI;
         for I in 1..N loop
            AII(I,I) := AII(I,I)-1.0;
         end loop;
         for I in 1..N loop
            for J in 1..N loop
               Err := abs(AII(I,J));
               Avgerr := Avgerr+Err;
               if Err > Maxerr then
                  Maxerr := Avgerr;
               end if;
            end loop; -- I
         end loop; -- J
         Avgerr := Avgerr/(Real(N)*Real(N));
         Put_Line("N="&Integer'Image(N)&", norm="&Real'Image(Norm1)
                  &", mindiag="&Real'Image(Mindiag));
         Put_Line("     avgerr="&Real'Image(Avgerr)
                  &", maxerr="&Real'Image(Maxerr));

         -- Now, uniformly distributed random number 0.0..1,0
         for I in 1..N loop
            for J in 1..N loop
               A(I,J) := Udrnrt;
            end loop; -- J
         end loop; -- I
         Avgerr := 0.0;
         Maxerr := 0.0;
         Mindiag := 1.0e20;
         AI := Inverse(A);
         AII := A*AI;
         for I in 1..N loop
            AII(I,I) := AII(I,I)-1.0;
         end loop;
         for I in 1..N loop
            for J in 1..N loop
               Err := abs(AII(I,J));
               Avgerr := Avgerr+Err;
               if Err > Maxerr then
                  Maxerr := Avgerr;
               end if;
            end loop; -- I
         end loop; -- J
         Avgerr := Avgerr/(Real(N)*Real(N));
         Put_Line("random  norm="&Real'Image(Norm1)
                  &", mindiag="&Real'Image(Mindiag));
         Put_Line("     avgerr="&Real'Image(Avgerr)
                  &", maxerr="&Real'Image(Maxerr));
         New_Line;
      end; -- declare
   end loop; -- NI
   Put_Line("hilbert_inverse ends");
end Hilbert_Inverse;