From d0bce68cfe7f45fc417fb197772c03ebba4af902 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Fri, 11 Jan 2019 15:48:05 +0100 Subject: 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. --- ell_relax.c | 41 ++++++++++++++++++----------------------- 1 file 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]--; } } -- cgit v1.2.3