-- test_sparse.adb

with Real_Arrays;
use Real_Arrays;
with Sparse;
with Ada.Text_Io;
use Ada.Text_Io;

procedure  Test_Sparse is
  Nrow : Integer := 6;
  package A is new Sparse(Nrow);
  Val : Real;
  X : Real_Vector(0..Nrow-1);
  Y : Real_Vector(0..Nrow-1);
  YC : Real_Vector(0..Nrow-1);
  AA : Real_Matrix(0..Nrow-1,0..Nrow);
begin
    Put_Line("test_sparse_ada.out nrow="&Integer'Image(Nrow));
    Put_Line("test put(0,4,3.14) and get in A");
    A.put(0,4,3.14);
    val := A.get(0,4);
    Put_Line("test put in row, 3.14 should be "&Real'Image(Val));
    Put_Line("test put(5,0,6.28) and get in A");
    A.put(5,0,6.28);
    val := A.get(5,0);
    Put_Line("test put in row, 6.28 should be"&Real'Image(Val));
    A.write_all;
    A.clear;

    for I in 0..Nrow-1 loop
      for J in I-1..I+1 loop
        val := 1.0+real(i*i+(j+1)*(j+1));
        if I>=0 and J>=0 and i<Nrow and j<Nrow then
          A.add(I, J, Val);
        end if;
        A.add(i, nrow, 8.0+4.0*Real(I));
      end loop;
    end loop;
    Put_Line("initial set of equations in A");
    A.write_all;
    A.Export(AA);
    A.GetRHS(Y);      -- for testing
    Put_Line("calling sparse.simeq");
    A.Simeq(X); -- modifies, destroys, A
    for I in 0..Nrow-1 loop
      Put_Line("X("&Integer'Image(I)&")="&Real'Image(X(I)));
    end loop;
    Put_Line("A after simeq");
    A.Write_All;

    Put_Line(" ");
    Put_Line("test multiply and simeq");
    A.Clear;
    Put_Line("cleared");
    A.Import(AA); -- reset equations
    Put_Line("imported");
    A.Write_All;
    A.Multiply(X,YC); -- check solution
    for I in 0..Nrow-1 loop
      if abs(Y(I)-Yc(I))> 0.00001 then
        Put_Line("either multiply or simeq failed");
        Put_Line("i="&Integer'Image(I)&", Y[i]="&Real'Image(Y(I))&
             " yet found "&Real'Image(Yc(I)));
      end if;
    end loop;
  Put_Line("test_sparse finished");
end Test_Sparse;