SYSTEM  Dirichlet;
(* Compute the Dirichlet function on a matrix.*)
(* David Kotz 1991 *)
CONST N = 50;   (* remember that P=N^2 *)
TYPE  matrix = ARRAY [N],[N] OF REAL;

CONFIGURATION  grid [N],[N];
CONNECTION  	 	    (* non-wrapped mesh *)
               left:  grid[i,j]  <->  grid[i, j-1].right;
               up:    grid[i,j]  <->  grid[i-1, j].down;

SCALAR iter         : INTEGER;
       Niter : integer;
       printmatrix : boolean;
       A : matrix;
VECTOR val : REAL;                  (* grid version of matrix A and then B *)
      
(* One iteration of the Dirchlet function *)
PROCEDURE Iterate;
VECTOR vleft,vright,vup,vdown : REAL;
BEGIN
  PARALLEL             (* can't just add  [1..N-2], [1..N-2] specifier  here *)
    PROPAGATE.left(val, vright);
    PROPAGATE.right(val, vleft);
    PROPAGATE.up(val, vdown);
    PROPAGATE.down(val, vup);
    IF (0 < dim1 < N-1)  AND  (0 < dim2 < N-1)         (* exclude boundaries *)
    THEN val := (vleft + vright + vup + vdown) / 4.;
    END;
  ENDPARALLEL;
  IF printmatrix THEN STORE(val,A); END;
END Iterate;

(* Print the matrix out *)
PROCEDURE PrintOut;
SCALAR i,j: CARDINAL;
BEGIN
  IF printmatrix THEN
    FOR i:=0 TO N-1 DO
      FOR j:=0 TO N-1 DO WriteFixPt(A[i,j], 10,2) END;
      WriteLn;
    END;
    WriteLn;
  END;
END PrintOut;

(* Initialize the matrix *)
PROCEDURE Init;
BEGIN
  PARALLEL
    IF (0 < dim1 < N-1)  AND  (0 < dim2 < N-1)
    THEN val := 0.0;		(* internal points are zero *)
    ELSE	  (* one case for each edge *)
    	   IF (DIM1 = 0) THEN val := float(DIM2+1);
    	   ELSIF (DIM2 = 0) THEN val := float(DIM1+1);
    	   ELSIF (DIM1 = N-1) THEN val := float(N-DIM2);
    	   ELSIF (DIM2 = N-1) THEN val := float(N-DIM1);
    	   END;
    END;
  ENDPARALLEL;
END Init;

BEGIN
  WriteString("Number of iterations: ");
  ReadInt(Niter);
  WriteString("Print matrix? ");
  ReadBool(printmatrix);

  Init;
  IF printmatrix THEN PrintOut; END;

  FOR iter := 1 TO Niter DO
      Iterate;
      IF printmatrix THEN PrintOut; END;
  END;
END Dirichlet.