SYSTEM  systolic_array;
(* compute the matrix product "c := a * b" *)
CONST max    = 5;
TYPE  matrix = ARRAY [1..max],[1..max] OF REAL;

CONFIGURATION  grid [max],[max];
CONNECTION     left:   grid[i,j]  ->  grid[i,(j-1) MOD max].left;
               up:     grid[i,j]  ->  grid[(i-1) MOD max,j].up;
               shiftA: grid[i,j]  ->  grid[i,(j-i) MOD max].shiftA;
               shiftB: grid[i,j]  ->  grid[(i-j) MOD max,j].shiftB;

SCALAR i,j         : INTEGER;
       a,b,c       : matrix;


PROCEDURE matrix_mult(SCALAR VAR a,b,c : matrix);
(* c := a * b *)
SCALAR k: INTEGER;
VECTOR ra,rb,rc : REAL;
BEGIN
  LOAD (ra,a);
  LOAD (rb,b);
  PARALLEL
    PROPAGATE.shiftA(ra);
    PROPAGATE.shiftB(rb);
    rc := ra * rb;
    FOR k := 2 TO max DO
      PROPAGATE.left(ra);
      PROPAGATE.up(rb);
      rc := rc + ra * rb;
    END;
  ENDPARALLEL;
  STORE(rc,c);
END matrix_mult;


PROCEDURE out(SCALAR VAR a: matrix);
SCALAR i,j: CARDINAL;
BEGIN
  FOR i:=1 TO max DO
    FOR j:=1 TO max DO WriteFixPt(a[i,j], 10,2) END;
    WriteLn
  END;
  WriteLn
END out;

BEGIN
  (* preset input matrices "a" and "b" (or read them from a file instead) *)
  FOR i:=1 TO max DO
    FOR j:=1 TO max DO
      a[i,j]:=FLOAT(10*i+j);
      b[i,j]:=FLOAT(i+j);
   END;
  END;
  out(a); out(b);
  matrix_mult(a,b,c);
  out(c);                (* print result matrix "c" *)
END systolic_array.