Add and use a new timer API.

5 files changed, 166 insertions, 175 deletions

diff --git a/common.h b/common.h
index 0a36ca9..03b33f6 100644
--- a/common.h
+++ b/common.h
@@ -24,16 +24,18 @@
 #define MIN(x, y) ((x) > (y) ? (y) : (x))
 #define ARRAY_ELEMS(arr) (sizeof(arr) / sizeof(*arr))
 
+#include <stdio.h>
 #include <stdlib.h>
 #include <stdint.h>
-#include <sys/time.h>
-static inline int64_t gettime(void)
-{
-    struct timeval tv;
-    gettimeofday(&tv, NULL);
-    return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec;
-}
 
 #define FD_STENCIL_MAX 4
 
+#define mg2di_assert(x)                                 \
+do {                                                    \
+    if (!(x)) {                                         \
+        fprintf(stderr, "Assertion " #x " failed\n");   \
+        abort();                                        \
+    }                                                   \
+} while (0)
+
 #endif /* MG2D_COMMON_H */
diff --git a/ell_grid_solve.c b/ell_grid_solve.c
index e792b47..a01f124 100644
--- a/ell_grid_solve.c
+++ b/ell_grid_solve.c
@@ -37,6 +37,7 @@
 #include "mg2d_constants.h"
 #include "ndarray.h"
 #include "residual_calc.h"
+#include "timer.h"
 
 static const double relax_factors[FD_STENCIL_MAX] = {
     [0] = 1.0 / 5,
@@ -118,9 +119,8 @@ static void residual_calc(EGSContext *ctx)
 {
     EGSInternal *priv = ctx->priv;
     const double *diff_coeffs[MG2D_DIFF_COEFF_NB];
-    int64_t start;
 
-    start = gettime();
+    mg2di_timer_start(&ctx->timer_res_calc);
 
     for (int i = 0; i < ARRAY_ELEMS(diff_coeffs); i++)
         diff_coeffs[i] = NDPTR2D(ctx->diff_coeffs[i], priv->residual_calc_offset[1], priv->residual_calc_offset[0]);
@@ -134,8 +134,7 @@ static void residual_calc(EGSContext *ctx)
                         ctx->rhs->stride[0],
                         diff_coeffs, ctx->diff_coeffs[0]->stride[0], priv->fd_factors);
 
-    ctx->time_res_calc += gettime() - start;
-    ctx->count_res++;
+    mg2di_timer_stop(&ctx->timer_res_calc);
 }
 
 static void boundaries_apply_fixval(double *dst, const ptrdiff_t dst_stride[2],
@@ -200,9 +199,9 @@ static void boundaries_apply(EGSContext *ctx, int init)
     static const enum MG2DBCType bnd_type_order[] = { MG2D_BC_TYPE_FIXVAL, MG2D_BC_TYPE_REFLECT, MG2D_BC_TYPE_FALLOFF };
     EGSInternal *priv = ctx->priv;
     const ptrdiff_t strides[2] = { 1, ctx->u->stride[0] };
-    int64_t start, start_bnd;
 
-    start = gettime();
+    mg2di_timer_start(&ctx->timer_bnd);
+
     for (int order_idx = 0; order_idx < ARRAY_ELEMS(bnd_type_order); order_idx++) {
         for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) {
             MG2DBoundary *bnd = ctx->boundaries[i];
@@ -222,31 +221,28 @@ static void boundaries_apply(EGSContext *ctx, int init)
             switch (bnd->type) {
             case MG2D_BC_TYPE_FIXVAL:
                 if (init && !priv->bnd_zero[i]) {
-                    start_bnd = gettime();
+                    mg2di_timer_start(&ctx->timer_bnd_fixval);
                     boundaries_apply_fixval(dst, dst_strides, bnd->val,
                                             bnd->val_stride, size_boundary);
-                    ctx->time_bnd_fixval += gettime() - start_bnd;
-                    ctx->count_bnd_fixval++;
+                    mg2di_timer_stop(&ctx->timer_bnd_fixval);
                 }
                 break;
             case MG2D_BC_TYPE_REFLECT:
-                start_bnd = gettime();
+                mg2di_timer_start(&ctx->timer_bnd_reflect);
                 boundaries_apply_reflect(dst, dst_strides, size_boundary);
-                ctx->time_bnd_reflect += gettime() - start_bnd;
-                ctx->count_bnd_reflect++;
+                mg2di_timer_stop(&ctx->timer_bnd_reflect);
                 break;
             case MG2D_BC_TYPE_FALLOFF:
-                start_bnd = gettime();
+                mg2di_timer_start(&ctx->timer_bnd_falloff);
                 boundaries_apply_falloff(dst, dst_strides, size_boundary, ctx);
-                ctx->time_bnd_falloff += gettime() - start_bnd;
-                ctx->count_bnd_falloff++;
+                mg2di_timer_stop(&ctx->timer_bnd_falloff);
                 break;
             }
         }
     }
 
     /* fill in the corner ghosts */
-    start_bnd = gettime();
+    mg2di_timer_start(&ctx->timer_bnd_corners);
     for (int pos_y = 0; pos_y < 2; pos_y++) {
         enum MG2DBoundaryLoc loc_y = mg2d_bnd_id(1, pos_y);
         MG2DBoundary        *bnd_y = ctx->boundaries[loc_y];
@@ -277,11 +273,8 @@ static void boundaries_apply(EGSContext *ctx, int init)
                 }
         }
     }
-    ctx->time_bnd_corners += gettime() - start_bnd;
-    ctx->count_bnd_corners++;
-
-    ctx->time_boundaries += gettime() - start;
-    ctx->count_boundaries++;
+    mg2di_timer_stop(&ctx->timer_bnd_corners);
+    mg2di_timer_stop(&ctx->timer_bnd);
 }
 
 #if HAVE_EXTERNAL_ASM
@@ -316,14 +309,11 @@ static int residual_add_task(void *arg, unsigned int job_idx, unsigned int threa
 static int solve_relax_step(EGSContext *ctx)
 {
     EGSRelaxContext *r = ctx->solver_data;
-    int64_t start;
 
-    start = gettime();
 
+    mg2di_timer_start(&r->timer_correct);
     tp_execute(ctx->tp, ctx->domain_size[1], residual_add_task, ctx);
-
-    r->time_correct += gettime() - start;
-    r->count_correct++;
+    mg2di_timer_stop(&r->timer_correct);
 
     boundaries_apply(ctx, 0);
     residual_calc(ctx);
@@ -583,10 +573,9 @@ static int solve_exact(EGSContext *ctx)
     EGSInternal *priv = ctx->priv;
     EGSExactContext  *ec = ctx->solver_data;
     EGSExactInternal *e = &priv->e;
-    int64_t start;
     int ret;
 
-    start = gettime();
+    mg2di_timer_start(&ec->timer_mat_construct);
 
     if (!e->mat_filled) {
         memset(e->mat, 0, SQR(e->N) * sizeof(*e->mat));
@@ -602,29 +591,27 @@ static int solve_exact(EGSContext *ctx)
            e->rhs[mat_row_idx] = e->rhs_factor[mat_row_idx] * (*NDPTR2D(ctx->rhs, idx0, idx1)) + e->rhs_mat[mat_row_idx];
        }
 
-    ec->time_mat_construct += gettime() - start;
-    ec->count_mat_construct++;
-
-    start = gettime();
+    mg2di_timer_stop(&ec->timer_mat_construct);
+    mg2di_timer_start(&ec->timer_bicgstab);
 
     ret = mg2di_bicgstab_solve(e->bicgstab, e->mat, e->rhs, e->x);
     if (ret >= 0)
         ec->bicgstab_iterations += ret;
 
-    ec->time_bicgstab_solve += gettime() - start;
-    ec->count_bicgstab_solve++;
+    mg2di_timer_stop(&ec->timer_bicgstab);
 
     if (ret < 0) {
         char equed = 'N';
         double cond, ferr, berr, rpivot;
 
-        start = gettime();
+        mg2di_timer_start(&ec->timer_lu_solve);
 
         ret = LAPACKE_dgesvx(LAPACK_COL_MAJOR, 'N', 'N', e->N, 1,
                              e->mat, e->N, e->mat_f, e->N, e->ipiv, &equed, NULL, NULL,
                              e->rhs, e->N, e->x, e->N, &cond, &ferr, &berr, &rpivot);
         if (ret == 0)
             ret = LAPACKE_dgetri(LAPACK_COL_MAJOR, e->N, e->mat_f, e->N, e->ipiv);
+        mg2di_timer_stop(&ec->timer_lu_solve);
         if (ret != 0) {
             mg2di_log(&ctx->logger, MG2D_LOG_ERROR,
                       "Error solving the linear system: %d\n", ret);
@@ -635,21 +622,17 @@ static int solve_exact(EGSContext *ctx)
                   "condition number %16.16g; forward error %16.16g backward error %16.16g\n",
                   e->N, e->N, cond, ferr, berr);
 
-        ec->time_lu_solve += gettime() - start;
-        ec->count_lu_solve++;
-
         ret = mg2di_bicgstab_init(e->bicgstab, e->mat_f, e->x);
         if (ret < 0)
             return ret;
     }
 
-    start = gettime();
+    mg2di_timer_start(&ec->timer_export);
 
     for (size_t idx1 = 0; idx1 < ctx->domain_size[1]; idx1++)
         memcpy(ctx->u->data + idx1 * ctx->u->stride[0], e->x + idx1 * ctx->domain_size[0], ctx->domain_size[0] * sizeof(*e->x));
 
-    ec->time_export += gettime() - start;
-    ec->count_export++;
+    mg2di_timer_stop(&ec->timer_export);
 
     boundaries_apply(ctx, 0);
     residual_calc(ctx);
@@ -659,10 +642,9 @@ static int solve_exact(EGSContext *ctx)
 
 int mg2di_egs_solve(EGSContext *ctx)
 {
-    int64_t start;
     int ret;
 
-    start = gettime();
+    mg2di_timer_start(&ctx->timer_solve);
 
     switch (ctx->solver_type) {
     case EGS_SOLVER_RELAXATION: ret = solve_relax_step(ctx); break;
@@ -670,7 +652,7 @@ int mg2di_egs_solve(EGSContext *ctx)
     default:                    ret = -EINVAL;
     }
 
-    ctx->time_total += gettime() - start;
+    mg2di_timer_stop(&ctx->timer_solve);
 
     return ret;
 }
@@ -678,16 +660,15 @@ int mg2di_egs_solve(EGSContext *ctx)
 int mg2di_egs_init(EGSContext *ctx, int flags)
 {
     EGSInternal *priv = ctx->priv;
-    int64_t start, start_init;
-    int ret;
-
-    start = gettime();
+    int ret = 0;
 
-    start_init = gettime();
+    mg2di_timer_start(&ctx->timer_solve);
+    mg2di_timer_start(&ctx->timer_init);
 
     if (ctx->step[0] <= DBL_EPSILON || ctx->step[1] <= DBL_EPSILON) {
         mg2di_log(&ctx->logger, 0, "Spatial step size too small\n");
-        return -EINVAL;
+        ret = -EINVAL;
+        goto finish;
     }
 
     if (ctx->solver_type == EGS_SOLVER_RELAXATION) {
@@ -718,7 +699,7 @@ int mg2di_egs_init(EGSContext *ctx, int flags)
     if (!ctx->tp) {
         ret = tp_init(&priv->tp_internal, 1);
         if (ret < 0)
-            return ret;
+            goto finish;
         ctx->tp = priv->tp_internal;
     }
 
@@ -735,7 +716,8 @@ int mg2di_egs_init(EGSContext *ctx, int flags)
         default:
             mg2di_log(&ctx->logger, 0, "Invalid finite difference stencil: %zd\n",
                       ctx->fd_stencil);
-            return -EINVAL;
+            ret = -EINVAL;
+            goto finish;
         }
 
         if (e->scratch_lines_allocated < nb_threads) {
@@ -745,8 +727,10 @@ int mg2di_egs_init(EGSContext *ctx, int flags)
 
             e->scratch_lines = calloc(e->N_ghosts * nb_threads,
                                       sizeof(*e->scratch_lines));
-            if (!e->scratch_lines)
-                return -ENOMEM;
+            if (!e->scratch_lines) {
+                ret = -ENOMEM;
+                goto finish;
+            }
             e->scratch_lines_allocated = nb_threads;
         }
     }
@@ -782,17 +766,19 @@ int mg2di_egs_init(EGSContext *ctx, int flags)
 
     ret = mg2di_residual_calc_init(priv->rescalc);
     if (ret < 0)
-        return ret;
+        goto finish;
 
-    ctx->time_init += gettime() - start_init;
-    ctx->count_init++;
+finish:
+    mg2di_timer_stop(&ctx->timer_init);
 
-    boundaries_apply(ctx, 1);
-    residual_calc(ctx);
+    if (ret >= 0) {
+        boundaries_apply(ctx, 1);
+        residual_calc(ctx);
+    }
 
-    ctx->time_total += gettime() - start;
+    mg2di_timer_stop(&ctx->timer_solve);
 
-    return 0;
+    return ret;
 }
 
 static int arrays_alloc(EGSContext *ctx, const size_t domain_size[2])
@@ -902,16 +888,28 @@ EGSContext *mg2di_egs_alloc(enum EGSType type, size_t domain_size[2])
         goto fail;
 
     switch (type) {
-    case EGS_SOLVER_RELAXATION:
-        ctx->solver_data = calloc(1, sizeof(EGSRelaxContext));
-        if (!ctx->solver_data)
+    case EGS_SOLVER_RELAXATION: {
+        EGSRelaxContext *r;
+        r = calloc(1, sizeof(EGSRelaxContext));
+        if (!r)
             goto fail;
+        ctx->solver_data = r;
+        mg2di_timer_init(&r->timer_correct);
         break;
-    case EGS_SOLVER_EXACT:
-        ctx->solver_data = calloc(1, sizeof(EGSExactContext));
-        if (!ctx->solver_data)
+        }
+    case EGS_SOLVER_EXACT: {
+        EGSExactContext *e;
+        e = calloc(1, sizeof(EGSExactContext));
+        if (!e)
             goto fail;
+        ctx->solver_data = e;
+
+        mg2di_timer_init(&e->timer_mat_construct);
+        mg2di_timer_init(&e->timer_bicgstab);
+        mg2di_timer_init(&e->timer_lu_solve);
+        mg2di_timer_init(&e->timer_export);
         break;
+        }
     default: goto fail;
     }
     ctx->solver_type = type;
@@ -931,6 +929,15 @@ EGSContext *mg2di_egs_alloc(enum EGSType type, size_t domain_size[2])
     if (!ctx->priv->rescalc)
         goto fail;
 
+    mg2di_timer_init(&ctx->timer_bnd);
+    mg2di_timer_init(&ctx->timer_bnd_fixval);
+    mg2di_timer_init(&ctx->timer_bnd_falloff);
+    mg2di_timer_init(&ctx->timer_bnd_reflect);
+    mg2di_timer_init(&ctx->timer_bnd_corners);
+    mg2di_timer_init(&ctx->timer_res_calc);
+    mg2di_timer_init(&ctx->timer_init);
+    mg2di_timer_init(&ctx->timer_solve);
+
     return ctx;
 fail:
     mg2di_egs_free(&ctx);
diff --git a/ell_grid_solve.h b/ell_grid_solve.h
index d46edd9..8f984ce 100644
--- a/ell_grid_solve.h
+++ b/ell_grid_solve.h
@@ -44,6 +44,7 @@
 #include "mg2d_boundary.h"
 #include "mg2d_constants.h"
 #include "ndarray.h"
+#include "timer.h"
 
 enum EGSType {
     /**
@@ -84,20 +85,15 @@ typedef struct EGSRelaxContext {
      */
     double relax_multiplier;
 
-    int64_t count_correct;
-    int64_t time_correct;
+    Timer timer_correct;
 } EGSRelaxContext;
 
 typedef struct EGSExactContext {
-    int64_t count_mat_construct;
-    int64_t time_mat_construct;
-    int64_t count_bicgstab_solve;
+    Timer   timer_mat_construct;
+    Timer   timer_bicgstab;
     int64_t bicgstab_iterations;
-    int64_t time_bicgstab_solve;
-    int64_t count_lu_solve;
-    int64_t time_lu_solve;
-    int64_t count_export;
-    int64_t time_export;
+    Timer   timer_lu_solve;
+    Timer   timer_export;
 } EGSExactContext;
 
 typedef struct EGSContext {
@@ -195,21 +191,14 @@ typedef struct EGSContext {
     NDArray *diff_coeffs[MG2D_DIFF_COEFF_NB];
 
     /* timings */
-    int64_t time_boundaries;
-    int64_t count_boundaries;
-    int64_t time_bnd_fixval;
-    int64_t count_bnd_fixval;
-    int64_t time_bnd_falloff;
-    int64_t count_bnd_falloff;
-    int64_t time_bnd_reflect;
-    int64_t count_bnd_reflect;
-    int64_t time_bnd_corners;
-    int64_t count_bnd_corners;
-    int64_t time_res_calc;
-    int64_t count_res;
-    int64_t time_init;
-    int64_t count_init;
-    int64_t time_total;
+    Timer   timer_bnd;
+    Timer   timer_bnd_fixval;
+    Timer   timer_bnd_falloff;
+    Timer   timer_bnd_reflect;
+    Timer   timer_bnd_corners;
+    Timer   timer_res_calc;
+    Timer   timer_init;
+    Timer   timer_solve;
 } EGSContext;
 
 #define EGS_INIT_FLAG_SAME_DIFF_COEFFS (1 << 0)
diff --git a/meson.build b/meson.build
index 6d3b356..0c55f2e 100644
--- a/meson.build
+++ b/meson.build
@@ -12,6 +12,7 @@ lib_src = [
     'mg2d.c',
     'ndarray.c',
     'residual_calc.c',
+    'timer.c',
     'transfer.c',
 ]
 lib_obj = []
diff --git a/mg2d.c b/mg2d.c
index a22a0a5..6639879 100644
--- a/mg2d.c
+++ b/mg2d.c
@@ -34,6 +34,7 @@
 #include "mg2d_boundary_internal.h"
 #include "mg2d_constants.h"
 #include "ndarray.h"
+#include "timer.h"
 #include "transfer.h"
 
 typedef struct MG2DLevel {
@@ -52,11 +53,11 @@ typedef struct MG2DLevel {
 
     /* timings */
     int64_t count_cycles;
-    int64_t time_solve;
-    int64_t time_prolong;
-    int64_t time_restrict;
-    int64_t time_correct;
-    int64_t time_reinit;
+    Timer   timer_solve;
+    Timer   timer_prolong;
+    Timer   timer_restrict;
+    Timer   timer_correct;
+    Timer   timer_reinit;
 } MG2DLevel;
 
 struct MG2DInternal {
@@ -74,10 +75,8 @@ struct MG2DInternal {
 
     int cpuflags;
 
-    int64_t time_solve;
-    int64_t count_solve;
-    int64_t time_levels_init;
-    int64_t count_levels_init;
+    Timer  timer_solve;
+    Timer   timer_levels_init;
 };
 
 static void log_callback(MG2DLogger *log, int level, const char *fmt, va_list vl)
@@ -113,16 +112,13 @@ static int coarse_correct_task(void *arg, unsigned int job_idx, unsigned int thr
 static int mg_relax_step(MG2DContext *ctx, MG2DLevel *level, const char *step_desc)
 {
     double res_old;
-    int64_t start;
     int ret;
 
     res_old = level->solver->residual_max;
 
-    start = gettime();
-
+    mg2di_timer_start(&level->timer_solve);
     ret = mg2di_egs_solve(level->solver);
-
-    level->time_solve += gettime() - start;
+    mg2di_timer_stop(&level->timer_solve);
 
     if (ret < 0)
         return ret;
@@ -140,7 +136,7 @@ static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level)
     /* on the refined levels, use zero as the initial guess for the
      * solution (correction for the upper level) */
     if (level->depth > 0) {
-        int64_t start = gettime();
+        mg2di_timer_start(&level->timer_reinit);
 
         memset(level->solver->u->data, 0, sizeof(*level->solver->u->data) * level->solver->u->stride[0] *
                 level->solver->domain_size[1]);
@@ -149,9 +145,8 @@ static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level)
         ret = mg2di_egs_init(level->solver, level->egs_init_flags);
         if (ret < 0)
             return ret;
+        mg2di_timer_stop(&level->timer_reinit);
         level->egs_init_flags |= EGS_INIT_FLAG_SAME_DIFF_COEFFS;
-
-        level->time_reinit += gettime() - start;
     }
 
     res_old = level->solver->residual_max;
@@ -168,7 +163,6 @@ static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level)
 
     for (int i = 0; i < ctx->nb_cycles; i++) {
         double res_prev;
-        int64_t start;
 
         /* pre-restrict relaxation */
         for (int j = 0; j < ctx->nb_relax_pre; j++) {
@@ -178,12 +172,12 @@ static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level)
         }
 
         /* restrict the residual as to the coarser-level rhs */
-        start = gettime();
+        mg2di_timer_start(&level->timer_restrict);
         ret = mg2di_gt_transfer(level->transfer_restrict, level->child->solver->rhs,
                                 level->solver->residual);
+        mg2di_timer_stop(&level->timer_restrict);
         if (ret < 0)
             return ret;
-        level->time_restrict += gettime() - start;
 
         /* solve on the coarser level */
         ret = mg_solve_subgrid(ctx, level->child);
@@ -191,27 +185,25 @@ static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level)
             return ret;
 
         /* prolongate the coarser-level correction */
-        start = gettime();
+        mg2di_timer_start(&level->timer_prolong);
         ret = mg2di_gt_transfer(level->transfer_prolong, level->prolong_tmp,
                                 level->child->solver->u);
+        mg2di_timer_stop(&level->timer_prolong);
         if (ret < 0)
             return ret;
-        level->time_prolong += gettime() - start;
 
         /* apply the correction */
-        start = gettime();
-
+        mg2di_timer_start(&level->timer_correct);
         tp_execute(ctx->priv->tp, level->solver->domain_size[1], coarse_correct_task, level);
-
-        level->time_correct += gettime() - start;
+        mg2di_timer_stop(&level->timer_correct);
 
         /* re-init the current-level solver (re-calc the residual) */
         res_prev = level->solver->residual_max;
-        start = gettime();
+        mg2di_timer_start(&level->timer_reinit);
         ret = mg2di_egs_init(level->solver, 0);
+        mg2di_timer_stop(&level->timer_reinit);
         if (ret < 0)
             return ret;
-        level->time_reinit += gettime() - start;
 
         log_relax_step(ctx, level, "correct", res_prev, level->solver->residual_max);
 
@@ -488,7 +480,6 @@ int mg2d_solve(MG2DContext *ctx)
     MG2DInternal *priv = ctx->priv;
     MG2DLevel    *root;
     EGSContext *s_root;
-    int64_t time_start, start;
     double res_orig, res_prev, res_cur;
     int ret;
 
@@ -507,24 +498,19 @@ int mg2d_solve(MG2DContext *ctx)
     root   = priv->root;
     s_root = root->solver;
 
-    time_start = gettime();
-
-    start = gettime();
+    mg2di_timer_start(&priv->timer_solve);
 
+    mg2di_timer_start(&priv->timer_levels_init);
     ret = mg_levels_init(ctx);
+    mg2di_timer_stop(&priv->timer_levels_init);
     if (ret < 0)
-        return ret;
-
-    priv->time_levels_init += gettime() - start;
-    priv->count_levels_init++;
-
-    start = gettime();
+        goto finish;
 
+    mg2di_timer_start(&root->timer_reinit);
     ret = mg2di_egs_init(s_root, 0);
+    mg2di_timer_stop(&root->timer_reinit);
     if (ret < 0)
-        return ret;
-
-    root->time_reinit += gettime() - start;
+        goto finish;
 
     res_orig = s_root->residual_max;
     res_prev = res_orig;
@@ -542,10 +528,7 @@ int mg2d_solve(MG2DContext *ctx)
 
             //mg2di_ndarray_copy(priv->u, s_root->u);
 
-            priv->time_solve += gettime() - time_start;
-            priv->count_solve++;
-
-            return 0;
+            goto finish;
         }
 
         mg2di_log(&priv->logger, MG2D_LOG_VERBOSE,
@@ -564,6 +547,8 @@ int mg2d_solve(MG2DContext *ctx)
               "Maximum number of iterations (%d) reached\n", ctx->maxiter);
 fail:
     mg2di_log(&priv->logger, MG2D_LOG_ERROR, "The solver failed to converge\n");
+finish:
+    mg2di_timer_stop(&priv->timer_solve);
     return ret;
 }
 
@@ -607,6 +592,12 @@ static MG2DLevel *mg_level_alloc(enum EGSType type, const size_t domain_size)
     if (!level->solver)
         goto fail;
 
+    mg2di_timer_init(&level->timer_solve);
+    mg2di_timer_init(&level->timer_prolong);
+    mg2di_timer_init(&level->timer_restrict);
+    mg2di_timer_init(&level->timer_correct);
+    mg2di_timer_init(&level->timer_reinit);
+
     return level;
 fail:
     mg_level_free(&level);
@@ -738,6 +729,9 @@ MG2DContext *mg2d_solver_alloc(size_t domain_size)
     ctx->log_level     = MG2D_LOG_INFO;
     ctx->nb_threads    = 1;
 
+    mg2di_timer_init(&priv->timer_levels_init);
+    mg2di_timer_init(&priv->timer_solve);
+
     return ctx;
 fail:
     mg2d_solver_free(&ctx);
@@ -788,8 +782,8 @@ void mg2d_print_stats(MG2DContext *ctx, const char *prefix)
         prefix = "";
 
     mg2di_log(&priv->logger, MG2D_LOG_VERBOSE, "%s%ld solves; %g s total time; %g ms avg per call; %g avg cycles per solve\n",
-              prefix, priv->count_solve, priv->time_solve / 1e6, priv->time_solve / 1e3 / priv->count_solve,
-              (double)level->count_cycles / priv->count_solve);
+              prefix, priv->timer_solve.nb_runs, priv->timer_solve.time_nsec / 1e9, priv->timer_solve.time_nsec / 1e6 / priv->timer_solve.nb_runs,
+              (double)level->count_cycles / priv->timer_solve.nb_runs);
 
     while (level) {
         char buf[1024], *p;
@@ -798,70 +792,68 @@ void mg2d_print_stats(MG2DContext *ctx, const char *prefix)
         EGSRelaxContext *r = NULL;
         EGSExactContext *e = NULL;
 
-        int64_t level_total = level->time_solve + level->time_prolong + level->time_restrict +
-                              level->time_correct + level->time_reinit;
+        int64_t level_total = level->timer_solve.time_nsec + level->timer_prolong.time_nsec + level->timer_restrict.time_nsec +
+                              level->timer_correct.time_nsec + level->timer_reinit.time_nsec;
 
         if (level->solver->solver_type == EGS_SOLVER_RELAXATION)
             r = level->solver->solver_data;
         else if (level->solver->solver_type == EGS_SOLVER_EXACT)
             e = level->solver->solver_data;
 
-        level_total = level->time_solve + level->time_prolong + level->time_restrict +
-                              level->time_correct + level->time_reinit;
-
         levels_total += level_total;
 
         p = buf;
 
         ret = snprintf(p, sizeof(buf) - (p - buf),
                        "%2.2f%% level %d: %ld cycles %g s total time %g ms avg per call",
-                       level_total * 100.0 / priv->time_solve, level->depth, level->count_cycles,
-                       level_total / 1e6, level_total / 1e3 / level->count_cycles);
+                       level_total * 100.0 / priv->timer_solve.time_nsec, level->depth, level->count_cycles,
+                       level_total / 1e9, level_total / 1e6 / level->count_cycles);
         if (ret > 0)
             p += ret;
 
         ret = snprintf(p, sizeof(buf) - (p - buf),
                        "||%2.2f%% solve %2.2f%% reinit ",
-                       level->time_solve  * 100.0 / level_total,
-                       level->time_reinit * 100.0 / level_total);
+                       level->timer_solve.time_nsec  * 100.0 / level_total,
+                       level->timer_reinit.time_nsec * 100.0 / level_total);
         if (ret > 0)
             p += ret;
 
         if (level->child) {
             ret = snprintf(p, sizeof(buf) - (p - buf),
                            "%2.2f%% prolong %2.2f%% restrict %2.2f%% correct",
-                           level->time_prolong  * 100.0 / level_total,
-                           level->time_restrict * 100.0 / level_total,
-                           level->time_correct  * 100.0 / level_total);
+                           level->timer_prolong.time_nsec  * 100.0 / level_total,
+                           level->timer_restrict.time_nsec * 100.0 / level_total,
+                           level->timer_correct.time_nsec  * 100.0 / level_total);
             if (ret > 0)
                 p += ret;
         }
 
         ret = snprintf(p, sizeof(buf) - (p - buf),
                        "||%2.2f%% init %2.2f%% residual %2.2f%% boundaries (%2.2f%% fixval %2.2f%% reflect %2.2f%% falloff %2.2f%% corners)",
-                       level->solver->time_init       * 100.0 / level->solver->time_total,
-                       level->solver->time_res_calc   * 100.0 / level->solver->time_total,
-                       level->solver->time_boundaries * 100.0 / level->solver->time_total,
-                       level->solver->time_bnd_fixval * 100.0 / level->solver->time_total,
-                       level->solver->time_bnd_reflect * 100.0 / level->solver->time_total,
-                       level->solver->time_bnd_falloff * 100.0 / level->solver->time_total,
-                       level->solver->time_bnd_corners * 100.0 / level->solver->time_total);
+                       level->solver->timer_init.time_nsec        * 100.0 / level->solver->timer_solve.time_nsec,
+                       level->solver->timer_res_calc.time_nsec    * 100.0 / level->solver->timer_solve.time_nsec,
+                       level->solver->timer_bnd.time_nsec         * 100.0 / level->solver->timer_solve.time_nsec,
+                       level->solver->timer_bnd_fixval.time_nsec  * 100.0 / level->solver->timer_solve.time_nsec,
+                       level->solver->timer_bnd_reflect.time_nsec * 100.0 / level->solver->timer_solve.time_nsec,
+                       level->solver->timer_bnd_falloff.time_nsec * 100.0 / level->solver->timer_solve.time_nsec,
+                       level->solver->timer_bnd_corners.time_nsec * 100.0 / level->solver->timer_solve.time_nsec);
         if (ret > 0)
             p += ret;
 
         if (r) {
             ret = snprintf(p, sizeof(buf) - (p - buf),
                            " %2.2f%% correct",
-                           r->time_correct * 100.0 / level->solver->time_total);
+                           r->timer_correct.time_nsec * 100.0 / level->solver->timer_solve.time_nsec);
             if (ret > 0)
                 p += ret;
         } else if (e) {
             ret = snprintf(p, sizeof(buf) - (p - buf),
                            " %2.2f%% const %2.2f%% bicgstab (%ld; %g it/slv) %2.2f%% lu (%ld) %2.2f%% export",
-                           e->time_mat_construct * 100.0 / level->solver->time_total,
-                           e->time_bicgstab_solve     * 100.0 / level->solver->time_total, e->count_bicgstab_solve, (double)e->bicgstab_iterations / e->count_bicgstab_solve,
-                           e->time_lu_solve     * 100.0 / level->solver->time_total, e->count_lu_solve,
-                           e->time_export     * 100.0 / level->solver->time_total);
+                           e->timer_mat_construct.time_nsec   * 100.0 / level->solver->timer_solve.time_nsec,
+                           e->timer_bicgstab.time_nsec        * 100.0 / level->solver->timer_solve.time_nsec,
+                           e->timer_bicgstab.nb_runs, (double)e->bicgstab_iterations / e->timer_bicgstab.nb_runs,
+                           e->timer_lu_solve.time_nsec        * 100.0 / level->solver->timer_solve.time_nsec, e->timer_lu_solve.nb_runs,
+                           e->timer_export.time_nsec          * 100.0 / level->solver->timer_solve.time_nsec);
             if (ret > 0)
                 p += ret;
         }
@@ -872,14 +864,14 @@ void mg2d_print_stats(MG2DContext *ctx, const char *prefix)
 
     mg2di_log(&priv->logger, MG2D_LOG_VERBOSE,
               "%s%2.2f%% levels init %g s total time %g ms avg per call\n",
-              prefix, priv->time_levels_init * 100.0 / priv->time_solve,
-              priv->time_levels_init / 1e6, priv->time_levels_init / 1e3 / priv->count_levels_init);
+              prefix, priv->timer_levels_init.time_nsec * 100.0 / priv->timer_solve.time_nsec,
+              priv->timer_levels_init.time_nsec / 1e9, priv->timer_levels_init.time_nsec / 1e6 / priv->timer_levels_init.nb_runs);
 
-    other = priv->time_solve - levels_total - priv->time_levels_init;
+    other = priv->timer_solve.time_nsec - levels_total - priv->timer_levels_init.time_nsec;
     mg2di_log(&priv->logger, MG2D_LOG_VERBOSE,
               "%s%2.2f%% other %g s total time %g ms avg per call\n",
-              prefix, other * 100.0 / priv->time_solve,
-              other / 1e6, other / 1e3 / priv->count_solve);
+              prefix, other * 100.0 / priv->timer_solve.time_nsec,
+              other / 1e9, other / 1e6 / priv->timer_solve.nb_runs);
 
 }