-- psimeq_dbg.java  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_dbg(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 real_arrays; use real_arrays;
with Ada.Text_IO;
use  Ada.Text_IO;

function psimeq_dbg ( A : real_matrix ; Y : real_vector ) return real_vector is
  NP : Integer :=4; -- number of processors, modify to suit
  N : integer := A'length(1);     -- number of equations
  Part : Integer := N/NP;         -- number of rows per task
  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;


  task type slave is
    entry Barr1(Me : Integer);
    entry Barr2(KC : Integer);
    entry Barr3(KR : Integer);
    entry Barr4(K3 : Integer);
  end slave;
  task body slave is
    myid  : Integer;
    K_Row : Integer;
  begin
    Put_Line("slave running");
    Accept Barr1(Me : Integer) do
      myid := Me;
    end Barr1;
    -- setup working array, one time
    Put_Line("slave " & Integer'Image(myid) & " setup ");

    -- 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);
        Put_Line(Integer'Image(Myid)&" B("&Integer'Image(I)&
                 ","&Integer'Image(J)&")="&Real'Image(B(I,J)));
      end loop;
      B(I,N+1) := Y(I);
      Put_Line(Integer'Image(Myid)&" B("&Integer'Image(I)&","&
               Integer'Image(N+1)&")="&Real'Image(B(I,N+1)));
    end loop;
    Put_Line(Integer'Image(Myid)&" built working data structure "&
             Integer'Image(Srow(Myid))&".."&Integer'Image(Srow(myid+1)-1));


    -- start column processing find max pivot, wait, reduce. wait
    for K_Col in 1..N loop
      Accept Barr2(KC : Integer) do
        if KC /= K_COL then
          Put_Line("sync ERROR KC="&Integer'Image(KC)&" /= K_Col="&
                   Integer'Image(K_Col));
        end if;
      end Barr2;
      Put_Line("slave "&Integer'Image(myid)&" finding pivot "&
                           Integer'Image(K_Col));

      -- 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;
      Put_Line(Integer'Image(Myid)&" spivot="&Integer'Image(spivot(myid))&
               ", smax="&Real'Image(smax_pivot(myid)));

      Accept Barr3(KR : Integer) do
        K_Row := KR; -- pivot row, do not reduce
      end Barr3;
      Put_Line("slave "&Integer'Image(myid)&" reducing about pivot"
                           &Integer'Image(K_Row));
      --  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);
            Put_Line(Integer'Image(Myid)&" reduce B("&Integer'Image(I)&
                       ","&Integer'Image(J)&")="&Real'Image(B(I,J)));
          end loop;
        end if;
      end loop;
      --  finished inner reduction
      Put_Line(Integer'Image(Myid)&" finished reduction for K_Col="&
               Integer'Image(K_Col));


      Accept Barr4(K3 : Integer) do
        null; -- critical section available
      end Barr4;
      Put_Line("slave " & Integer'Image(myid) & " looping to next K_Col");
    end loop; -- N times, then terminate
    Put_Line("slave " & Integer'Image(myid) & " done ");
  end slave;

  G : array(1..NP) of slave; -- Start NP tasks running
                             -- task wait at Barr1
begin -- psimeq_dbg
  Put_Line("psimeq_dbg.adb running");
  if A'length(1)/=A'length(2) or A'length(1)/= Y'length then
    raise array_index_error ;
  end if ;

  Put_Line("psimeq_dbg debug starting");
  Put_Line("psimeq_dbg input data, n="&Integer'Image(N));
  for I in 1..N loop
    for J in 1..N loop
      Put_Line("A("&Integer'Image(I)&","&Integer'Image(J)&")="&
               Real'Image(A(I,J)));
    end loop;
    Put_Line("Y("&Integer'Image(I)&")="&Real'Image(Y(i)));
  end loop;

  -- set up row range for each task
  Srow(1) := 1;
  for I in 2..NP loop
    Srow(I) := Srow(I-1)+part;
    Put_Line("task "&Integer'Image(I-1)&" range="&Integer'Image(Srow(I-1))&
             ".."&Integer'Image(Srow(I)-1));
  end loop;
  Srow(NP+1) := N+1; -- may not be plus part
  Put_Line("task "&Integer'Image(NP)&" range="&Integer'Image(Srow(NP))&
           ".."&Integer'Image(Srow(NP+1)-1));

  -- 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).Barr1(T);
  end loop;

  for K_Col in 1..N loop  -- column to find max pivot
    Put_Line("task find their max pivot ");
    for T in 1..NP loop
      G(T).Barr2(K_Col);
    end loop;
    Put_Line("task found their max pivot ");

    -- 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;
    Put_Line("update_row I_pivot="&Integer'Image(I_Pivot)&
             ", abs_pivot="&Real'Image(Abs_Pivot));
    -- have pivot, record for X at end
    K_Row := I_pivot;
    prow(K_Col) := K_Row;
    frow(K_Row) := 1;
    Put_Line("update prow("&Integer'Image(K_Col)&")="&Integer'Image(K_Row));
    -- check for near singular
    if abs_pivot < 1.0E-10 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);
        Put_Line("pivrow B("&Integer'Image(K_Row)&","&Integer'Image(J)&
                 ")="&Real'Image(B(K_Row,J)));
      end loop;
      Put_Line("update reduced about pivot k_row="&Integer'Image(k_row));
    end if;

    Put_Line("task may now reduce their rows ");
    for t in 1..NP loop
      G(t).Barr3(K_row);
    end loop;


    Put_Line("task_loop starting stage 3 ");
    for t in 1..NP loop
      G(t).Barr4(t);
    end loop;

  end loop; -- K_Col

  Put_Line("master build X for return");
  --  build  X  for return, unscrambling rows
  for I in 1..N loop
    X(i) := B(Prow(i),N+1);
    Put_Line("prow("&Integer'Image(I)&")="&Integer'Image(prow(I)));
  end loop;
  return X;
end psimeq_dbg;