Jacobi

/*
 * Blocked Jacobi solver: one iteration step
 */
double relax_jacobi (double *u, double *utmp, unsigned sizex, unsigned sizey)
{
    double sum=0.0;
    int nbx, bx, nby, by;

    nbx = NB;
    bx = sizex/nbx;
    nby = NB;
    by = sizey/nby;
    for (int ii=0; ii<nbx; ii++){
        double sum2 = 0.0;
        #pragma omp parallel for schedule(guided, 2) firstprivate(sum2)
        for (int jj=0; jj<nby; jj++){
           for (int i=1+ii*bx; i<=min((ii+1)*bx, sizex-2); i++){
                for (int j=1+jj*by; j<=min((jj+1)*by, sizey-2); j++) {
                   double diff=0.0;
                   utmp[i*sizey+j]= 0.25 * (u[ i*sizey     + (j-1) ]+  // left
                                             u[ i*sizey     + (j+1) ]+  // right
                                             u[ (i-1)*sizey + j     ]+  // top
                                             u[ (i+1)*sizey + j     ]); // bottom
                    diff = utmp[i*sizey+j] - u[i*sizey + j];
                    sum2 += diff * diff;
                }
           }
           #pragma omp critical
           sum+=sum2;
        }
    }

    return sum;
}