Add a new separate module for grid transfers/interpolation.

1 files changed, 64 insertions, 0 deletions

diff --git a/transfer_interp_template.c b/transfer_interp_template.c
new file mode 100644
index 0000000..65ae98b
--- /dev/null
+++ b/transfer_interp_template.c
@@ -0,0 +1,64 @@
+/*
+ * Polynomial interpolation template
+ * Copyright 2019 Anton Khirnov <anton@khirnov.net>
+ *
+ * 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 <http://www.gnu.org/licenses/>.
+ */
+
+#if CONTIGUOUS
+# define CONT cont
+#else
+# define CONT generic
+#endif
+
+#define JOIN3(a, b, c) a ## _ ## b ## _ ## c
+#define FUNC2(name, cont, stencil) JOIN3(name, cont, stencil)
+#define FUNC(name) FUNC2(name, CONT, STENCIL)
+
+static void
+FUNC(interp_transfer_line)(double *dst, ptrdiff_t dst_len,
+                           const double *src, ptrdiff_t src_stride,
+                           const ptrdiff_t *idx_x, const double *fact_x, const double *fact_y
+#if !CONTIGUOUS
+                           , ptrdiff_t dst_stride0, ptrdiff_t src_stride0
+# define SSTRIDE1 src_stride0
+# define DSTRIDE1 dst_stride0
+#else
+# define SSTRIDE1 1
+# define DSTRIDE1 1
+#endif
+                           )
+{
+    for (size_t x = 0; x < dst_len; x++) {
+        const double *src_data = src + SSTRIDE1 * idx_x[x];
+
+        double val = 0.0;
+
+        for (ptrdiff_t idx0 = 0; idx0 < STENCIL; idx0++) {
+            double tmp = 0.0;
+            for (ptrdiff_t idx1 = 0; idx1 < STENCIL; idx1++)
+                tmp += src_data[idx0 * src_stride + idx1 ] * fact_x[STENCIL * x + idx1];
+            val += tmp * fact_y[idx0];
+        }
+
+        dst[x * DSTRIDE1] = val;
+    }
+}
+
+#undef SSTRIDE1
+#undef DSTRIDE1
+#undef CONT
+#undef JOIN3
+#undef FUNC2
+#undef FUNC