sweep1d.c

#include "array.h"
#include "poisson.h"


double
sweep1d(dmatrix a, dmatrix f, dmatrix b, Decomposition oned)
{
  /* jacobi sweep through the local array */

   int             i, j;
   double          error = 0.0;
   double          resid = 0.0;


   for (i = 1; i <= oned.lx; i++)
   {
      for (j = 1; j <= oned.ly; j++)
      {
         resid = INDEX(a, i-1, j) + INDEX(a, i, j+1) 
            + INDEX(a, i+1, j) + INDEX(a, i, j-1) - 4*INDEX(a, i, j)
               - INDEX(f, i, j);
         error += fabs(resid);
         INDEX(b, i, j) = INDEX(a, i, j) + 0.25*resid;
      }
   }

   return(error);
}