Switch to external ndarray library.

10 files changed, 101 insertions, 442 deletions

diff --git a/ell_grid_solve.c b/ell_grid_solve.c
index fd6e455..dad8861 100644
--- a/ell_grid_solve.c
+++ b/ell_grid_solve.c
@@ -27,6 +27,7 @@
 
 #include <lapacke.h>
 #include <mpi.h>
+#include <ndarray.h>
 
 #include "bicgstab.h"
 #include "common.h"
@@ -37,7 +38,6 @@
 #include "mg2d_boundary.h"
 #include "mg2d_boundary_internal.h"
 #include "mg2d_constants.h"
-#include "ndarray.h"
 #include "residual_calc.h"
 #include "timer.h"
 
@@ -169,7 +169,7 @@ static void residual_calc(EGSContext *ctx, int export_res)
     mg2di_timer_start(&ctx->timer_res_calc);
 
     for (int i = 0; i < ARRAY_ELEMS(diff_coeffs); i++)
-        diff_coeffs[i] = NDPTR2D(priv->diff_coeffs[i], offset[0], offset[1]);
+        diff_coeffs[i] = NDA_PTR2D(priv->diff_coeffs[i], offset[0], offset[1]);
 
     for (int bnd_idx = 0; bnd_idx < 4; bnd_idx++)
         if (ctx->boundaries[bnd_idx]->type == MG2D_BC_TYPE_REFLECT &&
@@ -177,9 +177,9 @@ static void residual_calc(EGSContext *ctx, int export_res)
             reflect_flags |= 1 << bnd_idx;
 
     mg2di_residual_calc(priv->rescalc, size, &ctx->residual_max,
-                        NDPTR2D(dst,      offset[0], offset[1]), dst->stride[0],
-                        NDPTR2D(ctx->u,   offset[0], offset[1]), ctx->u->stride[0],
-                        NDPTR2D(ctx->rhs, offset[0], offset[1]), ctx->rhs->stride[0],
+                        NDA_PTR2D(dst,      offset[0], offset[1]), dst->stride[0],
+                        NDA_PTR2D(ctx->u,   offset[0], offset[1]), ctx->u->stride[0],
+                        NDA_PTR2D(ctx->rhs, offset[0], offset[1]), ctx->rhs->stride[0],
                         diff_coeffs, priv->diff_coeffs[0]->stride[0],
                         export_res ? 0.0 : 1.0, export_res ? 1.0 : priv->r.relax_factor,
                         reflect_flags, FD_STENCIL_MAX);
@@ -580,7 +580,7 @@ static void mat_fill_row(EGSContext *ctx, double *scratch_line,
     EGSInternal *priv = ctx->priv;
     EGSExactInternal *e = &priv->e;
 
-    const ptrdiff_t idx_src_dc  = NDIDX2D(priv->diff_coeffs[0], idx0, idx1);
+    const ptrdiff_t idx_src_dc  = NDA_IDX2D(priv->diff_coeffs[0], idx0, idx1);
     const ptrdiff_t mat_row_idx = idx1 * ctx->domain_size[0] + idx0;
 
     int is_bnd[4], boundary;
@@ -785,7 +785,7 @@ static int solve_exact(EGSContext *ctx)
     for (ptrdiff_t idx1 = 0; idx1 < ctx->domain_size[1]; idx1++)
        for (ptrdiff_t idx0 = 0; idx0 < ctx->domain_size[0]; idx0++) {
            const ptrdiff_t mat_row_idx = idx1 * ctx->domain_size[0] + idx0;
-           e->rhs[mat_row_idx] = e->rhs_factor[mat_row_idx] * (*NDPTR2D(ctx->rhs, idx0, idx1)) + e->rhs_mat[mat_row_idx];
+           e->rhs[mat_row_idx] = e->rhs_factor[mat_row_idx] * (*NDA_PTR2D(ctx->rhs, idx0, idx1)) + e->rhs_mat[mat_row_idx];
        }
 
     mg2di_timer_stop(&ec->timer_mat_construct);
@@ -865,7 +865,7 @@ static int init_diff_coeffs_task(void *arg, unsigned int job_idx, unsigned int t
     const double factor = a->fact;
 
     for (ptrdiff_t idx1 = 0; idx1 < dst->shape[1]; idx1++)
-        *NDPTR2D(dst, idx1, job_idx) = *NDPTR2D(src, idx1, job_idx) * factor;
+        *NDA_PTR2D(dst, idx1, job_idx) = *NDA_PTR2D(src, idx1, job_idx) * factor;
 
     return 0;
 }
@@ -1126,8 +1126,8 @@ static int arrays_alloc(EGSContext *ctx, const size_t domain_size[2])
     EGSInternal *priv = ctx->priv;
     const DomainComponent *dc = &priv->dg->components[priv->local_component];
 
-    Slice slice_interior[2];
-    Slice slice_exterior[2];
+    NDASlice slice_interior[2];
+    NDASlice slice_exterior[2];
 
     size_t size_alloc[2];
     size_t size_dc[2];
@@ -1155,73 +1155,73 @@ static int arrays_alloc(EGSContext *ctx, const size_t domain_size[2])
     slice_interior[1].end   = slice_interior[1].start + dc->interior.size[0];
     slice_interior[1].step  = 1;
 
-    ret = mg2di_ndarray_alloc(&priv->u_base,      2, size_alloc, NDARRAY_ALLOC_ZERO);
+    ret = ndarray_alloc(&priv->u_base,      2, size_alloc, NDARRAY_ALLOC_ZERO);
     if (ret < 0)
         return ret;
-    ret = mg2di_ndarray_alloc(&priv->u_next_base, 2, size_alloc, NDARRAY_ALLOC_ZERO);
+    ret = ndarray_alloc(&priv->u_next_base, 2, size_alloc, NDARRAY_ALLOC_ZERO);
     if (ret < 0)
         return ret;
 
-    ret = mg2di_ndarray_slice(&ctx->u,       priv->u_base,      slice_interior);
+    ret = ndarray_slice(&ctx->u,       priv->u_base,      slice_interior);
     if (ret < 0)
         return ret;
-    ret = mg2di_ndarray_slice(&priv->u_next, priv->u_next_base, slice_interior);
+    ret = ndarray_slice(&priv->u_next, priv->u_next_base, slice_interior);
     if (ret < 0)
         return ret;
 
-    ret = mg2di_ndarray_slice(&ctx->u_exterior,       priv->u_base,      slice_exterior);
+    ret = ndarray_slice(&ctx->u_exterior,       priv->u_base,      slice_exterior);
     if (ret < 0)
         return ret;
-    ret = mg2di_ndarray_slice(&priv->u_next_exterior, priv->u_next_base, slice_exterior);
+    ret = ndarray_slice(&priv->u_next_exterior, priv->u_next_base, slice_exterior);
     if (ret < 0)
         return ret;
 
-    ret = mg2di_ndarray_alloc(&priv->rhs_base, 2, size_alloc, 0);
+    ret = ndarray_alloc(&priv->rhs_base, 2, size_alloc, 0);
     if (ret < 0)
         return ret;
 
-    ret = mg2di_ndarray_slice(&ctx->rhs, priv->rhs_base, slice_interior);
+    ret = ndarray_slice(&ctx->rhs, priv->rhs_base, slice_interior);
     if (ret < 0)
         return ret;
 
-    ret = mg2di_ndarray_alloc(&priv->residual_base, 2, size_alloc,
+    ret = ndarray_alloc(&priv->residual_base, 2, size_alloc,
                               NDARRAY_ALLOC_ZERO);
     if (ret < 0)
         return ret;
 
-    ret = mg2di_ndarray_slice(&ctx->residual, priv->residual_base, slice_interior);
+    ret = ndarray_slice(&ctx->residual, priv->residual_base, slice_interior);
     if (ret < 0)
         return ret;
 
     size_dc[0]    = size_alloc[0] * MG2D_DIFF_COEFF_NB;
     size_dc[1]    = size_alloc[1];
 
-    ret = mg2di_ndarray_alloc(&priv->diff_coeffs_public_base, 2,
+    ret = ndarray_alloc(&priv->diff_coeffs_public_base, 2,
                               size_dc, NDARRAY_ALLOC_ZERO);
     if (ret < 0)
         return ret;
-    ret = mg2di_ndarray_alloc(&priv->diff_coeffs_base, 2,
+    ret = ndarray_alloc(&priv->diff_coeffs_base, 2,
                               size_dc, NDARRAY_ALLOC_ZERO);
     if (ret < 0)
         return ret;
 
     for (int i = 0; i < ARRAY_ELEMS(ctx->diff_coeffs); i++) {
-        const Slice slice_dc[2] = { SLICE(i * size_alloc[0], (i + 1) * size_alloc[0], 1), SLICE_NULL };
+        const NDASlice slice_dc[2] = { NDASLICE(i * size_alloc[0], (i + 1) * size_alloc[0], 1), NDASLICE_NULL };
         NDArray *tmp;
 
-        ret = mg2di_ndarray_slice(&tmp, priv->diff_coeffs_public_base, slice_dc);
+        ret = ndarray_slice(&tmp, priv->diff_coeffs_public_base, slice_dc);
         if (ret < 0)
             return ret;
-        ret = mg2di_ndarray_slice(&ctx->diff_coeffs[i], tmp, slice_interior);
-        mg2di_ndarray_free(&tmp);
+        ret = ndarray_slice(&ctx->diff_coeffs[i], tmp, slice_interior);
+        ndarray_free(&tmp);
         if (ret < 0)
             return ret;
 
-        ret = mg2di_ndarray_slice(&tmp, priv->diff_coeffs_base, slice_dc);
+        ret = ndarray_slice(&tmp, priv->diff_coeffs_base, slice_dc);
         if (ret < 0)
             return ret;
-        ret = mg2di_ndarray_slice(&priv->diff_coeffs[i], tmp, slice_interior);
-        mg2di_ndarray_free(&tmp);
+        ret = ndarray_slice(&priv->diff_coeffs[i], tmp, slice_interior);
+        ndarray_free(&tmp);
         if (ret < 0)
             return ret;
     }
@@ -1384,28 +1384,28 @@ void mg2di_egs_free(EGSContext **pctx)
 
     exact_arrays_free(ctx);
 
-    mg2di_ndarray_free(&ctx->u);
-    mg2di_ndarray_free(&ctx->u_exterior);
-    mg2di_ndarray_free(&ctx->priv->u_base);
-    mg2di_ndarray_free(&ctx->priv->u_next);
-    mg2di_ndarray_free(&ctx->priv->u_next_exterior);
-    mg2di_ndarray_free(&ctx->priv->u_next_base);
+    ndarray_free(&ctx->u);
+    ndarray_free(&ctx->u_exterior);
+    ndarray_free(&ctx->priv->u_base);
+    ndarray_free(&ctx->priv->u_next);
+    ndarray_free(&ctx->priv->u_next_exterior);
+    ndarray_free(&ctx->priv->u_next_base);
 
-    mg2di_ndarray_free(&ctx->rhs);
-    mg2di_ndarray_free(&ctx->priv->rhs_base);
+    ndarray_free(&ctx->rhs);
+    ndarray_free(&ctx->priv->rhs_base);
 
-    mg2di_ndarray_free(&ctx->residual);
-    mg2di_ndarray_free(&ctx->priv->residual_base);
+    ndarray_free(&ctx->residual);
+    ndarray_free(&ctx->priv->residual_base);
 
     for (int i = 0; i < ARRAY_ELEMS(ctx->diff_coeffs); i++) {
-        mg2di_ndarray_free(&ctx->diff_coeffs[i]);
+        ndarray_free(&ctx->diff_coeffs[i]);
     }
-    mg2di_ndarray_free(&ctx->priv->diff_coeffs_public_base);
+    ndarray_free(&ctx->priv->diff_coeffs_public_base);
 
     for (int i = 0; i < ARRAY_ELEMS(ctx->priv->diff_coeffs); i++) {
-        mg2di_ndarray_free(&ctx->priv->diff_coeffs[i]);
+        ndarray_free(&ctx->priv->diff_coeffs[i]);
     }
-    mg2di_ndarray_free(&ctx->priv->diff_coeffs_base);
+    ndarray_free(&ctx->priv->diff_coeffs_base);
 
     for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++)
         mg2di_bc_free(&ctx->boundaries[i]);
diff --git a/ell_grid_solve.h b/ell_grid_solve.h
index 095c591..1e00393 100644
--- a/ell_grid_solve.h
+++ b/ell_grid_solve.h
@@ -41,12 +41,12 @@
 #include <threadpool.h>
 
 #include <mpi.h>
+#include <ndarray.h>
 
 #include "components.h"
 #include "log.h"
 #include "mg2d_boundary.h"
 #include "mg2d_constants.h"
-#include "ndarray.h"
 #include "timer.h"
 
 enum EGSType {
diff --git a/meson.build b/meson.build
index 4599633..a08656a 100644
--- a/meson.build
+++ b/meson.build
@@ -11,7 +11,6 @@ lib_src = [
     'ell_grid_solve.c',
     'log.c',
     'mg2d.c',
-    'ndarray.c',
     'residual_calc.c',
     'timer.c',
     'transfer.c',
@@ -25,12 +24,13 @@ cc = meson.get_compiler('c')
 libm = cc.find_library('m', required : false)
 libcblas   = cc.find_library('blas')
 liblapacke = cc.find_library('lapacke')
+libndarray = cc.find_library('ndarray')
 
 dep_tp  = declare_dependency(link_args : '-lthreadpool')
 
 dep_mpi = dependency('mpi', language: 'c')
 
-deps = [dep_tp, libm, libcblas, liblapacke, dep_mpi]
+deps = [dep_tp, libm, libcblas, liblapacke, libndarray, dep_mpi]
 
 cdata = configuration_data()
 cdata.set10('ARCH_X86',          false)
diff --git a/mg2d.c b/mg2d.c
index 3255ce7..08a91df 100644
--- a/mg2d.c
+++ b/mg2d.c
@@ -26,6 +26,7 @@
 #include <threadpool.h>
 
 #include <mpi.h>
+#include <ndarray.h>
 
 #include "common.h"
 #include "components.h"
@@ -36,7 +37,6 @@
 #include "mg2d_boundary.h"
 #include "mg2d_boundary_internal.h"
 #include "mg2d_constants.h"
-#include "ndarray.h"
 #include "timer.h"
 #include "transfer.h"
 
@@ -505,7 +505,7 @@ static int restrict_diff_coeffs(MG2DContext *ctx, MG2DLevel *l)
         NDArray *a_dc = l->depth > 0 ? l->solver->diff_coeffs[dc_idx] : priv->diff_coeffs_tmp[dc_idx];
 
         if (!l->depth)
-            mg2di_ndarray_copy(a_dc, priv->root->solver->diff_coeffs[dc_idx]);
+            ndarray_copy(a_dc, priv->root->solver->diff_coeffs[dc_idx]);
 
         for (int bnd_idx = 0; bnd_idx < ARRAY_ELEMS(dc->boundaries); bnd_idx++) {
             const int ci = mg2d_bnd_coord_idx(bnd_idx);
@@ -636,12 +636,12 @@ static int diff_coeffs_fixup(MG2DContext *ctx, MG2DLevel *l)
                 const ptrdiff_t stride_dst = a_dst->stride[ci];
                 const size_t      size_dst = l->solver->domain_size[!ci];
 
-                ret = mg2di_ndarray_slice(&dc_src, priv->dc_bnd_val[i][bnd_idx],
-                                          &SLICE(FD_STENCIL_MAX, -FD_STENCIL_MAX, 1));
+                ret = ndarray_slice(&dc_src, priv->dc_bnd_val[i][bnd_idx],
+                                          &NDASLICE(FD_STENCIL_MAX, -FD_STENCIL_MAX, 1));
                 if (ret < 0)
                     goto fail;
 
-                ret = mg2di_ndarray_alloc(&dc_dst, 1, &size_dst, 0);
+                ret = ndarray_alloc(&dc_dst, 1, &size_dst, 0);
                 if (ret < 0)
                     goto fail;
 
@@ -675,8 +675,8 @@ static int diff_coeffs_fixup(MG2DContext *ctx, MG2DLevel *l)
 
 fail:
                 mg2di_gt_free(&gt_bnd);
-                mg2di_ndarray_free(&dc_src);
-                mg2di_ndarray_free(&dc_dst);
+                ndarray_free(&dc_src);
+                ndarray_free(&dc_dst);
                 if (ret < 0)
                     return ret;
             }
@@ -788,26 +788,26 @@ static int mg_levels_init(MG2DContext *ctx)
     prev = NULL;
 
     if (priv->u) {
-        ret = mg2di_ndarray_copy(cur->solver->u,   priv->u);
+        ret = ndarray_copy(cur->solver->u,   priv->u);
         if (ret < 0)
             return ret;
-        mg2di_ndarray_free(&priv->u);
+        ndarray_free(&priv->u);
         ctx->u = cur->solver->u->data;
         ctx->u_stride = cur->solver->u->stride[0];
     }
 
     if (priv->rhs) {
-        ret = mg2di_ndarray_copy(cur->solver->rhs, priv->rhs);
+        ret = ndarray_copy(cur->solver->rhs, priv->rhs);
         if (ret < 0)
             return ret;
-        mg2di_ndarray_free(&priv->rhs);
+        ndarray_free(&priv->rhs);
         ctx->rhs = cur->solver->rhs->data;
         ctx->rhs_stride = cur->solver->rhs->stride[0];
     }
 
     if (ctx->diff_coeffs[0]->data == priv->diff_coeffs_tmp[0]->data) {
         for (int i = 0; i < ARRAY_ELEMS(ctx->diff_coeffs); i++) {
-            mg2di_ndarray_copy(cur->solver->diff_coeffs[i], priv->diff_coeffs_tmp[i]);
+            ndarray_copy(cur->solver->diff_coeffs[i], priv->diff_coeffs_tmp[i]);
             ctx->diff_coeffs[i]->data   = cur->solver->diff_coeffs[i]->data;
             ctx->diff_coeffs[i]->stride = cur->solver->diff_coeffs[i]->stride[0];
         }
@@ -1139,7 +1139,7 @@ int mg2d_solve(MG2DContext *ctx)
             mg2di_log(&priv->logger, MG2D_LOG_INFO, "converged on iteration %d, residual %g\n",
                       i, res_cur);
 
-            //mg2di_ndarray_copy(priv->u, s_root->u);
+            //ndarray_copy(priv->u, s_root->u);
 
             goto finish;
         }
@@ -1214,11 +1214,11 @@ static void mg_level_free(MG2DLevel **plevel)
     free(level->prolong_recvcounts);
     free(level->prolong_recvdispl);
 
-    mg2di_ndarray_free(&level->restrict_dst);
-    mg2di_ndarray_free(&level->prolong_dst);
+    ndarray_free(&level->restrict_dst);
+    ndarray_free(&level->prolong_dst);
 
-    mg2di_ndarray_free(&level->prolong_tmp);
-    mg2di_ndarray_free(&level->prolong_tmp_base);
+    ndarray_free(&level->prolong_tmp);
+    ndarray_free(&level->prolong_tmp_base);
     mg2di_egs_free(&level->solver);
 
     mg2di_gt_free(&level->transfer_restrict);
@@ -1247,12 +1247,12 @@ static MG2DLevel *mg_level_alloc(const DomainGeometry *dg,
 
     level->mpi_comm = comm;
 
-    ret = mg2di_ndarray_alloc(&level->prolong_tmp_base, 2,
+    ret = ndarray_alloc(&level->prolong_tmp_base, 2,
                               (size_t [2]){dc->interior.size[1] + 1, dc->interior.size[0] + 1}, 0);
     if (ret < 0)
         goto fail;
 
-    ret = mg2di_ndarray_slice(&level->prolong_tmp, level->prolong_tmp_base, (Slice [2]){ SLICE(0, -1, 1), SLICE(0, -1, 1) });
+    ret = ndarray_slice(&level->prolong_tmp, level->prolong_tmp_base, (NDASlice [2]){ NDASLICE(0, -1, 1), NDASLICE(0, -1, 1) });
     if (ret < 0)
         goto fail;
 
@@ -1359,12 +1359,12 @@ static int mg_interdomain_setup(MG2DContext *ctx, MG2DLevel *level)
     if (dg_fine->nb_components == 1)
         goto finish;
 
-    ret = mg2di_ndarray_alloc(&level->restrict_dst, 2,
+    ret = ndarray_alloc(&level->restrict_dst, 2,
                               (size_t [2]){ restrict_components[lc].size[1], restrict_components[lc].size[0] }, 0);
     if (ret < 0)
         goto finish;
 
-    ret = mg2di_ndarray_alloc(&level->prolong_dst, 2,
+    ret = ndarray_alloc(&level->prolong_dst, 2,
                               (size_t [2]){ prolong_components[lc].size[1], prolong_components[lc].size[0] }, 0);
     if (ret < 0)
         goto finish;
@@ -1600,14 +1600,14 @@ static MG2DContext *solver_alloc(DomainGeometry *dg, unsigned int local_componen
         }
     }
 
-    ret = mg2di_ndarray_alloc(&priv->u, 2, (size_t [2]){ domain_size[1], domain_size[0] },
+    ret = ndarray_alloc(&priv->u, 2, (size_t [2]){ domain_size[1], domain_size[0] },
                               NDARRAY_ALLOC_ZERO);
     if (ret < 0)
         goto fail;
     ctx->u        = priv->u->data;
     ctx->u_stride = priv->u->stride[0];
 
-    ret = mg2di_ndarray_alloc(&priv->rhs, 2, (size_t [2]){ domain_size[1], domain_size[0] },
+    ret = ndarray_alloc(&priv->rhs, 2, (size_t [2]){ domain_size[1], domain_size[0] },
                               NDARRAY_ALLOC_ZERO);
     if (ret < 0)
         goto fail;
@@ -1616,17 +1616,17 @@ static MG2DContext *solver_alloc(DomainGeometry *dg, unsigned int local_componen
 
     for (int i = 0; i < ARRAY_ELEMS(ctx->diff_coeffs); i++) {
         MG2DDiffCoeffs *dc;
-        const Slice slice[2] = { SLICE(FD_STENCIL_MAX, -FD_STENCIL_MAX, 1),
-                                 SLICE(FD_STENCIL_MAX, -FD_STENCIL_MAX, 1) };
+        const NDASlice slice[2] = { NDASLICE(FD_STENCIL_MAX, -FD_STENCIL_MAX, 1),
+                                    NDASLICE(FD_STENCIL_MAX, -FD_STENCIL_MAX, 1) };
 
-        ret = mg2di_ndarray_alloc(&priv->diff_coeffs_base[i], 2,
+        ret = ndarray_alloc(&priv->diff_coeffs_base[i], 2,
                                   (size_t [2]){ domain_size[1] + 2 * FD_STENCIL_MAX,
                                                 domain_size[0] + 2 * FD_STENCIL_MAX },
                                    NDARRAY_ALLOC_ZERO);
         if (ret < 0)
             goto fail;
 
-        ret = mg2di_ndarray_slice(&priv->diff_coeffs_tmp[i], priv->diff_coeffs_base[i], slice);
+        ret = ndarray_slice(&priv->diff_coeffs_tmp[i], priv->diff_coeffs_base[i], slice);
         if (ret < 0)
             goto fail;
 
@@ -1641,7 +1641,7 @@ static MG2DContext *solver_alloc(DomainGeometry *dg, unsigned int local_componen
         for (int bnd_idx = 0; bnd_idx < ARRAY_ELEMS(dc->boundaries); bnd_idx++) {
             const int ci = mg2d_bnd_coord_idx(bnd_idx);
 
-            ret = mg2di_ndarray_alloc(&priv->dc_bnd_val[i][bnd_idx], 1,
+            ret = ndarray_alloc(&priv->dc_bnd_val[i][bnd_idx], 1,
                                       (size_t [1]){ dg->domain_size[!ci] + 2 * FD_STENCIL_MAX },
                                       NDARRAY_ALLOC_ZERO);
             if (ret < 0)
@@ -1820,15 +1820,15 @@ void mg2d_solver_free(MG2DContext **pctx)
 
     mg2di_gt_free(&ctx->priv->transfer_init);
 
-    mg2di_ndarray_free(&ctx->priv->u);
-    mg2di_ndarray_free(&ctx->priv->rhs);
+    ndarray_free(&ctx->priv->u);
+    ndarray_free(&ctx->priv->rhs);
     for (int i = 0; i < ARRAY_ELEMS(ctx->priv->diff_coeffs_base); i++) {
         for (int j = 0; j < ARRAY_ELEMS(ctx->priv->dc_bnd_val[i]); j++)
-            mg2di_ndarray_free(&ctx->priv->dc_bnd_val[i][j]);
+            ndarray_free(&ctx->priv->dc_bnd_val[i][j]);
         free(ctx->diff_coeffs[i]);
 
-        mg2di_ndarray_free(&ctx->priv->diff_coeffs_tmp[i]);
-        mg2di_ndarray_free(&ctx->priv->diff_coeffs_base[i]);
+        ndarray_free(&ctx->priv->diff_coeffs_tmp[i]);
+        ndarray_free(&ctx->priv->diff_coeffs_base[i]);
     }
 
     mg2di_dg_free(&ctx->priv->dg);
@@ -2019,13 +2019,13 @@ int mg2d_init_guess(MG2DContext *ctx, const double *src,
             return ret;
     }
 
-    ret = mg2di_ndarray_wrap(&a_src, 2, (size_t [2]){ src_size[1], src_size[0] }, src,
+    ret = ndarray_wrap(&a_src, 2, (size_t [2]){ src_size[1], src_size[0] }, src,
                              (ptrdiff_t [2]){ src_stride, 1 });
     if (ret < 0)
         return ret;
 
     ret = mg2di_gt_transfer(priv->transfer_init, priv->u ? priv->u : priv->root->solver->u, a_src);
 
-    mg2di_ndarray_free(&a_src);
+    ndarray_free(&a_src);
     return ret;
 }
diff --git a/ndarray.c b/ndarray.c
deleted file mode 100644
index dee985c..0000000
--- a/ndarray.c
+++ /dev/null
@@ -1,272 +0,0 @@
-/*
- * N-dimensional strided array
- * Copyright 2018 Anton Khirnov <anton@khirnov.net>
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-#include <errno.h>
-#include <limits.h>
-#include <stddef.h>
-#include <stdlib.h>
-#include <string.h>
-
-#include "ndarray.h"
-
-struct NDArrayInternal {
-    size_t     alloc_size;
-    void      *data;
-
-    size_t    *shape;
-    ptrdiff_t *stride;
-};
-
-int mg2di_ndarray_alloc(NDArray **result, unsigned int dims,
-                        const size_t * const size, unsigned int flags)
-{
-    NDArray *a;
-    int ret;
-
-    if (!dims || dims >= INT_MAX)
-        return -EINVAL;
-
-    a = calloc(1, sizeof(*a));
-    if (!a)
-        return -ENOMEM;
-
-    a->priv = calloc(1, sizeof(*a->priv));
-    if (!a->priv) {
-        free(a);
-        return -ENOMEM;
-    }
-
-    a->priv->shape = calloc(dims, sizeof(*a->priv->shape));
-    if (!a->priv->shape) {
-        ret = -ENOMEM;
-        goto fail;
-    }
-    memcpy(a->priv->shape, size, dims * sizeof(*a->priv->shape));
-
-    a->priv->stride = calloc(dims, sizeof(*a->priv->stride));
-    if (!a->priv->stride) {
-        ret = -ENOMEM;
-        goto fail;
-    }
-
-    a->priv->stride[dims - 1] = 1;
-    for (int i = (int)dims - 2; i >= 0; i--)
-        a->priv->stride[i] = size[i + 1] * a->priv->stride[i + 1];
-
-    a->priv->alloc_size = a->priv->stride[0] * size[0] * sizeof(*a->data);
-    ret = posix_memalign((void**)&a->priv->data, 64, a->priv->alloc_size);
-    if (ret != 0) {
-        ret = -ret;
-        goto fail;
-    }
-
-    if (flags & NDARRAY_ALLOC_ZERO)
-        memset(a->priv->data, 0, a->priv->alloc_size);
-
-    a->dims   = dims;
-
-    a->data   = a->priv->data;
-    a->shape  = a->priv->shape;
-    a->stride = a->priv->stride;
-
-    *result = a;
-
-    return 0;
-fail:
-    mg2di_ndarray_free(&a);
-    return ret;
-}
-
-int mg2di_ndarray_slice(NDArray **result, NDArray *src,
-                        const Slice * const slice)
-{
-    NDArray *a;
-    ptrdiff_t offset;
-    int ret;
-
-    a = calloc(1, sizeof(*a));
-    if (!a)
-        return -ENOMEM;
-
-    a->priv = calloc(1, sizeof(*a->priv));
-    if (!a->priv) {
-        free(a);
-        return -ENOMEM;
-    }
-
-    a->priv->shape = calloc(src->dims, sizeof(*a->priv->shape));
-    if (!a->priv->shape) {
-        ret = -ENOMEM;
-        goto fail;
-    }
-    memcpy(a->priv->shape, src->priv->shape, src->dims * sizeof(*a->priv->shape));
-
-    a->priv->stride = calloc(src->dims, sizeof(*a->priv->stride));
-    if (!a->priv->stride) {
-        ret = -ENOMEM;
-        goto fail;
-    }
-    memcpy(a->priv->stride, src->priv->stride, src->dims * sizeof(*a->priv->shape));
-
-    offset = 0;
-    for (unsigned int i = 0; i < src->dims; i++) {
-        const Slice *s = &slice[i];
-        ptrdiff_t start, end, step;
-
-        if (s->start == PTRDIFF_MAX) {
-            start = 0;
-            end   = src->shape[i];
-            step  = 1;
-        } else {
-            start = s->start;
-            end   = s->end;
-            step  = s->step;
-
-            if (step < 1) {
-                ret = step ? -ENOSYS : -EINVAL;
-                goto fail;
-            }
-
-            if (start < 0)
-                start = src->shape[i] + start;
-            if (end <= 0)
-                end = src->shape[i] + end;
-
-            if (start >= end || end > src->shape[i]) {
-                ret = -EINVAL;
-                goto fail;
-            }
-        }
-
-        offset += src->stride[i] * start;
-
-        a->priv->stride[i] *= step;
-        a->priv->shape[i]   = (end - start) / step;
-    }
-
-    a->dims   = src->dims;
-
-    a->data   = src->data + offset;
-    a->shape  = a->priv->shape;
-    a->stride = a->priv->stride;
-
-    *result = a;
-    return 0;
-
-fail:
-    mg2di_ndarray_free(&a);
-    return ret;
-}
-
-int mg2di_ndarray_wrap(NDArray **result, unsigned int dims,
-                       const size_t * const size, double *data,
-                       const ptrdiff_t * const stride)
-{
-    NDArray *a;
-    int ret;
-
-    if (!dims || dims >= INT_MAX)
-        return -EINVAL;
-
-    a = calloc(1, sizeof(*a));
-    if (!a)
-        return -ENOMEM;
-
-    a->priv = calloc(1, sizeof(*a->priv));
-    if (!a->priv) {
-        free(a);
-        return -ENOMEM;
-    }
-
-    a->priv->shape = calloc(dims, sizeof(*a->priv->shape));
-    if (!a->priv->shape) {
-        ret = -ENOMEM;
-        goto fail;
-    }
-    memcpy(a->priv->shape, size, dims * sizeof(*a->priv->shape));
-
-    a->priv->stride = calloc(dims, sizeof(*a->priv->stride));
-    if (!a->priv->stride) {
-        ret = -ENOMEM;
-        goto fail;
-    }
-    memcpy(a->priv->stride, stride, dims * sizeof(*a->priv->shape));
-
-    a->dims   = dims;
-
-    a->data   = data;
-    a->shape  = a->priv->shape;
-    a->stride = a->priv->stride;
-
-    *result = a;
-
-    return 0;
-fail:
-    mg2di_ndarray_free(&a);
-    return ret;
-}
-
-void mg2di_ndarray_free(NDArray **pa)
-{
-    NDArray *a = *pa;
-
-    if (!a)
-        return;
-
-    if (a->priv) {
-        free(a->priv->data);
-        free(a->priv->shape);
-        free(a->priv->stride);
-    }
-    free(a->priv);
-
-    free(a);
-    *pa = NULL;
-}
-
-static int copy_axis(NDArray *dst, const NDArray *src, unsigned int axis,
-                     ptrdiff_t offset_dst, ptrdiff_t offset_src)
-{
-    if (dst->shape[axis] != src->shape[axis])
-        return -EINVAL;
-
-    if (axis == dst->dims - 1) {
-        if (dst->stride[axis] == 1 && src->stride[axis] == 1)
-            memcpy(dst->data + offset_dst, src->data + offset_src, sizeof(*dst->data) * dst->shape[axis]);
-        else {
-            for (size_t idx = 0; idx < dst->shape[axis]; idx++)
-                dst->data[offset_dst + idx * dst->stride[axis]] = src->data[offset_src + idx * src->stride[axis]];
-        }
-        return 0;
-    }
-
-    for (size_t idx = 0; idx < dst->shape[axis]; idx++)
-        copy_axis(dst, src, axis + 1, offset_dst + idx * dst->stride[axis], offset_src + idx * src->stride[axis]);
-
-    return 0;
-}
-
-int mg2di_ndarray_copy(NDArray *dst, const NDArray *src)
-{
-    const unsigned int dims = src->dims;
-
-    if (dims != dst->dims)
-        return -EINVAL;
-
-    return copy_axis(dst, src, 0, 0, 0);
-}
diff --git a/ndarray.h b/ndarray.h
deleted file mode 100644
index db58254..0000000
--- a/ndarray.h
+++ /dev/null
@@ -1,69 +0,0 @@
-/*
- * N-dimensional strided array
- * Copyright 2018 Anton Khirnov <anton@khirnov.net>
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-#ifndef MG2D_NDARRAY_H
-#define MG2D_NDARRAY_H
-
-#include <stddef.h>
-#include <stdint.h>
-
-typedef struct Slice {
-    ptrdiff_t start;
-    ptrdiff_t end;
-    ptrdiff_t step;
-} Slice;
-
-#define SLICE(start, end, step) ((Slice){ start, end, step })
-#define SLICE_NULL ((Slice){ PTRDIFF_MAX })
-
-typedef struct NDArrayInternal NDArrayInternal;
-
-typedef struct NDArray {
-    NDArrayInternal *priv;
-
-    unsigned int     dims;
-
-    double          *data;
-    const size_t    *shape;
-    const ptrdiff_t *stride;
-} NDArray;
-
-#define NDARRAY_ALLOC_ZERO (1 << 0)
-
-#define NDIDX1D(arr, x)       ((arr)->stride[0] * (x))
-#define NDIDX2D(arr, x, y)    ((arr)->stride[0] * (y) + (arr)->stride[1] * (x))
-#define NDIDX3D(arr, x, y, z) ((arr)->stride[0] * (z) + (arr)->stride[1] * (y) + (arr)->stride[0] * (x))
-
-#define NDPTR1D(arr, x)       ((arr)->data + NDIDX1D(arr, x))
-#define NDPTR2D(arr, x, y)    ((arr)->data + NDIDX2D(arr, x, y))
-#define NDPTR3D(arr, x, y, z) ((arr)->data + NDIDX2D(arr, x, y, z))
-
-int mg2di_ndarray_alloc(NDArray **result, unsigned int dims,
-                        const size_t * const size, unsigned int flags);
-void mg2di_ndarray_free(NDArray **a);
-
-int mg2di_ndarray_wrap(NDArray **result, unsigned int dims,
-                       const size_t * const size, double *data,
-                       const ptrdiff_t * const strides);
-
-int mg2di_ndarray_slice(NDArray **result, NDArray *src,
-                        const Slice * const slice);
-
-int mg2di_ndarray_copy(NDArray *dst, const NDArray *src);
-
-#endif // MG2D_ARRAY_H
diff --git a/relax_mpi_test.c b/relax_mpi_test.c
index 5e69b49..39dff1b 100644
--- a/relax_mpi_test.c
+++ b/relax_mpi_test.c
@@ -5,13 +5,13 @@
 #include <string.h>
 
 #include <mpi.h>
+#include <ndarray.h>
 
 #include "components.h"
 #include "ell_grid_solve.h"
 #include "log.h"
 #include "mg2d_boundary.h"
 #include "mg2d_constants.h"
-#include "ndarray.h"
 
 #define DOMAIN_SIZE 1.0
 #define FD_STENCIL 2
@@ -163,23 +163,23 @@ int main(int argc, char **argv)
     for (size_t y = 0; y < ctx->domain_size[1]; y++) {
         const double y_coord = (y + dc->interior.start[1]) * ctx->step[1];
 
-        memset(NDPTR2D(ctx->u, 0, y), 0, sizeof(*ctx->u->data) * ctx->domain_size[0]);
+        memset(NDA_PTR2D(ctx->u, 0, y), 0, sizeof(*ctx->u->data) * ctx->domain_size[0]);
 
         for (size_t x = 1; x < ctx->domain_size[0]; x++) {
             const double x_coord = x * ctx->step[0];
 
-            *NDPTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_02], x, y) =  1.0;
-            *NDPTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_20], x, y) =  1.0;
-            *NDPTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_11], x, y) =  1.0;
+            *NDA_PTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_02], x, y) =  1.0;
+            *NDA_PTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_20], x, y) =  1.0;
+            *NDA_PTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_11], x, y) =  1.0;
 
-            *NDPTR2D(ctx->rhs, x, y) = sol_dxx(x_coord, y_coord) + sol_dyy(x_coord, y_coord) + sol_dxy(x_coord, y_coord);
+            *NDA_PTR2D(ctx->rhs, x, y) = sol_dxx(x_coord, y_coord) + sol_dyy(x_coord, y_coord) + sol_dxy(x_coord, y_coord);
         }
 
-        memset(NDPTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_00], 0, y), 0,
+        memset(NDA_PTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_00], 0, y), 0,
                sizeof(*ctx->diff_coeffs[0]->data) * ctx->domain_size[0]);
-        memset(NDPTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_01], 0, y), 0,
+        memset(NDA_PTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_01], 0, y), 0,
                sizeof(*ctx->diff_coeffs[0]->data) * ctx->domain_size[0]);
-        memset(NDPTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_10], 0, y), 0,
+        memset(NDA_PTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_10], 0, y), 0,
                sizeof(*ctx->diff_coeffs[0]->data) * ctx->domain_size[0]);
     }
 
diff --git a/relax_test.c b/relax_test.c
index c2dd468..3d2fab2 100644
--- a/relax_test.c
+++ b/relax_test.c
@@ -3,12 +3,12 @@
 #include <stdio.h>
 #include <stdint.h>
 #include <string.h>
+#include <ndarray.h>
 
 #include "ell_grid_solve.h"
 #include "log.h"
 #include "mg2d_boundary.h"
 #include "mg2d_constants.h"
-#include "ndarray.h"
 
 #define ARRAY_ELEMS(x) (sizeof(x) / sizeof(*x))
 #define DOMAIN_SIZE 1.0
@@ -103,23 +103,23 @@ int main(int argc, char **argv)
     for (size_t y = 0; y < ctx->domain_size[1]; y++) {
         const double y_coord = y * ctx->step[1];
 
-        memset(NDPTR2D(ctx->u, 0, y), 0, sizeof(*ctx->u->data) * ctx->domain_size[0]);
+        memset(NDA_PTR2D(ctx->u, 0, y), 0, sizeof(*ctx->u->data) * ctx->domain_size[0]);
 
         for (size_t x = 0; x < ctx->domain_size[0]; x++) {
             const double x_coord = x * ctx->step[0];
 
-            *NDPTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_02], x, y) =  1.0;
-            *NDPTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_20], x, y) =  1.0;
-            *NDPTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_11], x, y) =  1.0;
+            *NDA_PTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_02], x, y) =  1.0;
+            *NDA_PTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_20], x, y) =  1.0;
+            *NDA_PTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_11], x, y) =  1.0;
 
-            *NDPTR2D(ctx->rhs, x, y) = sol_dxx(x_coord, y_coord) + sol_dyy(x_coord, y_coord) + sol_dxy(x_coord, y_coord);
+            *NDA_PTR2D(ctx->rhs, x, y) = sol_dxx(x_coord, y_coord) + sol_dyy(x_coord, y_coord) + sol_dxy(x_coord, y_coord);
         }
 
-        memset(NDPTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_00], 0, y), 0,
+        memset(NDA_PTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_00], 0, y), 0,
                sizeof(*ctx->diff_coeffs[0]->data) * ctx->domain_size[0]);
-        memset(NDPTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_01], 0, y), 0,
+        memset(NDA_PTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_01], 0, y), 0,
                sizeof(*ctx->diff_coeffs[0]->data) * ctx->domain_size[0]);
-        memset(NDPTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_10], 0, y), 0,
+        memset(NDA_PTR2D(ctx->diff_coeffs[MG2D_DIFF_COEFF_10], 0, y), 0,
                sizeof(*ctx->diff_coeffs[0]->data) * ctx->domain_size[0]);
     }
 
diff --git a/transfer.c b/transfer.c
index a791467..7fa41ee 100644
--- a/transfer.c
+++ b/transfer.c
@@ -21,6 +21,7 @@
 #include <errno.h>
 #include <limits.h>
 #include <math.h>
+#include <ndarray.h>
 #include <stddef.h>
 #include <stdint.h>
 #include <stdlib.h>
@@ -28,7 +29,6 @@
 
 #include "common.h"
 #include "cpu.h"
-#include "ndarray.h"
 #include "transfer.h"
 
 typedef struct GridTransfer {
diff --git a/transfer.h b/transfer.h
index 5c434e3..705de55 100644
--- a/transfer.h
+++ b/transfer.h
@@ -19,10 +19,10 @@
 #ifndef MG2D_TRANSFER_H
 #define MG2D_TRANSFER_H
 
+#include <ndarray.h>
 #include <stddef.h>
 #include <threadpool.h>
 
-#include "ndarray.h"
 
 enum GridTransferOperator {
     GRID_TRANSFER_LAGRANGE_1,