resource jacobi
  import main
  op toprow([*] real)
  op bottomrow([*] real)
  op compute(cap jacobi up, cap jacobi down, cap main coord)
body jacobi(int id, int N, int PR, int S, real l, real t, real r, real b)
  real grid[2][0:S+1,0:N+1]
  int cur = 1
  int nxt = 2
  real diff    # max difference

  # initialize grids
  grid = ([2] ([S+2] ([N+2] 0.0)))
  for [ i = 0 to S+1 ] {  # left and right
    grid[cur][i,0] = l; grid[nxt][i,0] = l
    grid[cur][i,N+1] = r; grid[nxt][i,N+1] = r
  }
  if (id == 1) {  # top for process 1
    for [ j = 0 to N+1 ] {
      grid[cur][0,j] = t; grid[nxt][0,j] = t
    }
  }
  if (id == PR) {  # bottom for process PR
    for [ j = 0 to N+1 ] {
      grid[cur][S+1,j] = b; grid[nxt][S+1,j] = b
    }
  }

  proc compute(up, down, coord)  {
    while (true) {
      # compute new values for grid points
      diff = 0.0
      for [ i = 1 to S, j = 1 to N ] {
	grid[nxt][i,j] = (grid[cur][i-1,j] + grid[cur][i,j-1] + \
			   grid[cur][i+1,j] + grid[cur][i,j+1]) / 4
	diff = max(diff, abs(grid[nxt][i,j] - grid[cur][i,j]))
      }
      # replace old values by new ones, and
      # exchange top and bottom rows with neighbors
      cur :=: nxt
      send up.bottomrow(grid[cur][1,*])
      send down.toprow(grid[cur][S,*])
      if (id != 1) { receive toprow(grid[cur][0,*]) }
      if (id != PR) {
	receive bottomrow(grid[cur][S+1,*])
      }
      # check for termination
      if (coord.terminate(diff)) { exit }
    }
  }

  final {
    for [ i = 1 to S ] {
      for [ j = 1 to N ] { printf("%5.3f ",grid[cur][i,j]) }
      write()
    }
  }
end

resource main
  import jacobi
  op terminate(real diff) returns bool  ans
body main()
  int N = 8    # grid dimension
  int PR = 2   # number of processes
  int S        # strip size
  # read command-line arguments, if present
  getarg(1, N); getarg(2, PR); S = N/PR
  if ((N mod PR) != 0) {
    write("N must be a multiple of PR"); stop(1)
  }
  # read boundary values and epsilon
  real l = 1.0
  real t = 2.1
  real r = 3.3
  real b = 4.9
  real epsilon = 0.3

  # create instances of jacobi()
  cap jacobi  jcap[1:PR]
  if (PR == 1) {
      jcap[1] = create jacobi(1,N,PR,S,l,t,r,b)
  } else {
      jcap[1] = create jacobi(1,N,PR,S,l,t,r,0.0)
      for [ i = 2 to PR-1 ] {
	jcap[i] = create jacobi(i,N,PR,S,l,0.0,r,0.0)
      }
      jcap[PR] = create jacobi(PR,N,PR,S,l,0.0,r,b)
  }

  # start the computation
  if (PR == 1) {
      send jcap[1].compute(noop,noop,myresource())
  } else {
      send jcap[1].compute(noop,jcap[2],myresource())
      for [ i = 2 to PR-1 ] {
	send jcap[i].compute(jcap[i-1],
			     jcap[i+1],myresource())
      }
      send jcap[PR].compute(jcap[PR-1],
			    noop,myresource())
  }

  # do termination checks until convergence
  int iters = 0
  while (true) {
    iters++
    # wait for all processes to call terminate, then
    # service invocation with largest value for diff
    in terminate(diff) returns ans st ?terminate == PR by -diff ->
	   ans = diff <= epsilon
	   for [ i = 1 to PR-1 ] {
	     in terminate(diff2) returns ans2 ->
	       ans2 = ans
	     ni
	   }
	   if (ans) { exit }
    ni
  }

  final {  # print results
    write("convergence after", iters, "iterations\n")
    for [ i = 1 to PR ] { destroy jcap[i] }
  }
end