-- sparse.adb     Sparse storage of a matrix
--                specifically designed for solving PDE
--                simultaneous equations.
--                zero based subscripting
--                column(nrow) is right hand side

with Ada.Text_Io;
use Ada.Text_Io;

package body Sparse is

  -- new cell'(j, val, Next)  cell stuff, abstraction
  -- Cell operations, uses private part of package spec

  procedure SetVal(C : in out Link; Aj : Integer; Aval : real) is
                   -- creates Aj if not in list
    P : Link;
  begin
    if C = null then
      C := new Cell'(Aj, Aval, null);
      return;
    end if;
    P := C; -- walk the list
    while True loop
      if Aj = P.J then
          P.Val := Aval;
          return;
      end if;
      if Aj < P.J then -- we passed it, insert before
        P.Next := new Cell'(P.J, P.Val, P.Next);
        P.J := Aj;
        P.Val := Aval;
        return ;
      end if;
      if P.Next = null then -- insert at end
        P.Next := new Cell'(Aj, Aval, null);
        return;
      end if;
      P := P.Next; -- move along list
    end loop; -- while
  end Setval; -- Aval to Sj in Cell C

  procedure AddVal(C : in out Link; Aj : Integer; Aval : real) is
                   -- creates Aj if not in list
    P : Link;
  begin
    if C = null then
      C := new Cell'(Aj, Aval, null);
      return;
    end if;
    P := C; -- walk the list
    while True loop
      if Aj = P.J then
          P.Val := P.Val + Aval;
          return;
      end if;
      if Aj < P.J then -- we passed it, insert before
        P.Next := new Cell'(P.J, P.Val, P.Next);
        P.J := Aj;
        P.Val := Aval;
        return ;
      end if;
      if P.Next = null then -- insert at end
        P.Next := new Cell'(Aj, Aval, null);
        return;
      end if;
      P := P.Next; -- move along list
    end loop; -- while
  end Addval; -- Aval added to Sj in Cell C

  function GetVal(C : Link; Aj : Integer) return Real is
    P : Link;
  begin
    P := C; -- walk the list
    while P /= null loop
      if Aj = P.J then
        return P.Val;
      elsif Aj < P.J then -- we passed it
        return 0.0;
      end if;
      P := P.Next; -- move along list
    end loop; -- while
    return 0.0;
  end Getval; -- Aval to Sj in Cell C



  -- exported via package spec

  procedure clear is
  begin
    for I in 0..Nrow-1 loop
      A(I) := null;
    end loop;
  end Clear;

  procedure Write_All is
    P : Link;
  begin
    Put_Line("Sparse Write_All");
    for I in 0..Nrow-1 loop
      P := A(I);
      while P /= null loop
        Put_Line("i="&Integer'Image(I)&", j="&Integer'Image(P.J)&
                   ", val="&Real'Image(P.Val));
        P := P.Next;
      end loop;
    end loop;
  end Write_All;

  procedure Import(M : in Real_Matrix) is
  begin
    if M'Length(1) < Nrow or M'Length(2) < Nrow+1 then
       Put_Line("Sparse Import error need at least"&Integer'Image(Nrow)&" rows by "&
                Integer'Image(Nrow+1)&" columns");
    end if;
    for I in 0..Nrow-1 loop
      for J in 0..Nrow loop
        if M(I+M'First(1),J+M'First(2)) /= 0.0 then
          setval(A(I),J,M(I+M'First(1),J+M'First(2)));
        end if;
      end loop;
    end loop;
  end Import;

  procedure Export(M : out Real_Matrix) is
    P : Link;
  begin
    for I in M'Range(1) loop
      for J in M'Range(2) loop
        M(I,J) := 0.0;
      end loop;
    end loop;
    if M'Length(1) < Nrow or M'Length(2) < Nrow+1 then
      Put_Line("Sparse Export error, matrix too small, must be at least "&
               Integer'Image(Nrow)&" by "&Integer'Image(Nrow+1));
    end if;
    for I in 0..Nrow-1 loop
      P := A(I);
      while P /= null loop
        M(I+M'First(1),P.J+M'First(2)) := P.Val;
        P := P.Next;
      end loop;
    end loop;
  end Export;

  function  Get(I:Integer; J:Integer) return Real is -- zero for no such J
  begin
    if I < 0 or I > Nrow then
      Put_Line("Sparse Get error i="&Integer'Image(I)&", j="&Integer'Image(J));
      return 0.0;
    end if;
    return GetVal(A(I),J);
  end Get;

  procedure Put(I:Integer; J:Integer; Val:real) is
  begin
     if I < 0 or I > Nrow then
       Put_Line("Sparse Put error i="&Integer'Image(I)&", j="&Integer'Image(J));
       return;
     end if;
     SetVal(A(I), J, Val);
  end Put;

  procedure Add(I:Integer; J:Integer; Val:real) is
  begin
    if I < 0 or I > Nrow then
      Put_Line("Sparse Add error i="&Integer'Image(I)&", j="&Integer'Image(J));
      return;
    end if;
    AddVal(A(I), J, Val);
  end Add;

  procedure  GetRHS(Y : out Real_Vector) is
  begin
    if Y'Length < Nrow then
       Put_Line("Sparse GetRHS error need at least"&Integer'Image(Nrow)&" rows");
    end if;
    for I in 0..Nrow-1 loop
      Y(I+Y'First) := GetVal(A(I),Nrow);
    end loop;
  end GetRHS;

  procedure SetRHS(Y:Real_Vector) is
  begin
    if Y'Length < Nrow then
      Put_Line("Sparse SetRHS error need at least"&Integer'Image(Nrow)&
               " elements");
      return;
    end if;
    for I in 0..Nrow-1 loop
      setval(A(I),Nrow,Y(I+Y'First));
    end loop;
  end SetRHS;

  procedure Multiply(X : Real_Vector; Y : out Real_Vector) is
    P : Link;
  begin
    if X'Length < Nrow or Y'Length < Nrow then
       Put_Line("Sparse Multiply error need at least"&Integer'Image(Nrow)&
                " elements");
       return;
    end if;
    for I in 0..Nrow-1 loop
      Y(I+Y'First) := 0.0;
    end loop;
    for I in 0..Nrow-1 loop
      P := A(I);
      while P /= null loop
        if P.J >= Nrow then
           exit;
        end if;
        Y(I) := Y(I) + P.Val*X(P.J+X'First);
        P := P.Next;
      end loop;
    end loop;
  end Multiply;

  -- simeq      solve simultaneous equations AX=Y
  -- solve real linear equations for X where Y = A * X
  -- method: Gauss-Jordan elimination using maximum pivot
  --    Translated to Sparse by : Jon Squire , 25 March 2009
  --    First written by Jon Squire December 1959 for IBM 650, translated to
  --    other languages  e.g. Fortran converted to Ada converted to C
  --    then converted to java, then to sparse
  procedure Simeq(X : out Real_Vector) is -- X returned solution
     Hold, I_Pivot : Integer;
     Pivot, Abs_Pivot : Real;
     N : Integer := Nrow;
     Row : array(0..N-1) of Integer; -- row interchange indices
  begin
    -- set up row interchange vectors
    for K in 0..N-1 loop
      Row(k) := k;
    end loop;
    --  begin main reduction loop
    for K in 0..N-1 loop
      -- find largest element for pivot
      Pivot := GetVal(A(Row(K)),K); -- B(row(k))(k);
      Abs_Pivot := abs(Pivot);
      I_Pivot := K;
      for I in K+1..N-1 loop
        -- if A(Row(I),K) = null then exit; end if;
        Pivot := GetVal(A(Row(I)),K); -- B(row(i))(k))
        if abs(Pivot) > Abs_Pivot then -- B(row(i))(k))
          I_Pivot := I;
          Abs_Pivot := abs(Pivot);
        end if;
      end loop;
      -- have pivot, interchange row indices
      Hold := Row(K);
      Row(K) := Row(I_Pivot);
      Row(I_Pivot) := Hold;
      -- check for near singular
      if Abs_Pivot < 1.0E-10 then
        Put_Line("singular at k="&Integer'Image(K)&" redundant row "&
                 Integer'Image(Row(K)));
      else
        -- can divide by pivot
        -- reduce about pivot
        -- for(int j=k+1; j<n+1; j++)
        -- {
        -- B(row(k))(j) = B(row(k))(j) / B(row(k))(k);
          Divide(Row(K),K); -- j loop in here  should use pivot
        -- }
        --  inner reduction loop
        for I in 0..N-1 loop
          if I /= K then
              -- for(int j=k+1; j<n+1; j++)
              -- {
              -- B(row(i))(j) = B(row(i))(j) - B(row(i))(k) * B(row(k))(j);
              Reduce(Row(I),K,Row(K)); -- j loop in here on both
                                       -- row(i) and row(k)
              -- }
          end if;
        end loop;
      end if;
      --  finished inner reduction
    end loop;
    --  end main reduction loop
    --  build  X  for return, unscrambling rows
    for I in 0..N-1 loop
      X(I+X'first) := GetVal(A(Row(I)),N); -- this could be direct if matrix shrunk
    end loop;
  end Simeq;

  procedure Divide(Krow:Integer; K:Integer) is -- works on a Krow j=k+1..n
    P : Link;
    Piv : Real;
  begin
    P := A(Krow);
    while P /= null loop -- find pivot, position to pivot
      if P.J = K then -- found, thus must exists
        Piv := P.Val;
        P := P.Next;
        exit;
      end if;
      P := P.Next;
    end loop;
    while P /= null loop
      P.Val := P.Val/Piv;
      P := P.Next;
    end loop;
  end Divide;

  procedure Reduce(Irow:Integer; K:Integer; Krow:Integer) is
    -- works on complete Irow  j=k+1..n
    -- using Krow j=k+1..n  where it exists
    -- B[row[i]][j] = B[row[i]][j] - B[row[i]][k] * B[row[k]][j];
    --                               Valik
    Pij, Pkj : link;
    Valik : Real := 0.0;
  begin
    Pkj := A(Krow);
    Pij := A(Irow);
    -- find Pik and Valik  constant for  while loop on j
    while Pij /= null loop
      if Pij.J = K then
        Valik := Pij.Val; -- have multiplier
        Pij := Pij.Next;
        exit;
      elsif Pij.J > K then
        return; -- no multiplier
      else
        Pij := Pij.Next;
      end if;
    end loop;
    if Valik = 0.0 then
      return; -- no multiplier
    end if;
    while Pkj /= null loop -- position Pkj
      if Pkj.J > K then
        exit;
      end if;
      Pkj := Pkj.Next;
    end loop;
    while Pij /= null and Pkj /= null loop
    -- find Pij.J = Pkj,J (both move in j)
      if Pij.J < Pkj.J then
        Pij := Pij.Next;
      elsif Pij.J = Pkj.J then
        Pij.Val := Pij.Val - Valik * Pkj.Val;
        Pij := Pij.Next;
        Pkj := Pkj.Next;
      elsif Pij.J > Pkj.J then -- must insert negative value, Pij.val = 0
        Pij.Next := new Cell'(Pij.J, Pij.Val, Pij.Next);
        Pij.J := Pkj.J;
        Pij.Val := -Valik*Pkj.Val;
        Pij := Pij.Next;
        Pkj := Pkj.Next;
      else
        Pkj := Pkj.Next;
      end if;
    end loop;
  end Reduce;
end Sparse;