Finish support for 4th order accuracy.

ABI and API break.

1 files changed, 32 insertions, 10 deletions

diff --git a/mg2d.c b/mg2d.c
index 9d98ec9..e4ba3b0 100644
--- a/mg2d.c
+++ b/mg2d.c
@@ -56,6 +56,8 @@ struct MG2DInternal {
 
     int64_t time_solve;
     int64_t count_solve;
+
+    double *boundaries_base[4];
 };
 
 static double findmax(double *arr, size_t len)
@@ -273,6 +275,23 @@ finish:
     return 0;
 }
 
+static void bnd_zero(EllRelaxBoundary *bdst, size_t nb_rows, size_t domain_size)
+{
+    for (size_t i = 0; i < nb_rows; i++) {
+        memset(bdst->val + i * bdst->val_stride - i, 0,
+               sizeof(*bdst->val) * (domain_size + 2 * i));
+    }
+}
+
+static void bnd_copy(EllRelaxBoundary *bdst, const double *src, ptrdiff_t src_stride,
+                     size_t nb_rows, size_t domain_size)
+{
+    for (size_t i = 0; i < nb_rows; i++) {
+        memcpy(bdst->val + i * bdst->val_stride - i, src + i * src_stride - i,
+               sizeof(*bdst->val) * (domain_size + 2 * i));
+    }
+}
+
 static int mg_levels_init(MG2DContext *ctx)
 {
     MG2DInternal *priv = ctx->priv;
@@ -301,13 +320,14 @@ static int mg_levels_init(MG2DContext *ctx)
             case MG2D_BC_TYPE_FIXVAL:
                 bdst->type = ELL_RELAX_BC_TYPE_FIXVAL;
                 if (!prev)
-                    memcpy(bdst->val, bsrc->val, sizeof(*bsrc->val) * bsrc->val_len);
+                    bnd_copy(bdst, bsrc->val, bsrc->val_stride, ctx->fd_stencil,
+                             cur->solver->domain_size[0]);
                 else
-                    memset(bdst->val, 0, sizeof(*bdst->val) * bdst->val_len);
+                    bnd_zero(bdst, ctx->fd_stencil, cur->solver->domain_size[0]);
                 break;
             case MG2D_BC_TYPE_FIXDIFF:
                 bdst->type = ELL_RELAX_BC_TYPE_FIXDIFF;
-                memset(bdst->val, 0, sizeof(*bdst->val) * bdst->val_len);
+                bnd_zero(bdst, ctx->fd_stencil, cur->solver->domain_size[0]);
                 break;
             default:
                 return -ENOSYS;
@@ -477,15 +497,18 @@ MG2DContext *mg2d_solver_alloc(size_t domain_size)
     priv->logger.opaque = ctx;
 
     for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) {
+        const size_t boundary_len_max = domain_size + 2 * (FD_STENCIL_MAX - 1);
+
         ctx->boundaries[i] = calloc(1, sizeof(*ctx->boundaries[i]));
         if (!ctx->boundaries[i])
             goto fail;
 
-        ret = posix_memalign((void**)&ctx->boundaries[i]->val, 32,
-                             sizeof(*ctx->boundaries[i]->val) * domain_size);
+        ret = posix_memalign((void**)&priv->boundaries_base[i], 32,
+                             sizeof(*ctx->boundaries[i]->val) * boundary_len_max * FD_STENCIL_MAX);
         if (ret != 0)
             goto fail;
-        ctx->boundaries[i]->val_len = domain_size;
+        ctx->boundaries[i]->val        = priv->boundaries_base[i] + FD_STENCIL_MAX - 1;
+        ctx->boundaries[i]->val_stride = boundary_len_max;
     }
 
     ret = mg_levels_alloc(ctx, domain_size);
@@ -532,14 +555,13 @@ void mg2d_solver_free(MG2DContext **pctx)
         ctx->priv->root = next;
     }
 
-    free(ctx->priv);
-
     for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) {
-        if (ctx->boundaries[i])
-            free(ctx->boundaries[i]->val);
+        free(ctx->priv->boundaries_base[i]);
         free(ctx->boundaries[i]);
     }
 
+    free(ctx->priv);
+
     free(ctx);
     *pctx = NULL;
 }