diff options
| -rw-r--r-- | ell_relax.c | 29 |
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++) { |