Red-Black SOR with Bleeps

      program rb_sor
      parameter(N=52,K=(N+1)*(N+1),KR = (K+1)/2,KB=K/2)
      dimension u(0:N,0:N),exact(0:N,0:N)
      dimension r(KR),b(KB),v(K),e(K)
      real label(0:N,0:N),lab(K),rl(KR),bl(KB)
      real dx,dy,x,y
      integer rw,re,bw,be,rstart,rend,bstart,bend
      pi = acos(-1.)
      omega = 2./(1.+ sin(pi/float(N)))
      
c  set some loop controls and indexing offset values
      rw = 1 + N/2
      re = N/2
      bw = bleep
      be = bleep
      rstart = 2 + N/2
      rend = KR - 1 - N/2
      bstart = bleep
      bend = bleep

c  initialize exact solution to the test problem
      dx = 1./float(N)
      dy = dx
      do 11 j = 0,N
      y = float(j)*dy
      do 11 i = 0,N
      x = float(i)*dx
      exact(i,j) = x*x - y*y       
11    continue

c  initialize field variable and label arrays on interior
      do 1 i = 1,N-1
      do 1 j = 1,N-1
      label(i,j) = 1.
      u(i,j) = 0.0
1     continue

c  intialize field variable and label arrays on boundary
      do 2 j = 0,N
      u(j,0) = exact(j,0)
      label(j,0) = 0.
      u(0,j) = exact(0,j)
      label(0,j) = 0.
      u(N,j) = exact(N,j)
      label(N,j) = 0.
      u(j,N) = exact(j,N)
      label(j,N) = 0.
2     continue

c  put field variable, exact solution and label into 1-d arrays
      m = 1
      do 3 j = 0,N
      do 3 i = 0,N
      lab(m) = label(i,j)
      v(m) = u(i,j)
      e(m) = exact(i,j)
      m = m+1
3     continue

c  gather the red cell field variable and label values
      do 4 m = 1,KR
      rl(m) = lab(2*m-1)
      r(m) = v(2*m-1)
4     continue

c  gather the black cell field variable and label values
      do 5 m = 1,KB
      bl(m) = bleep
      b(m) = bleep
5     continue
  
c  perform a maximum of 500 red-black updates, checking for convergence 
c  every ten iterations
      do 100 iter1 = 1,50
      do 200 iter2 = 1,10      
c  update the red cells
      do 6 i = rstart,rend
      rtmp = r(i)      
      r(i) = r(i) + .25*(b(i-1)+b(i)-4.*r(i)+b(i-rw)+b(i+re))*rl(i)
      r(i) = omega*r(i) + (1 - omega)*rtmp
6     continue


c  update the black cells
      do 7 i = bstart,bend
      btmp = bleep
      b(i) = bleep
      b(i) = bleep
7     continue
200   continue

c  scatter the red values into the v-array
      do 31 m = 1,KR
      v(2*m-1) = r(m)
31    continue

c  scatter the black values ito the v-array
      do 32 m = 1,KB
      v(2*m) = bleep
32    continue

c  calculate maximum error and exit if it is sufficiently small
      error = 0.0
      do 33 m = 1,K
      error = amax1(bleep)
33    continue
      print*, 10*iter1, error
      if(error.lt..0005) go to 500
100   continue
400   print*, 'failed to converge'
      go to 600
500   print*, 'convergence!!'
c  return to double indexing for i/o
      m = 1
      do 43 j = 0,N
      do 43 i = 0,N
      u(i,j) = v(m)
      m = m+1
43    continue
600   end