From 1e83fd63c30d433ee53769d4e4768feabf822ae2 Mon Sep 17 00:00:00 2001
From: Anton Khirnov <anton@khirnov.net>
Date: Tue, 9 Apr 2019 09:44:46 +0200
Subject: mg2d: add higher-order interpolation operators

---
 common.h   |  2 +-
 mg2d.c     | 13 ++++++-----
 transfer.c | 75 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 transfer.h |  3 +++
 4 files changed, 86 insertions(+), 7 deletions(-)

diff --git a/common.h b/common.h
index 7ecdf01..0a36ca9 100644
--- a/common.h
+++ b/common.h
@@ -34,6 +34,6 @@ static inline int64_t gettime(void)
     return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec;
 }
 
-#define FD_STENCIL_MAX 2
+#define FD_STENCIL_MAX 4
 
 #endif /* MG2D_COMMON_H */
diff --git a/mg2d.c b/mg2d.c
index b1c7fb7..11a8158 100644
--- a/mg2d.c
+++ b/mg2d.c
@@ -396,16 +396,17 @@ static int mg_levels_init(MG2DContext *ctx)
             enum GridTransferOperator op_interp;
             enum GridTransferOperator op_restrict;
 
-            if (ctx->fd_stencil == 1)
-                op_interp = GRID_TRANSFER_LAGRANGE_1;
-            else
-                op_interp = GRID_TRANSFER_LAGRANGE_3;
+            switch (ctx->fd_stencil) {
+            case 1: op_interp = GRID_TRANSFER_LAGRANGE_3; break;
+            case 2: op_interp = GRID_TRANSFER_LAGRANGE_5; break;
+            default: return -ENOSYS;
+            }
 
             if (cur->solver->domain_size[0] == 2 * (cur->child->solver->domain_size[0] - 1) + 1) {
                 if (ctx->fd_stencil == 1)
-                    op_restrict = GRID_TRANSFER_FW_1;
-                else
                     op_restrict = GRID_TRANSFER_FW_3;
+                else
+                    op_restrict = GRID_TRANSFER_FW_5;
             } else
                 op_restrict = op_interp;
 
diff --git a/transfer.c b/transfer.c
index 232d6e1..311c6a9 100644
--- a/transfer.c
+++ b/transfer.c
@@ -83,6 +83,26 @@ void mg2di_transfer_interp_line_cont_4_fma3(double *dst, ptrdiff_t dst_len,
 #undef CONTIGUOUS
 #undef STENCIL
 
+#define STENCIL 6
+#define CONTIGUOUS 0
+#include "transfer_interp_template.c"
+#undef CONTIGUOUS
+
+#define CONTIGUOUS 1
+#include "transfer_interp_template.c"
+#undef CONTIGUOUS
+#undef STENCIL
+
+#define STENCIL 8
+#define CONTIGUOUS 0
+#include "transfer_interp_template.c"
+#undef CONTIGUOUS
+
+#define CONTIGUOUS 1
+#include "transfer_interp_template.c"
+#undef CONTIGUOUS
+#undef STENCIL
+
 // generate the interpolation source indices and weights
 static int transfer_factors_lagrange(GridTransferContext *ctx, unsigned int dim)
 {
@@ -146,6 +166,16 @@ static int transfer_lagrange_init(GridTransferContext *ctx)
         }
 #endif
         break;
+    case GRID_TRANSFER_LAGRANGE_5:
+        priv->transfer_cont    = interp_transfer_line_cont_6;
+        priv->transfer_generic = interp_transfer_line_generic_6;
+        priv->stencil = 6;
+        break;
+    case GRID_TRANSFER_LAGRANGE_7:
+        priv->transfer_cont    = interp_transfer_line_cont_8;
+        priv->transfer_generic = interp_transfer_line_generic_8;
+        priv->stencil = 8;
+        break;
     default:
         return -EINVAL;
     }
@@ -312,11 +342,56 @@ static const GridTransfer transfer_fw_3 = {
     .transfer       = transfer_fw_3_transfer,
 };
 
+static int fw_5_transfer_task(void *arg, unsigned int job_idx, unsigned int thread_idx)
+{
+    static const double factors[] = { 3.0 / 256, -25.0 / 256, 75.0 / 128, 1.0, 75.0 / 128, -25.0 / 256, 3.0 / 256 };
+    static const size_t order = ARRAY_ELEMS(factors) / 2;
+    FWThreadData *td = arg;
+
+    const ptrdiff_t ss0 = td->src->stride[0];
+    const ptrdiff_t ss1 = td->src->stride[1];
+    const double    *sd = td->src->data + (job_idx * 2 - order) * ss0 - order * ss1;
+
+    double          *dd = td->dst->data + job_idx * td->dst->stride[0];
+    const ptrdiff_t ds1 = td->dst->stride[1];
+
+    for (size_t x = 0; x < td->dst->shape[1]; x++) {
+        double val = 0.0;
+
+        for (int i = 0; i < ARRAY_ELEMS(factors); i++)
+            for (int j = 0; j < ARRAY_ELEMS(factors); j++)
+                val += sd[i * ss0 + j * ss1] * factors[i] * factors[j];
+
+        *dd = 0.25 * val;
+        dd += ds1;
+        sd += ss1 * 2;
+    }
+
+    return 0;
+}
+static int transfer_fw_5_transfer(GridTransferContext *ctx,
+                                  NDArray *dst, const NDArray *src)
+{
+    FWThreadData td = { dst, src };
+
+    tp_execute(ctx->tp, ctx->dst.size[0], fw_5_transfer_task, &td);
+
+    return 0;
+}
+
+static const GridTransfer transfer_fw_5 = {
+    .priv_data_size = sizeof(GridTransferLagrange),
+    .transfer       = transfer_fw_5_transfer,
+};
+
 static const GridTransfer *transfer_ops[] = {
     [GRID_TRANSFER_LAGRANGE_1] = &transfer_lagrange,
     [GRID_TRANSFER_LAGRANGE_3] = &transfer_lagrange,
+    [GRID_TRANSFER_LAGRANGE_5] = &transfer_lagrange,
+    [GRID_TRANSFER_LAGRANGE_7] = &transfer_lagrange,
     [GRID_TRANSFER_FW_1]       = &transfer_fw_1,
     [GRID_TRANSFER_FW_3]       = &transfer_fw_3,
+    [GRID_TRANSFER_FW_5]       = &transfer_fw_5,
 };
 
 int mg2di_gt_transfer(GridTransferContext *ctx,
diff --git a/transfer.h b/transfer.h
index 424ca64..b2133f4 100644
--- a/transfer.h
+++ b/transfer.h
@@ -27,8 +27,11 @@
 enum GridTransferOperator {
     GRID_TRANSFER_LAGRANGE_1,
     GRID_TRANSFER_LAGRANGE_3,
+    GRID_TRANSFER_LAGRANGE_5,
+    GRID_TRANSFER_LAGRANGE_7,
     GRID_TRANSFER_FW_1,
     GRID_TRANSFER_FW_3,
+    GRID_TRANSFER_FW_5,
 };
 
 typedef struct RegularGrid2D {
-- 
cgit v1.2.3