From c4bf2e6d39784b3357ac6e0239220a3213d33e47 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Sat, 18 Jan 2020 22:03:33 +0100 Subject: Add weighted restriction. --- transfer.c | 194 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--- transfer.h | 4 ++ 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, }; /** -- cgit v1.2.3