Finish support for 4th order accuracy.

ABI and API break.

6 files changed, 85 insertions, 49 deletions

diff --git a/common.h b/common.h
index d030ea4..7ecdf01 100644
--- a/common.h
+++ b/common.h
@@ -34,4 +34,6 @@ static inline int64_t gettime(void)
     return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec;
 }
 
+#define FD_STENCIL_MAX 2
+
 #endif /* MG2D_COMMON_H */
diff --git a/ell_relax.c b/ell_relax.c
index d79921a..e56296e 100644
--- a/ell_relax.c
+++ b/ell_relax.c
@@ -25,7 +25,7 @@
 #include "ell_relax.h"
 #include "log.h"
 
-#define FD_STENCIL_MAX 2
+#define CFL_FACT (1.0 / 7.0)
 
 struct EllRelaxInternal {
     ptrdiff_t stride;
@@ -34,6 +34,8 @@ struct EllRelaxInternal {
     double *residual_base;
     double *diff_coeffs_base[ELL_RELAX_DIFF_COEFF_NB];
 
+    double *boundary_val_base[4];
+
     void (*derivatives_calc)(EllRelaxContext *ctx, double *dst, ptrdiff_t idx, ptrdiff_t stride);
 };
 
@@ -141,35 +143,27 @@ static void residual_calc(EllRelaxContext *ctx)
     ctx->count_res++;
 }
 
-static void boundaries_apply_fixval(double *dst, ptrdiff_t dst_stride,
-                                    const double *src,
-                                    size_t boundary_size, ptrdiff_t boundary_stride,
-                                    int is_upper)
+static void boundaries_apply_fixval(double *dst, const ptrdiff_t dst_stride[2],
+                                    const double *src, ptrdiff_t src_stride,
+                                    size_t boundary_size)
 {
-    if (!is_upper)
-        boundary_stride *= -1;
-
-    for (size_t i = 0; i < boundary_size; i++) {
-        for (int j = 0; j <= FD_STENCIL_MAX; j++)
-            dst[j * boundary_stride] = *src;
-        dst += dst_stride;
-        src++;
+    for (int j = 0; j < FD_STENCIL_MAX; j++) {
+        for (ptrdiff_t i = -j; i < (ptrdiff_t)boundary_size + j; i++)
+            dst[i * dst_stride[0]] = src[i];
+        dst += dst_stride[1];
+        src += src_stride;
     }
 }
 
-static void boundaries_apply_fixdiff(double *dst, ptrdiff_t dst_stride,
-                                    const double *src,
-                                    size_t boundary_size, ptrdiff_t boundary_stride,
-                                    int is_upper)
+static void boundaries_apply_fixdiff(double *dst, const ptrdiff_t dst_stride[2],
+                                     const double *src, ptrdiff_t src_stride,
+                                     size_t boundary_size)
 {
-    if (!is_upper)
-        boundary_stride *= -1;
-
     for (size_t i = 0; i < boundary_size; i++) {
         for (int j = 1; j <= FD_STENCIL_MAX; j++)
-            dst[boundary_stride * j] = dst[-boundary_stride * j];
+            dst[dst_stride[1] * j] = dst[-dst_stride[1] * j];
 
-        dst += dst_stride;
+        dst += dst_stride[0];
     }
 }
 
@@ -188,15 +182,17 @@ static void boundaries_apply(EllRelaxContext *ctx)
         const size_t size_offset        = ctx->domain_size[!si];
 
         double *dst = ctx->u + boundary_def[i].is_upper * ((size_offset - 1) * stride_offset);
+        const ptrdiff_t dst_strides[] = { stride_boundary,
+            (boundary_def[i].is_upper ? 1 : -1) * stride_offset };
 
         switch (ctx->boundaries[i].type) {
         case ELL_RELAX_BC_TYPE_FIXVAL:
-            boundaries_apply_fixval(dst, stride_boundary, ctx->boundaries[i].val, size_boundary,
-                                    stride_offset, boundary_def[i].is_upper);
+            boundaries_apply_fixval(dst, dst_strides, ctx->boundaries[i].val,
+                                    ctx->boundaries[i].val_stride, size_boundary);
             break;
         case ELL_RELAX_BC_TYPE_FIXDIFF:
-            boundaries_apply_fixdiff(dst, stride_boundary, ctx->boundaries[i].val, size_boundary,
-                                     stride_offset, boundary_def[i].is_upper);
+            boundaries_apply_fixdiff(dst, dst_strides, ctx->boundaries[i].val,
+                                     ctx->boundaries[i].val_stride, size_boundary);
             break;
         }
     }
@@ -241,7 +237,7 @@ static void boundaries_apply(EllRelaxContext *ctx)
 int mg2di_ell_relax_step(EllRelaxContext *ctx)
 {
     EllRelaxInternal *priv = ctx->priv;
-    const double cfl_fac = (1.0 / 5.0) * ctx->step[0] * ctx->step[1];
+    const double cfl_fac = CFL_FACT * ctx->step[0] * ctx->step[1];
     int64_t start;
 
     start = gettime();
@@ -344,12 +340,16 @@ static int ell_relax_arrays_alloc(EllRelaxContext *ctx, const size_t domain_size
     priv->stride = stride;
 
     for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) {
-        const size_t boundary_size = domain_size[boundary_def[i].stride_idx];
-        ret = posix_memalign((void**)&ctx->boundaries[i].val, 32,
-                             sizeof(*ctx->boundaries[i].val) * boundary_size);
+        const size_t boundary_len0    = domain_size[boundary_def[i].stride_idx];
+        const size_t boundary_len_max = boundary_len0 + 2 * (FD_STENCIL_MAX - 1);
+
+        ret = posix_memalign((void**)&priv->boundary_val_base[i], 32,
+                             sizeof(*priv->boundary_val_base[i]) * boundary_len_max * FD_STENCIL_MAX);
         if (ret != 0)
             return -ret;
-        ctx->boundaries[i].val_len = boundary_size;
+
+        ctx->boundaries[i].val        = priv->boundary_val_base[i] + FD_STENCIL_MAX - 1;
+        ctx->boundaries[i].val_stride = boundary_len_max;
     }
 
     return 0;
@@ -399,10 +399,10 @@ void mg2di_ell_relax_free(EllRelaxContext **pctx)
     free(ctx->priv->residual_base);
     for (int i = 0; i < ARRAY_ELEMS(ctx->priv->diff_coeffs_base); i++)
         free(ctx->priv->diff_coeffs_base[i]);
-    free(ctx->priv);
+    for (int i = 0; i < ARRAY_ELEMS(ctx->priv->boundary_val_base); i++)
+        free(ctx->priv->boundary_val_base[i]);
 
-    for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++)
-        free(ctx->boundaries[i].val);
+    free(ctx->priv);
 
     free(ctx);
     *pctx = NULL;
diff --git a/ell_relax.h b/ell_relax.h
index 0df42d5..0e32015 100644
--- a/ell_relax.h
+++ b/ell_relax.h
@@ -127,14 +127,20 @@ typedef struct EllRelaxBoundary {
     /**
      * For type = ELL_RELAX_BC_TYPE_FIXVAL:
      *     Values of the unknown function on the boundary.
+     *     The number of boundary layers is equal to fd_stencil.
+     *     The first boundary layer has the number of points equal to the
+     *     corresponding domain_size. Each subsequent boundary layer has one
+     *     more boundary point at each end of the domain.
      *
      * Ignored otherwise.
      */
     double *val;
     /**
-     * Number of elements in val.
+     * Number of elements between rows in val. I.e. if val[0] is the first
+     * boundary point, then val[val_stride - 1] is the first boundary point in
+     * the second row and so on.
      */
-    size_t  val_len;
+    size_t  val_stride;
 } EllRelaxBoundary;
 
 typedef struct EllRelaxContext {
diff --git a/libmg2d.v b/libmg2d.v
index 50ee8ca..271d7eb 100644
--- a/libmg2d.v
+++ b/libmg2d.v
@@ -1,4 +1,4 @@
-LIBMG2D_0 {
+LIBMG2D_1 {
     global: mg2d_*;
     local: *;
 };
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;
 }
diff --git a/mg2d.h b/mg2d.h
index a0c0c70..4ccb5da 100644
--- a/mg2d.h
+++ b/mg2d.h
@@ -108,14 +108,20 @@ typedef struct MG2DBoundary {
     /**
      * For type = MG2D_BC_TYPE_FIXVAL:
      *     Values of the unknown function on the boundary.
+     *     The number of boundary layers is equal to fd_stencil.
+     *     The first boundary layer has the number of points equal to the
+     *     corresponding domain_size. Each subsequent boundary layer has one
+     *     more boundary point at each end of the domain.
      *
      * Ignored otherwise.
      */
     double *val;
     /**
-     * Number of elements in val.
+     * Number of elements between rows in val. I.e. if val[0] is the first
+     * boundary point, then val[val_stride - 1] is the first boundary point in
+     * the second row and so on.
      */
-    size_t  val_len;
+    size_t  val_stride;
 } MG2DBoundary;
 
 typedef struct MG2DContext {