----------------------------------------------------------------------
-- Aufgabe: partielle Differentialgleichung wie zB Laplace,
--                 diskretisiert ueber Gitter: 4*P(i,j) =   P(i-1,j) + P(i+1,j)
--                                                        + P(i,j-1) + P(i,j+1)
--                                                        - F(i,j)
-- Methode: tasks fuer Gitterpunkte, asynchrone Iteration, shared variables
--                 Hauptprogramm ueberwacht Konvergenz
-- Modifikation von Barnes95,p.443...:
--    Terminierung erst wenn Fehlerschranke simultan unterschritten
--
----------------------------------------------------------------------

with Stringpack; use Stringpack;
procedure Grid2 is
   N: constant := 5;
   subtype Full_Grid is Integer range 0 .. N;
   subtype Grid is Full_Grid range 1 .. N-1;   -- without border
   type Real is digits 7;  -- declare our own type
   type Matrix is array(Integer range <>, Integer range <>) of Real;
   pragma Atomic_Components(Matrix);
   P: Matrix(Full_Grid, Full_Grid)  := (others=>(others=>0.0));
   Delta_P: Matrix(Grid, Grid) := (others=>(others=>1.0));
   Tolerance: constant Real := 0.0001;
   Error_Limit: constant Real := Tolerance * Real(N-1)**2;
   Converged, NoChange: Boolean := False;
   pragma Atomic(Converged); pragma Atomic(NoChange);
   Error_Sum: Real;

   I : Positive := 1;  -- for counting interim output


   function F(I, J: Grid) return Real is   -- keep it simple
   begin return 0.0; end F;

   task type Iterator is
      entry Start(I, J: in Grid);
      entry Begin_Verify_Result;
      entry End_Verify_Result;
   end;

      Process: array (Grid, Grid) of Iterator;

      task body Iterator is
         I, J: Grid;
         New_P: Real;
      begin

         accept Start(I, J: in Grid) do
            Iterator.I := Start.I;
            Iterator.J := Start.J;
         end Start;

         loop
            New_P := 0.25 * (P(I-1, J) + P(I+1, J) + P(I, J-1)
			     + P(I, J+1) -  F(I, J));
            Delta_P(I, J) := New_P - P(I, J);
            P(I, J) := New_P;
            if Converged
	    then accept Begin_Verify_Result  do  -- do it again!
	       New_P := 0.25 * (P(I-1, J) + P(I+1, J) + P(I, J-1)
				+ P(I, J+1) -  F(I, J));
	       Delta_P(I, J) := New_P - P(I, J);
	       P(I, J) := New_P;
	    end Begin_Verify_Result;
	    accept End_Verify_Result;  -- synchro
	    exit when NoChange;
            end if;
         end loop;

      end Iterator;

begin            -- main subprogram, Iterator tasks now active

   for I in Full_Grid loop  --init P;
      P(I,0) := 1.0;
   end loop;

   for I in Grid loop
      for J in Grid loop
	 Process(I, J).Start(I, J);  -- tell them who they are
      end loop;
   end loop;

   loop
      Error_Sum := 0.0;
      for I in Grid loop
	 for J in Grid loop
	    Error_Sum := Error_Sum + Delta_P(I, J)**2;
	 end loop;
      end loop;
      Converged := Error_Sum < Error_Limit;
      if Converged
      then                   -- make sure that error bounds hold simultaneously
	 Error_Sum := 0.0;
	 for I in Grid loop
	    for J in Grid loop
	       Process(I,J).Begin_Verify_Result;
	       Error_Sum := Error_Sum + Delta_P(I, J)**2;
	    end loop;
	 end loop;
	 Converged := False;
	 NoChange := Error_Sum < Error_Limit;
	 for I in Grid loop
	    for J in Grid loop
	       Process(I,J).End_Verify_Result;  -- synchro
	    end loop;
	 end loop;
      end if;
      if NoChange
      then exit;
      else
	 if I < 5 then
	    I := I + 1;
	    Print("--------------------------------------------------");
	    for I in Full_Grid loop        -- output state
	       Outbuffer := Varnull;
	       for J in Full_Grid loop
		  Outbuffer := Outbuffer & Cvfs(Float(P(I,J)),1,3) & " ";
	       end loop;
	       Print;
	    end loop;
	 end if;
      end if;
   end loop;

   Print("results: ----------------------------------------------");
   for I in Full_Grid loop        -- output results
      Outbuffer := Varnull;
      for J in Full_Grid loop
	 Outbuffer := Outbuffer & Cvfs(Float(P(I,J)),1,3) & " ";
      end loop;
      Print;
   end loop;


end Grid2;