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

global data
  import size
  real grid[0:N+1,0:N+1] = ([N+2] ([N+2] 0.0))
  op print()
body data
  real  left, top, right, bottom
  left = 1.0; top = 2.1; right = 3.3; bottom = 4.9
  for [ i = 0 to N+1 ] { grid[i,0] = left }
  for [ i = 0 to N+1 ] { grid[0,i] = top }
  for [ i = 0 to N+1 ] { grid[i,N+1] = right }
  for [ i = 0 to N+1 ] { grid[N+1,i] = bottom }

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

resource gauss_seidel()
  import data, size
  real  diff
  real epsilon = .3            # convergence

  procedure update()  {
    real  newgrid
    diff = 0.0
    for [ i = 1 to N, j = 1 to N ] {
      newgrid = (grid[i-1,j] + grid[i,j-1] + grid[i+1,j] + grid[i,j+1]) / 4
      diff = max(diff, abs(newgrid - grid[i,j]))
      grid[i,j] = newgrid
    }
  }

  int iters = 0
  while (true) {
    iters++
    update()
    if (diff <= epsilon) { exit }
  }
  write("convergence after", iters, "iterations")
  print()
end