mg2d: add higher-order interpolation operators

1 files changed, 75 insertions, 0 deletions

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,