ell_relax: factor out the residual calculation kernel

1 files changed, 18 insertions, 11 deletions

diff --git a/ell_relax.c b/ell_relax.c
index 313af4e..aa791d8 100644
--- a/ell_relax.c
+++ b/ell_relax.c
@@ -103,6 +103,23 @@ derivatives_calc_s2(EllRelaxContext *ctx, double *dst, ptrdiff_t idx, ptrdiff_t
                                      / (144.0 * dx * dy);
 }
 
+static inline double residual_calc_kernel(EllRelaxContext *ctx, ptrdiff_t idx)
+{
+    EllRelaxInternal *priv = ctx->priv;
+
+    double u[ELL_RELAX_DIFF_COEFF_NB];
+    double res;
+
+    priv->derivatives_calc(ctx, u, idx, priv->stride);
+
+    res = -ctx->rhs[idx];
+    for (int i = 0; i < ARRAY_ELEMS(u); i++)
+        res += u[i] * ctx->diff_coeffs[i][idx];
+
+    return res;
+
+}
+
 static void residual_calc(EllRelaxContext *ctx)
 {
     EllRelaxInternal *priv = ctx->priv;
@@ -112,17 +129,7 @@ static void residual_calc(EllRelaxContext *ctx)
     for (int idx1 = 0; idx1 < ctx->domain_size[1]; idx1++)
         for (int idx0 = 0; idx0 < ctx->domain_size[0]; idx0++) {
             ptrdiff_t idx = idx1 * priv->stride + idx0;
-
-            double u[ELL_RELAX_DIFF_COEFF_NB];
-            double res;
-
-            priv->derivatives_calc(ctx, u, idx, priv->stride);
-
-            res = -ctx->rhs[idx];
-            for (int i = 0; i < ARRAY_ELEMS(u); i++)
-                res += u[i] * ctx->diff_coeffs[i][idx];
-
-            ctx->residual[idx] = res;
+            ctx->residual[idx] = residual_calc_kernel(ctx, idx);
         }
 
     for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) {