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.

1 files changed, 38 insertions, 26 deletions

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
                   );