! Red-black update for Poisson's equation
! =======================================

USE gnuplot   ! prepare for plotting
INTEGER nx,ny   ! declare two INTEGER variables
ASK nx,ny   ! ask array dimensions from the console user
ARRAY(1..nx,1..ny) OF REAL var,rhs   ! variable ARRAY dimensions are allowed

var=0; rhs=0
DO var(i,HI)=1 FOR ALL i   ! assign some boundary condition

LOOP FOR 100 TIMES   ! main iteration loop
  LOOP FOR parity=0 TO 1 AND i=LO+1 TO HI-1 AND j=LO+1 +
    (i+LO+1+parity) MOD 2 TO HI-1 BY 2
    ! this traverses the array in checkerboard fashion
    var(i,j)=[var(i,j+1)+var(i,j-1)+var(i+1,j)+var(i-1,j)+rhs(i,j)]/4
    ! this statement is unchanged from Gauss-Seidel
  REPEAT LOOP
REPEAT   ! "LOOP" is optional here
  
SPLOT var
READ ! wait on input and keep plot displayed until Enter is pressed

!(
  Variables are allowed as ARRAY dimensions. The bounds OF LOOP indices are
  automatically extracted from arrays the indices are used in, AND made
  available through the predefined LO AND HI variables.
!)