From ab72ad7bb19f08d78218d3558545f9f58e5b36e7 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Thu, 25 Jun 2020 20:51:02 +0200 Subject: Switch to external ndarray library. --- ell_grid_solve.c | 86 +++++++++--------- ell_grid_solve.h | 2 +- meson.build | 4 +- mg2d.c | 70 +++++++------- ndarray.c | 272 ------------------------------------------------------- ndarray.h | 69 -------------- relax_mpi_test.c | 18 ++-- relax_test.c | 18 ++-- transfer.c | 2 +- transfer.h | 2 +- 10 files changed, 101 insertions(+), 442 deletions(-) delete mode 100644 ndarray.c delete mode 100644 ndarray.h 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 #include +#include #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 #include +#include #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 #include +#include #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(>_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 - * - * 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 . - */ - -#include -#include -#include -#include -#include - -#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 - * - * 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 . - */ - -#ifndef MG2D_NDARRAY_H -#define MG2D_NDARRAY_H - -#include -#include - -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 #include +#include #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 #include #include +#include #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 #include #include +#include #include #include #include @@ -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 #include #include -#include "ndarray.h" enum GridTransferOperator { GRID_TRANSFER_LAGRANGE_1, -- cgit v1.2.3