#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); }