global size
  int N = 8    # grid dimension, default 8
body size
  getarg(1, N)
end

global data
  import size
  real grid[2][0:N+1,0:N+1]
  real epsilon = .3
  int nxt = 2
  int cur = 1
  op print()
body data
  grid = ([2] ([N+2] ([N+2] 0.0)))
  real l=1.0
  real t=2.1
  real r=3.3
  real b= 4.9
  for [ i = 0 to N+1 ] {
    grid[cur][i,0] = l; grid[nxt][i,0] = l
  }
  for [ i = 0 to N+1 ] {
    grid[cur][0,i] = t; grid[nxt][0,i] = t
  }
  for [ i = 0 to N+1 ] {
    grid[cur][i,N+1] = r; grid[nxt][i,N+1] = r
  }
  for [ i = 0 to N+1 ] {
    grid[cur][N+1,i] = b; grid[nxt][N+1,i] = b
  }

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

resource jacobi()
  import data, size
  real diff[1:N,1:N]    # differences

  procedure update(int i, int j)  {
    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[i,j] = abs(grid[nxt][i,j] - grid[cur][i,j])
  }

  procedure check_diffs(int i)  {
    for [ j = 2 to N ] {
      diff[i,1] = max(diff[i,1], diff[i,j])
    }
  }

  int iters = 0
  real maxdiff
  while (true) {    # iterate until convergence
    iters++
    co [ i = 1 to N, j = 1 to N ] update(i,j) oc
    nxt :=: cur    # swap roles of grids
    maxdiff = 0.0
    co [ i = 1 to N ] check_diffs(i) ->
	  maxdiff = max(maxdiff, diff[i,1])
    oc
    if (maxdiff <= epsilon) { exit }
  }
  write("convergence after", iters, "iterations")
  print()
end