ell_relax: do not calculate the residual at the fixed-val boundaries

Keep it at zero there to begin with, rather than enforce this later.
This will be useful in the following commits.

1 files changed, 18 insertions, 23 deletions

diff --git a/ell_relax.c b/ell_relax.c
index 042c4a5..a616e2e 100644
--- a/ell_relax.c
+++ b/ell_relax.c
@@ -47,6 +47,8 @@ struct EllRelaxInternal {
                                ptrdiff_t stride, const double *u, const double *rhs,
                                const double * const diff_coeffs[ELL_RELAX_DIFF_COEFF_NB],
                                const double *fd_factors);
+    ptrdiff_t residual_calc_offset;
+    size_t residual_calc_size[2];
     double fd_factors[ELL_RELAX_DIFF_COEFF_NB];
     double relax_factor;
 
@@ -205,13 +207,13 @@ static void residual_calc_task(void *arg, unsigned int job_idx, unsigned int thr
 {
     EllRelaxContext  *ctx = arg;
     EllRelaxInternal *priv = ctx->priv;
-    const ptrdiff_t offset = job_idx * priv->stride;
+    const ptrdiff_t offset = priv->residual_calc_offset + job_idx * priv->stride;
     const double *diff_coeffs[ELL_RELAX_DIFF_COEFF_NB];
 
     for (int i = 0; i < ARRAY_ELEMS(diff_coeffs); i++)
         diff_coeffs[i] = ctx->diff_coeffs[i] + offset;
 
-    priv->residual_calc_line(ctx->domain_size[0], ctx->residual + offset,
+    priv->residual_calc_line(priv->residual_calc_size[0], ctx->residual + offset,
                              priv->stride, ctx->u + offset, ctx->rhs + offset,
                              diff_coeffs, priv->fd_factors);
 }
@@ -223,25 +225,7 @@ static void residual_calc(EllRelaxContext *ctx)
 
     start = gettime();
 
-    tp_execute(ctx->tp, ctx->domain_size[1], residual_calc_task, ctx);
-
-    for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) {
-        const ptrdiff_t strides[2] = { 1, priv->stride };
-        const int si = boundary_def[i].stride_idx;
-        const ptrdiff_t stride_boundary = strides[si];
-        const ptrdiff_t stride_offset   = strides[!si];
-        const size_t size_boundary      = ctx->domain_size[si];
-        const size_t size_offset        = ctx->domain_size[!si];
-
-        double *res = ctx->residual + boundary_def[i].is_upper * ((size_offset - 1) * stride_offset);
-
-        if (ctx->boundaries[i].type == ELL_RELAX_BC_TYPE_FIXVAL) {
-            for (size_t j = 0; j < size_boundary; j++) {
-                *res = 0.0;
-                res += stride_boundary;
-            }
-        }
-    }
+    tp_execute(ctx->tp, priv->residual_calc_size[1], residual_calc_task, ctx);
 
     ctx->time_res_calc += gettime() - start;
     ctx->count_res++;
@@ -421,13 +405,24 @@ int mg2di_ell_relax_init(EllRelaxContext *ctx)
         ctx->tp = priv->tp_internal;
     }
 
+    priv->residual_calc_offset  = 0;
+    priv->residual_calc_size[0] = ctx->domain_size[0];
+    priv->residual_calc_size[1] = ctx->domain_size[1];
+    if (ctx->boundaries[ELL_RELAX_BOUNDARY_0L].type == ELL_RELAX_BC_TYPE_FIXVAL)
+        priv->residual_calc_offset += ctx->residual_stride;
+    if (ctx->boundaries[ELL_RELAX_BOUNDARY_1L].type == ELL_RELAX_BC_TYPE_FIXVAL)
+        priv->residual_calc_offset++;
+
     for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) {
-        if (ctx->boundaries[i].type == ELL_RELAX_BC_TYPE_FIXDIFF) {
+        EllRelaxBoundary *bnd = &ctx->boundaries[i];
+        if (bnd->type == ELL_RELAX_BC_TYPE_FIXDIFF) {
            for (int k = 0; k < ctx->domain_size[boundary_def[i].stride_idx]; k++)
-               if (ctx->boundaries[i].val[k] != 0.0) {
+               if (bnd->val[k] != 0.0) {
                    mg2di_log(&ctx->logger, 0, "Only zero boundary derivative supported\n");
                    return -ENOSYS;
                }
+        } else if (bnd->type == ELL_RELAX_BC_TYPE_FIXVAL) {
+            priv->residual_calc_size[!boundary_def[i].stride_idx]--;
         }
     }