summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAnton Khirnov <anton@khirnov.net>2020-01-18 22:03:33 +0100
committerAnton Khirnov <anton@khirnov.net>2020-01-26 20:15:47 +0100
commitc4bf2e6d39784b3357ac6e0239220a3213d33e47 (patch)
treeeea2e55fc2a8d1651002af67a7052a8b9e2a1dba
parent62c92eea4ec4859fff5002931e2d7d562b3deb5d (diff)
Add weighted restriction.fw_generic
-rw-r--r--transfer.c194
-rw-r--r--transfer.h4
2 files changed, 189 insertions, 9 deletions
diff --git a/transfer.c b/transfer.c
index a791467..aa26cd1 100644
--- a/transfer.c
+++ b/transfer.c
@@ -112,6 +112,22 @@ void mg2di_transfer_interp2d_line_cont_6_fma3(double *dst, ptrdiff_t dst_len,
#undef CONTIGUOUS
#undef STENCIL
+// compute the weight for src[src_idx]'s contribution in the src_len-point
+// lagrange interpolation of point dst
+static double lagrange_weight(const double dst, const double *src,
+ const size_t src_len, const ptrdiff_t src_idx)
+{
+ double ret = 1.0;
+ for (unsigned int i = 0; i < src_len; i++) {
+ if (i == src_idx)
+ continue;
+
+ ret *= (dst - src[i]) / (src[src_idx] - src[i]);
+ }
+
+ return ret;
+}
+
// generate the interpolation source indices and weights
static int transfer_factors_lagrange(GridTransferContext *ctx, unsigned int dim)
{
@@ -135,15 +151,8 @@ static int transfer_factors_lagrange(GridTransferContext *ctx, unsigned int dim)
for (int i = 0; i < priv->stencil; i++)
coord_src[i] = (idx_src + ctx->src.start[dim] + i) * step_src;
- for (int i = 0; i < priv->stencil; i++) {
- fact[i] = 1.0;
- for (int j = 0; j < priv->stencil; j++) {
- if (i == j)
- continue;
-
- fact[i] *= (coord_dst - coord_src[j]) / (coord_src[i] - coord_src[j]);
- }
- }
+ for (int i = 0; i < priv->stencil; i++)
+ fact[i] = lagrange_weight(coord_dst, coord_src, priv->stencil, i);
priv->idx[dim][idx_dst] = idx_src;
}
@@ -397,6 +406,169 @@ static const GridTransfer transfer_fw = {
.transfer = transfer_fw_transfer,
};
+typedef struct GridTransferFWGeneric {
+ // order of the coarse-to-fine Lagrange interpolation that is adjoint to
+ // this fine-to-coarse weighted restriction
+ unsigned int order;
+
+ // number of fine points used for weighting in each direction
+ unsigned int stencil[2];
+ double scale[2];
+ ptrdiff_t **idx;
+ double **fact;
+} GridTransferFWGeneric;
+
+static int transfer_fw_generic_transfer(GridTransferContext *ctx,
+ NDArray *a_dst, const NDArray *a_src)
+{
+ GridTransferFWGeneric *priv = ctx->priv;
+
+ const double scale = 1. / (priv->scale[0] * priv->scale[1]);
+
+ const ptrdiff_t ss0 = a_src->stride[0];
+ const ptrdiff_t ss1 = a_src->stride[1];
+
+ const ptrdiff_t ds1 = a_dst->stride[1];
+
+ for (size_t y = 0; y < a_dst->shape[0]; y++) {
+ double *dst = a_dst->data + a_dst->stride[0] * y;
+
+ for (size_t x = 0; x < a_dst->shape[1]; x++) {
+ const double *src = a_src->data + priv->idx[0][y] * ss0 + priv->idx[1][x] * ss1;
+
+ double val = 0.0;
+
+ for (ptrdiff_t idx0 = 0; idx0 < priv->stencil[0]; idx0++) {
+ double tmp = 0.0;
+ for (ptrdiff_t idx1 = 0; idx1 < priv->stencil[1]; idx1++)
+ tmp += src[idx0 * ss0 + idx1 * ss1] * priv->fact[1][priv->stencil[1] * x + idx1];
+ val += tmp * priv->fact[0][priv->stencil[0] * y + idx0];
+ }
+
+ dst[x * ds1] = val * scale;
+ }
+ }
+
+ return 0;
+}
+
+static int transfer_factors_fw(GridTransferContext *ctx, unsigned int dim)
+{
+ GridTransferFWGeneric *priv = ctx->priv;
+
+ const double step_dst = ctx->dst.step[dim];
+ const double step_src = ctx->src.step[dim];
+
+ double *coord_coarse;
+
+ coord_coarse = calloc(priv->order + 1, sizeof(*coord_coarse));
+ if (!coord_coarse)
+ return -ENOMEM;
+
+ for (ptrdiff_t idx_dst = 0; idx_dst < ctx->dst.size[dim]; idx_dst++) {
+ const ptrdiff_t idx_src = (ptrdiff_t)(priv->scale[dim] * (idx_dst + ctx->dst.start[dim])) - ctx->src.start[dim] - (priv->stencil[dim] / 2 - 1);
+ double *fact = priv->fact[dim] + priv->stencil[dim] * idx_dst;
+
+ for (int i = 0; i < priv->stencil[dim]; i++) {
+ const double coord_src = (idx_src + ctx->src.start[dim] + i) * step_src;
+ const ptrdiff_t idx_coarse = (ptrdiff_t)floor(coord_src / step_dst) - ctx->dst.start[dim] - ((priv->order + 1) / 2 - 1);
+
+ int idx = -1;
+
+ for (int j = 0; j < priv->order + 1; j++) {
+ coord_coarse[j] = (idx_coarse + j + ctx->dst.start[dim]) * step_dst;
+ if (idx_coarse + j == idx_dst)
+ idx = j;
+ }
+
+ //if (coord_dst > coord_src) {
+ // coord_coarse[0] = coord_dst - step_dst;
+ // coord_coarse[1] = coord_dst;
+ // idx = 1;
+ //} else {
+ // coord_coarse[0] = coord_dst;
+ // coord_coarse[1] = coord_dst + step_dst;
+ // idx = 0;
+ //}
+
+ fact[i] = idx >= 0 ? lagrange_weight(coord_src, coord_coarse, priv->order + 1, idx) : 0.0;
+ }
+
+ priv->idx[dim][idx_dst] = idx_src;
+ }
+
+ free(coord_coarse);
+
+ return 0;
+}
+
+static void transfer_fw_generic_free(GridTransferContext *ctx)
+{
+ GridTransferFWGeneric *priv = ctx->priv;
+
+ for (int dim = 0; dim < ctx->nb_dims; dim++) {
+ if (priv->idx)
+ free(priv->idx[dim]);
+ if (priv->fact)
+ free(priv->fact[dim]);
+ }
+ free(priv->idx);
+ free(priv->fact);
+}
+static int transfer_fw_generic_init(GridTransferContext *ctx)
+{
+ GridTransferFWGeneric *priv = ctx->priv;
+ int ret;
+
+ if (ctx->nb_dims != 2)
+ return -ENOSYS;
+
+ // destination must be coarser than source
+ for (int i = 0; i < ctx->nb_dims; i++) {
+ if (ctx->dst.step[i] <= ctx->src.step[i])
+ return -EINVAL;
+ }
+
+ switch (ctx->op) {
+ case GRID_TRANSFER_FW_1_GENERIC: priv->order = 1; break;
+ case GRID_TRANSFER_FW_3_GENERIC: priv->order = 3; break;
+ case GRID_TRANSFER_FW_5_GENERIC: priv->order = 5; break;
+ default:
+ return -ENOSYS;
+ }
+
+ priv->idx = calloc(ctx->nb_dims, sizeof(*priv->idx));
+ priv->fact = calloc(ctx->nb_dims, sizeof(*priv->fact));
+ if (!priv->idx || !priv->fact)
+ return -ENOMEM;
+
+ for (int dim = 0; dim < ctx->nb_dims; dim++) {
+ priv->scale[dim] = ctx->dst.step[dim] / ctx->src.step[dim];
+ priv->stencil[dim] = ceil((priv->order + 1) * priv->scale[dim]);
+
+ priv->idx[dim] = calloc(ctx->dst.size[dim], sizeof(*priv->idx[dim]));
+ if (!priv->idx[dim])
+ return -ENOMEM;
+
+ priv->fact[dim] = calloc(ctx->dst.size[dim], sizeof(*priv->fact[dim]) * priv->stencil[dim]);
+ if (!priv->fact[dim])
+ return -ENOMEM;
+
+ ret = transfer_factors_fw(ctx, dim);
+ if (ret < 0)
+ return ret;
+ }
+
+ return 0;
+}
+
+static const GridTransfer transfer_fw_generic = {
+ .priv_data_size = sizeof(GridTransferFWGeneric),
+ .init = transfer_fw_generic_init,
+ .free = transfer_fw_generic_free,
+ .transfer = transfer_fw_generic_transfer,
+};
+
static const GridTransfer *transfer_ops[] = {
[GRID_TRANSFER_LAGRANGE_1] = &transfer_lagrange,
[GRID_TRANSFER_LAGRANGE_3] = &transfer_lagrange,
@@ -405,6 +577,10 @@ static const GridTransfer *transfer_ops[] = {
[GRID_TRANSFER_FW_1] = &transfer_fw,
[GRID_TRANSFER_FW_2] = &transfer_fw,
[GRID_TRANSFER_FW_3] = &transfer_fw,
+
+ [GRID_TRANSFER_FW_1_GENERIC] = &transfer_fw_generic,
+ [GRID_TRANSFER_FW_3_GENERIC] = &transfer_fw_generic,
+ [GRID_TRANSFER_FW_5_GENERIC] = &transfer_fw_generic,
};
int mg2di_gt_transfer(GridTransferContext *ctx,
diff --git a/transfer.h b/transfer.h
index 5c434e3..20896e0 100644
--- a/transfer.h
+++ b/transfer.h
@@ -33,6 +33,10 @@ enum GridTransferOperator {
GRID_TRANSFER_FW_2,
GRID_TRANSFER_FW_3,
GRID_TRANSFER_FW_4,
+
+ GRID_TRANSFER_FW_1_GENERIC,
+ GRID_TRANSFER_FW_3_GENERIC,
+ GRID_TRANSFER_FW_5_GENERIC,
};
/**