-- test_svd_package_r.ada
-- Singular Value Decomposition
-- Given matrix A m by n, find matricies U, W, V such that
-- A = U * W * V**T  (possibly permuted) where:
-- W is an n by n diagonal matrix of non negative elements, "vector"
--   these are the singular values of A in decending order
-- U is an m by n column orthnormal matrix, U**T * U = I
-- V is an n by n orthonormal matrix, V * V**T = I and V**T * V= I
-- A**-1 = V * [diag(1/W(k)] * U**T
-- A * x = y solved by
--     x = V * [diag(1/W(k))] * (U**T * y)  via  SVD_SOLVE using U

with text_io; use text_io;
with real_arrays; use real_arrays;
with arrays_io;
with svd_package; use svd_package;
with random_numbers; use random_numbers;
with real_linear_equations; use real_linear_equations;
procedure test_svd_package_r is
  subtype real is float;
  package arr_io is new arrays_io(real, real_vector, real_matrix);
  use arr_io;  
  package int_io is new integer_io(integer);
  use int_io;

  n : integer;
  m : integer;
  debug : boolean := false; -- compile time switch
  G : Generator;
  Det : float;
  finished : exception;

  function norm1(A : Real_Matrix) return real is
    value : real := 0.0;
  begin
    for i in A'range(1) loop
      for j in A'range(2) loop
        if abs A(i,j) > value then
          value := abs A(i,j);
        end if;
      end loop;
    end loop;
    return value;
  end norm1;

  function norm1(A : Real_vector) return real is
    value : real := 0.0;
  begin
    for i in A'range loop
        if abs A(i) > value then
          value := abs A(i);
        end if;
    end loop;
    return value;
  end norm1;

begin
  put_line("test_svd_package_r, random matricies");
  n := 10;
  Reset(G);
  while n < 300 loop
    -- generate random data sets
    m := n;
    Put_Line("Test rank " & integer'image(n) & "  random matrix");
    declare
      A : real_matrix(1..m,1..n);
      AA : real_matrix(1..m,1..n);
      AI : real_matrix(1..n,1..m);
      IDENT : real_matrix(1..n,1..n);
      U : real_matrix(1..m,1..n);
      UT : real_matrix(1..m,1..n);
      V : real_matrix(1..n,1..n);
      VT : real_matrix(1..n,1..n);
      W : real_vector(1..n);
      X : real_vector(1..n);
      Y : real_vector(1..n);
      Z : real_vector(1..n);
      sum : real;
      INFO : integer;
    begin
      for i in A'range(1) loop
        for j in A'range(2) loop
          A(i,j) := float(Random(G));
        end loop;
      end loop;
      for i in Y'range loop
        Y(i) := real(i);
      end loop;
      if debug then
        put_line("A=  ");
        put(A);
      end if;

      Put(Norm1(A)); Put(", ");
      Put(Determinant(A)); Put_Line(" = Baseline Norm, Determinant");
      Put(Norm1(Y-A*Linear_Equations(A,Y)));
      put_line("= norm1 of Baseline Y-A*Linear_Equations(A,Y)");

      svd(A, U, V, W, INFO);
      VT := transpose(V);
      UT := transpose(U);
      if debug then
        put_line("Ada SVD INFO= " & integer'image(INFO));
        put_line("U=   ");
        put(U);
        put_line("V=   ");
        put(V);
        put_line("W=   ");
        put(W);
        new_line;
      end if;
      for i in 1..m loop
        for j in 1..n loop
          sum := 0.0;
          for k in 1..n loop
            sum := sum+U(i,k)*W(k)*VT(k,j);
          end loop;
          AA(i,j):= sum;
        end loop;
      end loop;
      put(norm1(AA-A));
      put_line("= norm1 of origonal matrix minus U * diag(W) * transpose(V)");
      if debug then
        put_line("AA=U*W*VT using loops");
        put(AA);
        new_line;
        -- AA := U * W * VT; -- does not work
        -- put_line("AA=U*W*VT using matrix equations");
        -- put(AA);
        -- new_line;
      end if;
      X := svd_solve(U, V, W, Y);
      Z := A * X;
      if debug then
        put_line("Y=  ");
        put(Y);
        put_line("X=  ");
        put(X);
        put_line("Z= ");
        put(Z);
        put_line("Z should equal Y, based on A");
        new_line;
      end if;
      put(norm1(Z-Y));
      put_line("= norm1 of origonal matrix * X - Y, solving equations");

      AI := svd_inverse(U, V, W);
      if debug then
        put_line("svd inverse AI");
        put(AI);
        new_line;
        put_line("A * AI  should be identity");
        IDENT := A * AI;
        put(IDENT);
        new_line;
        put_line("AI * A ");
        IDENT := AI * A;
        put(IDENT);
        new_line;
      end if;
      put(norm1(A*AI-Identity_Matrix(n)));
      put_line("= norm1 of origonal matrix * inverse - identity");
      put(norm1(AI*A-Identity_Matrix(n)));
      put_line("= norm1 of inverse * origonal matrix - identity");
      put(norm1(AA*AI-Identity_Matrix(n)));
      put_line("= norm1 of AA * AI - I");

      if debug then
        put_line("U * UT = I ");
        IDENT := U * transpose(U);
        put(IDENT);
        new_line;
        put_line("VT * V = I ");
        IDENT := VT * V;
        put(IDENT);
        new_line;
        put_line("V * VT = I ");
        IDENT := V * VT;
        put(IDENT);
        new_line;
      end if;
      put(norm1(U*UT-Identity_matrix(n)));
      put_line("= norm1 of U * UT - I");
      put(norm1(V*VT-Identity_matrix(n)));
      put_line("= norm1 of V * VT - I");
      put(norm1(VT*V-Identity_matrix(n)));
      put_line("= norm1 of VT * V - I");
    exception
      when Real_Arrays.Array_Index_error =>
        put_line("Array_Index_Error from package");
      when Matrix_Solution_error =>
        put_line("singular matrix");
        new_line;
    end; -- tests on each data set
    n := n + 10;
  end loop; -- all data sets finished
  put_line("test_svd_package_r completed");
exception
  when finished =>
    put_line("test_svd_package_r ends");
--  when others =>
--    put_line("FAILED  premature exit of test by unhandled exception");
end test_svd_package_r;