Add a function to measure timings for benchmarking.

4 files changed, 105 insertions, 1 deletions

diff --git a/ell_relax.c b/ell_relax.c
index 96038e8..d79921a 100644
--- a/ell_relax.c
+++ b/ell_relax.c
@@ -100,7 +100,9 @@ derivatives_calc_s2(EllRelaxContext *ctx, double *dst, ptrdiff_t idx, ptrdiff_t
 static void residual_calc(EllRelaxContext *ctx)
 {
     EllRelaxInternal *priv = ctx->priv;
+    int64_t start;
 
+    start = gettime();
     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;
@@ -134,6 +136,9 @@ static void residual_calc(EllRelaxContext *ctx)
             }
         }
     }
+
+    ctx->time_res_calc += gettime() - start;
+    ctx->count_res++;
 }
 
 static void boundaries_apply_fixval(double *dst, ptrdiff_t dst_stride,
@@ -172,7 +177,9 @@ static void boundaries_apply(EllRelaxContext *ctx)
 {
     EllRelaxInternal *priv = ctx->priv;
     const ptrdiff_t strides[2] = { 1, priv->stride };
+    int64_t start;
 
+    start = gettime();
     for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) {
         const int si = boundary_def[i].stride_idx;
         const ptrdiff_t stride_boundary = strides[si];
@@ -227,19 +234,25 @@ static void boundaries_apply(EllRelaxContext *ctx)
                 dst[-idx] = dst[idx];
             }
     }
+    ctx->time_boundaries += gettime() - start;
+    ctx->count_boundaries++;
 }
 
 int mg2di_ell_relax_step(EllRelaxContext *ctx)
 {
     EllRelaxInternal *priv = ctx->priv;
     const double cfl_fac = (1.0 / 5.0) * ctx->step[0] * ctx->step[1];
+    int64_t start;
 
+    start = gettime();
     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;
 
             ctx->u[idx] += cfl_fac * ctx->residual[idx];
         }
+    ctx->time_correct += gettime() - start;
+    ctx->count_correct++;
 
     boundaries_apply(ctx);
     residual_calc(ctx);
diff --git a/ell_relax.h b/ell_relax.h
index 74355e6..0df42d5 100644
--- a/ell_relax.h
+++ b/ell_relax.h
@@ -222,6 +222,15 @@ typedef struct EllRelaxContext {
      * Distance between neighbouring gridpoints along coord1.
      */
     ptrdiff_t diff_coeffs_stride;
+
+    /* timings */
+    int64_t count_correct;
+    int64_t time_correct;
+
+    int64_t time_boundaries;
+    int64_t count_boundaries;
+    int64_t time_res_calc;
+    int64_t count_res;
 } EllRelaxContext;
 
 /**
diff --git a/mg2d.c b/mg2d.c
index 668e8ef..9d98ec9 100644
--- a/mg2d.c
+++ b/mg2d.c
@@ -39,12 +39,23 @@ typedef struct MG2DLevel {
     ptrdiff_t prolong_tmp_stride;
 
     struct MG2DLevel *child;
+
+    /* timings */
+    int64_t count_cycles;
+    int64_t time_relax;
+    int64_t time_prolong;
+    int64_t time_restrict;
+    int64_t time_correct;
+    int64_t time_reinit;
 } MG2DLevel;
 
 struct MG2DInternal {
     MG2DLogger logger;
 
     MG2DLevel *root;
+
+    int64_t time_solve;
+    int64_t count_solve;
 };
 
 static double findmax(double *arr, size_t len)
@@ -176,23 +187,34 @@ static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level)
 
     /* handle coarsest grid */
     if (!level->child) {
+        int64_t start = gettime();
+
         for (int j = 0; j < 16; j++) {
             ret = mg2di_ell_relax_step(level->solver);
             if (ret < 0)
                 return ret;
         }
+
+        level->time_relax += gettime() - start;
+        level->count_cycles++;
+
         goto finish;
     }
 
     for (int i = 0; i < ctx->nb_cycles; i++) {
+        int64_t start;
+
         /* pre-restrict relaxation */
+        start = gettime();
         for (int j = 0; j < ctx->nb_relax_pre; j++) {
             ret = mg2di_ell_relax_step(level->solver);
             if (ret < 0)
                 return ret;
         }
+        level->time_relax += gettime() - start;
 
         /* restrict the residual as to the coarser-level rhs */
+        start = gettime();
         if (!is_power2(level->solver->domain_size[0] - 1)) {
             mg_interp_bilinear(level->child->solver, level->child->solver->rhs, level->child->solver->rhs_stride,
                                level->solver, level->solver->residual, level->solver->residual_stride);
@@ -200,7 +222,7 @@ static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level)
             mg_restrict_fw(level->child->solver, level->child->solver->rhs, level->child->solver->rhs_stride,
                            level->solver, level->solver->residual, level->solver->residual_stride);
         }
-
+        level->time_restrict += gettime() - start;
 
         /* solve on the coarser level */
         ret = mg_solve_subgrid(ctx, level->child);
@@ -208,6 +230,7 @@ static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level)
             return ret;
 
         /* prolongate the coarser-level correction */
+        start = gettime();
         if (!is_power2(level->solver->domain_size[0] - 1)) {
             mg_interp_bilinear(level->solver, level->prolong_tmp, level->prolong_tmp_stride,
                                level->child->solver, level->child->solver->u, level->child->solver->u_stride);
@@ -215,26 +238,35 @@ static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level)
             mg_prolongate(level->solver, level->prolong_tmp, level->prolong_tmp_stride,
                           level->child->solver, level->child->solver->u, level->child->solver->u_stride);
         }
+        level->time_prolong += gettime() - start;
 
         /* apply the correction */
+        start = gettime();
         for (size_t idx1 = 0; idx1 < level->solver->domain_size[1]; idx1++)
             for (size_t idx0 = 0; idx0 < level->solver->domain_size[0]; idx0++) {
                 const ptrdiff_t idx_dst = idx1 * level->solver->u_stride + idx0;
                 const ptrdiff_t idx_src = idx1 * level->prolong_tmp_stride + idx0;
                 level->solver->u[idx_dst] -= level->prolong_tmp[idx_src];
             }
+        level->time_correct += gettime() - start;
 
         /* re-init the current-level solver (re-calc the residual) */
+        start = gettime();
         ret = mg2di_ell_relax_init(level->solver);
         if (ret < 0)
             return ret;
+        level->time_reinit += gettime() - start;
 
         /* post-correct relaxation */
+        start = gettime();
         for (int j = 0; j < ctx->nb_relax_pre; j++) {
             ret = mg2di_ell_relax_step(level->solver);
             if (ret < 0)
                 return ret;
         }
+        level->time_relax += gettime() - start;
+
+        level->count_cycles++;
     }
 
 finish:
@@ -304,9 +336,12 @@ static int mg_levels_init(MG2DContext *ctx)
 int mg2d_solve(MG2DContext *ctx)
 {
     MG2DInternal *priv = ctx->priv;
+    int64_t time_start;
     double res_prev, res_cur;
     int ret;
 
+    time_start = gettime();
+
     ret = mg_levels_init(ctx);
     if (ret < 0)
         return ret;
@@ -327,6 +362,10 @@ int mg2d_solve(MG2DContext *ctx)
         if (res_cur < ctx->tol) {
             mg2di_log(&priv->logger, 0, "converged on iteration %d, residual %g\n",
                       i, res_cur);
+
+            priv->time_solve += gettime() - time_start;
+            priv->count_solve++;
+
             return 0;
         }
 
@@ -504,3 +543,38 @@ void mg2d_solver_free(MG2DContext **pctx)
     free(ctx);
     *pctx = NULL;
 }
+
+void mg2d_print_stats(MG2DContext *ctx, const char *prefix)
+{
+    MG2DInternal *priv = ctx->priv;
+    MG2DLevel *level = priv->root;
+
+    if (!prefix)
+        prefix = "";
+
+    mg2di_log(&priv->logger, 2, "%s%ld solves; %g s total time; %g ms avg per call\n",
+              prefix, priv->count_solve, priv->time_solve / 1e6, priv->time_solve / 1e3 / priv->count_solve);
+
+    while (level) {
+        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;
+
+        mg2di_log(&priv->logger, 2, "%s%2.2f%% level %d: %ld cycles %g s total time %g ms avg per call || "
+                  "%2.2f%% relax %2.2f%% prolong %2.2f%% restrict %2.2f%% correct %2.2f%% reinit || "
+                  "%2.2f%% residual %2.2f%% correct %2.2f%% boundaries\n",
+                  prefix, level_total * 100.0 / priv->time_solve, level->depth, level->count_cycles,
+                  level_total / 1e6, level_total / 1e3 / level->count_cycles,
+                  level->time_relax    * 100.0 / level_total,
+                  level->time_prolong  * 100.0 / level_total,
+                  level->time_restrict * 100.0 / level_total,
+                  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,
+                  level->solver->time_boundaries * 100.0 / relax_total
+                  );
+
+        level = level->child;
+    }
+}
diff --git a/mg2d.h b/mg2d.h
index 9a0a5a5..a0c0c70 100644
--- a/mg2d.h
+++ b/mg2d.h
@@ -241,4 +241,12 @@ int mg2d_solve(MG2DContext *ctx);
  */
 void mg2d_solver_free(MG2DContext **ctx);
 
+/**
+ * Print the solver timing statistics using the logging callback.
+ *
+ * @param prefix the string prepended to each output line, may be NULL to mean
+ *               none
+ */
+void mg2d_print_stats(MG2DContext *ctx, const char *prefix);
+
 #endif /* MG2D_H */