#include <stdlib.h>
# define N 800            /* for 2-d arrays with indices 0-250 in both dirs */
# define K N*N          /* for 1-d array to hold above 2-d array */
# define KR (K+1)/2 /* for 1-d array of red fd blocks */
# define KB K/2         /* for 1-d array of black fd blocks */

main()
{
int i,j,m,k=0;
float u[N][N];
float r[KR],b[KB],v[K];
int label[N][N],lab[K],rl[KR],bl[KB];
int rn,rs,bn,bs;

rn = (N+1)/2;
rs = N/2;
bn = rs;
bs = rn;

for(j=1;j<N-1;j++)      /* set interior values of label and u arrays */
{
for(i=1;i<N-1;i++)
   {
   label[i][j] = 1;
   u[i][j] = 7;
   }
}
for(j=0;j<N;j++) {u[j][0] = 15; label[j][0] = 0;}   /* set boundary values   */
for(j=0;j<N;j++) {u[0][j] = 15; label[0][j] = 0;}   /* in label and u arrays */
for(j=0;j<N;j++) {u[N-1][j] = 1; label[N-1][j] = 0;}
for(j=0;j<N;j++) {u[j][N-1] = 1; label[j][N-1] = 0;}

for(j=0;j<N;j++)
{
   for(i=0;i<N;i++)
   {
   lab[k] = label[i][j] ;    /* write 2-d label array into 1-d array */
   v[k] = u[i][j];               /* write 2-d u array into 1-d array */
   k++;
   }
}

for(k=0;k<KR;k++)
{
rl[k] = lab[2*k];         /* identify red cells in 1-d label array */
r[k] = v[2*k];              /* identify red cells in 1-d u array */
}
for(k=0;k<KB;k++)
{
bl[k] = lab[2*k+1];       /* identify black cells in 1-d label array */
b[k] = v[2*k+1];            /* identify black cells in 1-d u array */
}
k = 0;

while(k < 50) /* begin red-black iterations */
{
k++;

for(i=rn;i< KR-rn;i++) /* update all red cells using rl to mask boundary */
{
r[i] = r[i] + .25*(b[i-1]+b[i]-4*r[i]+b[i-rn]+b[i+rs])*rl[i];
}

for(i=0;i<KR;i++)
{
v[2*i] = r[i]; /* store red values in 1-d u-array */
}

m=0;
for(j=0;j<N;j++)
{
   for(i=0;i<N;i++)
   {
   u[i][j] = v[m];            /* write 1-d u-array into 2-d u-array */
   m++;                          /* in preparation for plotting red update */
   }
}

for(i=bs;i<KB-bs;i++) /* update all black cells using bl to mask boundary */
{
b[i] = b[i] + .25*(r[i+1]+r[i]-4*b[i]+r[i-bn]+r[i+bs])*bl[i];
}

for(i=0;i<KB;i++)
{
v[2*i+1]= b[i]; /* store black values in 1-d u-array */
}

m=0;
for(j=0;j<N;j++)
{
   for(i=0;i<N;i++)
   {                                    /* write 1-d u-array into 2-d u-array */
   u[i][j] = v[m];              /* in preparation for plotting black cells*/
   m++;
   }
}
/* printf("%f\n",u[20][20]);   */

}
return(0);
}