-- test_svd_package.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 real_arrays_io;
with svd_package; use svd_package;

procedure test_svd_package is
  subtype real is float;
  package arr_io is new real_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;
  Input : text_io.File_TYpe;
  Input_File_Name : string(1..260); --  "-" means read from Standard Input
  Input_file_name_length : Natural;
  Title : string(1..260);
  Title_length : natural;
  debug : boolean := false; -- turned on by $ on front of Input_file_name

  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, enter data file name:");
  get_line(input_file_name, input_file_name_length);
  if input_file_name_length = 0 then
    put_line("no file name provided, ending.");
    raise finished;
  end if;
  put_line(input_file_name(1..input_file_name_length));
  if input_file_name(1) = '$' then
    debug := true;
    input_file_name(1..input_file_name_length-1) :=
                            input_file_name(2..input_file_name_length);
    input_file_name_length := input_file_name_length-1;
  end if;
  -- test for '-'
  open(input, in_file, input_file_name(1..input_file_name_length));
  loop
    -- read data sets from input file, stop on end of file or "0"
    -- get n and m
    get(Input,n);
    if n = 0 then raise finished; end if;

    get(Input,m);
    get_line(Input, Title, Title_length);
    put_line(" n=" & integer'image(n) & "   m=" & integer'image(m) & "   " &
             Title(1..Title_length));
    declare
      A : real_matrix(1..m,1..n);
      AA : real_matrix(1..m,1..n);
      AI : real_matrix(1..n,1..m);
      AAA : real_matrix(1..m,1..n);
      U : real_matrix(1..m,1..n);
      UT : real_matrix(1..m,1..n);
      IDENT : real_matrix(1..n,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
      get(Input,A); -- read in matrix
      for i in Y'range loop
        Y(i) := real(i);
      end loop;
      if debug then
        put_line("A=  ");
        put(A);
      end if;
      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 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
  end loop; -- all data sets finished
  put_line("test_svd_package completed");
exception
  when Text_io.End_Error =>
    put_line("test_svd_package completed");
  when  DATA_ERROR | LAYOUT_ERROR =>
    put_line("bad input data, terminating at:");
    put_line(Title(1..Title_length));
  when finished =>
    put_line("test_svd_package ends");
  when STATUS_ERROR | MODE_ERROR | NAME_ERROR | USE_ERROR | DEVICE_ERROR =>
    put_line("Could not open: " & input_file_name(1..input_file_name_length));
--  when others =>
--    put_line("FAILED  premature exit of test by unhandled exception");
end test_svd_package;