----------------------------------------------------------------------
-- 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
-- Quelle:    Barnes95,p.443...
--                                   !!! problematic termination method !!!!
----------------------------------------------------------------------

with Stringpack; use Stringpack;
procedure Grid is
   N: constant := 5;
   subtype Full_Grid is Integer range 0 .. N;  -- grid with border
   subtype Grid is Full_Grid range 1 .. N-1;   -- grid 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);  -- shared!
   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: Boolean := False;
   pragma Atomic(Converged);  -- shared!
   Error_Sum: Real;


   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);
   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;
         exit when Converged;
      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  -- Konvergenzprüfung
      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;  -- problematic method !!!!!!!!!!!!
      exit when Converged;
   end loop;

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