--  bsimeq_2.adb  solve simultaneous equations shared memory, parallel tasks
-- solve real linear equations for X where Y = A * X
-- method: Gauss-Jordan elimination using maximum pivot
-- usage:  X = psimeq(A,Y);
--    Translated to java by : Jon Squire , 26 March 2003
--    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 made parallel October 2008
--    method : GAUSS-jordan elimination using maximum element for pivot.

with Ada.Exceptions; use Ada;
--with Ada.Numerics.Real_Arrays;
--use  Ada.Numerics.Real_Arrays;
with Real_Arrays;
use  Real_Arrays;
with Ada.Text_IO;
use  Ada.Text_IO;
with brAda.Synchronous_Barriers; use brAda.Synchronous_Barriers;

function bsimeq_2 (NP : Positive;  A : real_matrix ; Y : real_vector )
                 return real_vector is
  -- NP  number of processors, modify to suit
  N : integer := A'length(1);     -- number of equations
  Part : Integer := N/NP;         -- number of rows per task, or +1
  Remain : Integer;               -- balance rows
  X : real_vector(1..N);          -- result being computed
  B : real_matrix(1..N,1..N+1) ;  -- working matrix, part in each task
  Prow : array(1..n) of Integer ; -- pivot row if not -1, variable number
  Frow : array(1..n) of Integer ; -- finished row if not 0
  Srow : array(1..NP+1) of Integer; -- starting row for each task
  Smax_Pivot : Real_Vector(1..NP);  -- max pivot each task
  Spivot : array(1..NP) of Integer; -- row for max pivot each task
  K_Row : Integer := 1;             -- pivot row
  I_pivot : integer;                -- pivot indicies
  abs_pivot : real;                -- abs of pivot element
  matrix_data_error : exception;
  Barr1 : Synchronous_Barrier (Number_Waiting => NP);
  Barr2 : Synchronous_Barrier (Number_Waiting => NP);

  task type slave is
    entry Start(Me : Integer);
  end slave;
  task body slave is
    myid  : Integer;
    Released_Last : Boolean;
  begin

    Accept Start(Me : Integer) do
      myid := Me;
    end Start;
    -- setup working array, one time
    -- build working data structure
    for I in Srow(myid)..Srow(myid+1)-1 loop
      for J in 1..N loop
        B(I,J) := A(I,J);
      end loop;
      B(I,N+1) := Y(I);
    end loop;

    -- start column processing find max pivot, wait, reduce. wait
    for K_Col in 1..N loop
      -- find max of tasks max pivot
      -- set smax_pivot(myid)
      -- set spivot(myid)
      Spivot(myid) := Srow(myid);
      Smax_pivot(myid) := 0.0; -- singular caught elsewhere
      for I in Srow(myid)..Srow(myid+1)-1 loop
        if Frow(I)=0 and abs(B(I,K_Col)) > smax_pivot(myid) then
          spivot(myid) := I;
          smax_pivot(myid) := abs(B(I,K_Col));
        end if;
      end loop;

      Wait_For_Release (Barrier       => Barr1,
                        Released_Last => Released_Last);

      if Released_Last then

         -- find global pivot row and divide
         abs_pivot := smax_pivot(1);
         I_pivot := Spivot(1);
         for I in 2..NP loop
           if smax_pivot(I)>Abs_Pivot then
             I_pivot := Spivot(I);
             abs_pivot := smax_pivot(I);
           end if;
         end loop;

         -- have pivot, record for X at end
         K_Row := I_pivot;
         prow(K_Col) := K_Row;
         frow(K_Row) := 1;

         -- check for near singular
         if abs_pivot < 1.0E-12 then
           for J in K_Col+1..N+1 loop
             B(K_Row,j) := 0.0;
           end loop;
           Put_Line("redundant row (singular) "&Integer'Image(K_Row));
           -- singular, set row to zero, including solution
         else
           -- reduce about pivot
           for J in K_Col+1..N+1 loop
             B(K_Row,J) := B(K_Row,J) / B(K_Row,K_Col);
           end loop;
         end if;

      end if;

      Wait_For_Release (Barrier       => Barr2);

      --  inner reduction loop
      for I in Srow(myid)..Srow(myid+1)-1 loop
        if I /= K_Row then
          for J in K_Col+1..N+1 loop
            B(I,J) := B(I,J) - B(I,K_Col) * B(K_Row,J);
          end loop;
        end if;
      end loop;
      --  finished inner reduction
    end loop; -- N times, then terminate

   exception
      when E : others =>
         Put_Line ("Slave died: " & Exceptions.Exception_Information (E));
  end slave;

begin -- psimeq
  pragma Assert (A'length(1)=A'length(2) and A'length(1)= Y'length);

  declare
     G : array(1..NP) of slave; -- Start NP tasks running
  begin

     -- set up row range for each task
     Remain := N-Part*NP;
     Srow(1) := 1;
     for I in 2..NP loop
       Srow(I) := Srow(I-1)+Part;
       if I <= Remain+1 then
         Srow(I) := Srow(I)+1;
       end if;
     end loop;
     Srow(NP+1) := N+1; -- may not be plus part

     for I in 1..NP loop
        Put_Line("thread "&Integer'Image(I)&", from "&Integer'Image(Srow(I))&
                   " to "&Integer'Image(Srow(I+1)-1));
     end loop;

     -- set up pivot row and finished row
     for K in 1..N loop
       Prow(K) := -1;
       Frow(K) := 0;
     end loop;

     for T in 1..NP loop -- let tasks build working rows in B
       G(T).Start(T);
     end loop;
  end;

  -- Note, we cannot get here unless all the tasks have terminated,
  -- (since enclosed by an inner scope)
  -- therefore no need for a barrier here.

  --  build  X  for return, unscrambling rows
  for I in 1..N loop
    X(i) := B(Prow(i),N+1);
  end loop;
  return X;
end bsimeq_2;