From d519d82a3e4b32944b77b1ae26cfefa45ec29d71 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Sat, 26 Jan 2019 10:19:30 +0100 Subject: ell_relax -> ell_grid_solve Generalize the API to allow for multiple solver types. This is done in preparation for the exact linear system inversion solver. --- mg2d.c | 64 ++++++++++++++++++++++++++++++++++++++-------------------------- 1 file changed, 38 insertions(+), 26 deletions(-) (limited to 'mg2d.c') diff --git a/mg2d.c b/mg2d.c index 3406030..9216af4 100644 --- a/mg2d.c +++ b/mg2d.c @@ -27,7 +27,7 @@ #include "common.h" #include "cpu.h" -#include "ell_relax.h" +#include "ell_grid_solve.h" #include "log.h" #include "mg2d.h" #include "mg2d_boundary.h" @@ -39,7 +39,7 @@ typedef struct MG2DLevel { unsigned int depth; - EllRelaxContext *solver; + EGSContext *solver; double *prolong_tmp; ptrdiff_t prolong_tmp_stride; @@ -86,8 +86,8 @@ static void log_relax_step(MG2DContext *ctx, MG2DLevel *level, const char *step_ prefix, level->depth, step_desc, res_old, res_new, res_old / res_new); } -static void mg_restrict_inject(EllRelaxContext *coarse, double *dst, ptrdiff_t dst_stride, - EllRelaxContext *fine, const double *src, ptrdiff_t src_stride) +static void mg_restrict_inject(EGSContext *coarse, double *dst, ptrdiff_t dst_stride, + EGSContext *fine, const double *src, ptrdiff_t src_stride) { for (size_t j = 0; j < coarse->domain_size[1]; j++) for (size_t i = 0; i < coarse->domain_size[0]; i++) { @@ -101,8 +101,8 @@ static void mg_restrict_inject(EllRelaxContext *coarse, double *dst, ptrdiff_t d } } -static void mg_restrict_fw(EllRelaxContext *coarse, double *dst, ptrdiff_t dst_stride, - EllRelaxContext *fine, const double *src, ptrdiff_t src_stride) +static void mg_restrict_fw(EGSContext *coarse, double *dst, ptrdiff_t dst_stride, + EGSContext *fine, const double *src, ptrdiff_t src_stride) { for (size_t j = 0; j < coarse->domain_size[1]; j++) for (size_t i = 0; i < coarse->domain_size[0]; i++) { @@ -117,8 +117,8 @@ static void mg_restrict_fw(EllRelaxContext *coarse, double *dst, ptrdiff_t dst_s } } -static void mg_prolongate(EllRelaxContext *fine, double *dst, ptrdiff_t dst_stride, - EllRelaxContext *coarse, const double *src, ptrdiff_t src_stride) +static void mg_prolongate(EGSContext *fine, double *dst, ptrdiff_t dst_stride, + EGSContext *coarse, const double *src, ptrdiff_t src_stride) { // for simplicity, this code writes one point over the domain size; // this allows us to avoid treating the boundary layers specially @@ -196,8 +196,8 @@ static int mg_interp_bilinear_kernel(void *arg, unsigned int job_idx, unsigned i } static void mg_interp_bilinear(TPContext *tp, - EllRelaxContext *ctx_dst, double *dst, ptrdiff_t dst_stride, - EllRelaxContext *ctx_src, const double *src, ptrdiff_t src_stride) + EGSContext *ctx_dst, double *dst, ptrdiff_t dst_stride, + EGSContext *ctx_src, const double *src, ptrdiff_t src_stride) { InterpParamsBilin ip; ip.dst_domain_size[0] = ctx_dst->domain_size[0]; @@ -228,8 +228,8 @@ static void coarse_correct_task(void *arg, unsigned int job_idx, unsigned int th static void restrict_residual(MG2DContext *ctx, MG2DLevel *dst, MG2DLevel *src) { - EllRelaxContext *s_src = src->solver; - EllRelaxContext *s_dst = dst->solver; + EGSContext *s_src = src->solver; + EGSContext *s_dst = dst->solver; if (s_src->domain_size[0] == 2 * (s_dst->domain_size[0] - 1) + 1) { mg_restrict_fw(s_dst, s_dst->rhs, s_dst->rhs_stride, s_src, s_src->residual, s_src->residual_stride); @@ -242,8 +242,8 @@ static void restrict_residual(MG2DContext *ctx, MG2DLevel *dst, MG2DLevel *src) static void prolong_solution(MG2DContext *ctx, MG2DLevel *dst, MG2DLevel *src) { - EllRelaxContext *s_src = src->solver; - EllRelaxContext *s_dst = dst->solver; + EGSContext *s_src = src->solver; + EGSContext *s_dst = dst->solver; if (s_dst->domain_size[0] == 2 * (s_src->domain_size[0] - 1) + 1) { mg_prolongate(s_dst, dst->prolong_tmp, dst->prolong_tmp_stride, s_src, s_src->u, s_src->u_stride); @@ -264,7 +264,7 @@ static int mg_relax_step(MG2DContext *ctx, MG2DLevel *level, const char *step_de start = gettime(); - ret = mg2di_ell_relax_step(level->solver); + ret = mg2di_egs_solve(level->solver); level->time_relax += gettime() - start; @@ -288,7 +288,7 @@ static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level) level->solver->domain_size[1]); /* re-init the current-level solver (re-calc the residual) */ - ret = mg2di_ell_relax_init(level->solver); + ret = mg2di_egs_init(level->solver); if (ret < 0) return ret; } @@ -344,7 +344,7 @@ static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level) /* re-init the current-level solver (re-calc the residual) */ res_prev = level->solver->residual_max; start = gettime(); - ret = mg2di_ell_relax_init(level->solver); + ret = mg2di_egs_init(level->solver); if (ret < 0) return ret; level->time_reinit += gettime() - start; @@ -393,8 +393,8 @@ static void bnd_copy(MG2DBoundary *bdst, const double *src, ptrdiff_t src_stride static void restrict_diff_coeffs(MG2DContext *ctx, MG2DLevel *dst, MG2DLevel *src) { - EllRelaxContext *s_src = src->solver; - EllRelaxContext *s_dst = dst->solver; + EGSContext *s_src = src->solver; + EGSContext *s_dst = dst->solver; if (s_src->domain_size[0] == 2 * (s_dst->domain_size[0] - 1) + 1) { for (int i = 0; i < ARRAY_ELEMS(s_src->diff_coeffs); i++) { mg_restrict_inject(s_dst, s_dst->diff_coeffs[i], s_dst->diff_coeffs_stride, @@ -453,7 +453,10 @@ static int mg_levels_init(MG2DContext *ctx) cur->solver->cpuflags = priv->cpuflags; cur->solver->fd_stencil = ctx->fd_stencil; - cur->solver->relax_factor = ctx->cfl_factor; + if (cur->solver->solver_type == EGS_SOLVER_RELAXATION) { + EGSRelaxContext *r = cur->solver->solver_data; + r->relax_factor = ctx->cfl_factor; + } prev = cur; cur = cur->child; @@ -484,7 +487,7 @@ int mg2d_solve(MG2DContext *ctx) { MG2DInternal *priv = ctx->priv; MG2DLevel *root = priv->root; - EllRelaxContext *s_root = root->solver; + EGSContext *s_root = root->solver; int64_t time_start; double res_prev, res_cur; int ret; @@ -501,7 +504,7 @@ int mg2d_solve(MG2DContext *ctx) if (ret < 0) return ret; - ret = mg2di_ell_relax_init(s_root); + ret = mg2di_egs_init(s_root); if (ret < 0) return ret; @@ -544,7 +547,7 @@ static void mg_level_free(MG2DLevel **plevel) return; free(level->prolong_tmp); - mg2di_ell_relax_free(&level->solver); + mg2di_egs_free(&level->solver); free(level); *plevel = NULL; @@ -564,7 +567,8 @@ static MG2DLevel *mg_level_alloc(const size_t domain_size) goto fail; level->prolong_tmp_stride = domain_size + 1; - level->solver = mg2di_ell_relax_alloc((size_t [2]){domain_size, domain_size}); + level->solver = mg2di_egs_alloc(EGS_SOLVER_RELAXATION, + (size_t [2]){domain_size, domain_size}); if (!level->solver) goto fail; @@ -723,9 +727,17 @@ void mg2d_print_stats(MG2DContext *ctx, const char *prefix) prefix, priv->count_solve, priv->time_solve / 1e6, priv->time_solve / 1e3 / priv->count_solve); while (level) { + EGSRelaxContext *r = NULL; int64_t level_total = level->time_relax + level->time_prolong + level->time_restrict + level->time_correct + level->time_reinit; - int64_t relax_total = level->solver->time_correct + level->solver->time_res_calc + level->solver->time_boundaries; + int64_t relax_total = (r ? r->time_correct : 0) + level->solver->time_res_calc + level->solver->time_boundaries; + + if (level->solver->solver_type == EGS_SOLVER_RELAXATION) + r = level->solver->solver_data; + + level_total = level->time_relax + level->time_prolong + level->time_restrict + + level->time_correct + level->time_reinit; + relax_total = (r ? r->time_correct : 0) + level->solver->time_res_calc + level->solver->time_boundaries; levels_total += level_total; @@ -742,7 +754,7 @@ void mg2d_print_stats(MG2DContext *ctx, const char *prefix) level->time_correct * 100.0 / level_total, level->time_reinit * 100.0 / level_total, level->solver->time_res_calc * 100.0 / relax_total, - level->solver->time_correct * 100.0 / relax_total, + (r ? r->time_correct : 0) * 100.0 / relax_total, level->solver->time_boundaries * 100.0 / relax_total ); -- cgit v1.2.3