-- psimeq.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 real_arrays; use real_arrays;
-- Ada.Numerics.Real_Arrays where available
with Ada.Text_IO;
use  Ada.Text_IO;

function psimeq (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
  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
    Accept Barr1(Me : Integer) do
      myid := Me;
    end Barr1;
    -- 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;
      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; -- give Master time to reduce about pivot

      Accept Barr3(KR : Integer) do
        K_Row := KR; -- pivot row, do not reduce
      end Barr3;
      --  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
    Accept Barr4(K3 : Integer) do
      null; -- critical section available
    end Barr4;
  end slave;

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

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

  -- 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
    -- task find their max pivot
    for T in 1..NP loop
      G(T).Barr2(K_Col);
    end loop;

    -- 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;

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

    -- task may loop to find next pivot
  end loop; -- K_Col
  for t in 1..NP loop -- be sure all tasks are finished
    G(t).Barr4(t);
  end loop;

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