Mercurial > hg > sv-dependency-builds

annotate src/fftw-3.3.5/mpi/rdft2-rank-geq2-transposed.c @ 43:5ea0608b923f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Current zlib source
author Chris Cannam
date Tue, 18 Oct 2016 14:33:52 +0100
parents 2cd0e3b3e1fd
children
rev   line source
Chris@42 1 /*
Chris@42 2 * Copyright (c) 2003, 2007-14 Matteo Frigo
Chris@42 3 * Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology
Chris@42 4 *
Chris@42 5 * This program is free software; you can redistribute it and/or modify
Chris@42 6 * it under the terms of the GNU General Public License as published by
Chris@42 7 * the Free Software Foundation; either version 2 of the License, or
Chris@42 8 * (at your option) any later version.
Chris@42 9 *
Chris@42 10 * This program is distributed in the hope that it will be useful,
Chris@42 11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
Chris@42 12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
Chris@42 13 * GNU General Public License for more details.
Chris@42 14 *
Chris@42 15 * You should have received a copy of the GNU General Public License
Chris@42 16 * along with this program; if not, write to the Free Software
Chris@42 17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
Chris@42 18 *
Chris@42 19 */
Chris@42 20
Chris@42 21 /* Real-input (r2c) DFTs of rank >= 2, for the case where we are distributed
Chris@42 22 across the first dimension only, and the output is transposed both
Chris@42 23 in data distribution and in ordering (for the first 2 dimensions).
Chris@42 24
Chris@42 25 Conversely, real-output (c2r) DFTs where the input is transposed.
Chris@42 26
Chris@42 27 We don't currently support transposed-input r2c or transposed-output
Chris@42 28 c2r transforms. */
Chris@42 29
Chris@42 30 #include "mpi-rdft2.h"
Chris@42 31 #include "mpi-transpose.h"
Chris@42 32 #include "rdft.h"
Chris@42 33 #include "dft.h"
Chris@42 34
Chris@42 35 typedef struct {
Chris@42 36 solver super;
Chris@42 37 int preserve_input; /* preserve input even if DESTROY_INPUT was passed */
Chris@42 38 } S;
Chris@42 39
Chris@42 40 typedef struct {
Chris@42 41 plan_mpi_rdft2 super;
Chris@42 42
Chris@42 43 plan *cld1, *cldt, *cld2;
Chris@42 44 INT vn;
Chris@42 45 int preserve_input;
Chris@42 46 } P;
Chris@42 47
Chris@42 48 static void apply_r2c(const plan *ego_, R *I, R *O)
Chris@42 49 {
Chris@42 50 const P *ego = (const P *) ego_;
Chris@42 51 plan_rdft2 *cld1;
Chris@42 52 plan_dft *cld2;
Chris@42 53 plan_rdft *cldt;
Chris@42 54
Chris@42 55 /* RDFT2 local dimensions */
Chris@42 56 cld1 = (plan_rdft2 *) ego->cld1;
Chris@42 57 if (ego->preserve_input) {
Chris@42 58 cld1->apply(ego->cld1, I, I+ego->vn, O, O+1);
Chris@42 59 I = O;
Chris@42 60 }
Chris@42 61 else
Chris@42 62 cld1->apply(ego->cld1, I, I+ego->vn, I, I+1);
Chris@42 63
Chris@42 64 /* global transpose */
Chris@42 65 cldt = (plan_rdft *) ego->cldt;
Chris@42 66 cldt->apply(ego->cldt, I, O);
Chris@42 67
Chris@42 68 /* DFT final local dimension */
Chris@42 69 cld2 = (plan_dft *) ego->cld2;
Chris@42 70 cld2->apply(ego->cld2, O, O+1, O, O+1);
Chris@42 71 }
Chris@42 72
Chris@42 73 static void apply_c2r(const plan *ego_, R *I, R *O)
Chris@42 74 {
Chris@42 75 const P *ego = (const P *) ego_;
Chris@42 76 plan_rdft2 *cld1;
Chris@42 77 plan_dft *cld2;
Chris@42 78 plan_rdft *cldt;
Chris@42 79
Chris@42 80 /* IDFT local dimensions */
Chris@42 81 cld2 = (plan_dft *) ego->cld2;
Chris@42 82 if (ego->preserve_input) {
Chris@42 83 cld2->apply(ego->cld2, I+1, I, O+1, O);
Chris@42 84 I = O;
Chris@42 85 }
Chris@42 86 else
Chris@42 87 cld2->apply(ego->cld2, I+1, I, I+1, I);
Chris@42 88
Chris@42 89 /* global transpose */
Chris@42 90 cldt = (plan_rdft *) ego->cldt;
Chris@42 91 cldt->apply(ego->cldt, I, O);
Chris@42 92
Chris@42 93 /* RDFT2 final local dimension */
Chris@42 94 cld1 = (plan_rdft2 *) ego->cld1;
Chris@42 95 cld1->apply(ego->cld1, O, O+ego->vn, O, O+1);
Chris@42 96 }
Chris@42 97
Chris@42 98 static int applicable(const S *ego, const problem *p_,
Chris@42 99 const planner *plnr)
Chris@42 100 {
Chris@42 101 const problem_mpi_rdft2 *p = (const problem_mpi_rdft2 *) p_;
Chris@42 102 return (1
Chris@42 103 && p->sz->rnk > 1
Chris@42 104 && (!ego->preserve_input || (!NO_DESTROY_INPUTP(plnr)
Chris@42 105 && p->I != p->O))
Chris@42 106 && ((p->flags == TRANSPOSED_OUT && p->kind == R2HC
Chris@42 107 && XM(is_local_after)(1, p->sz, IB)
Chris@42 108 && XM(is_local_after)(2, p->sz, OB)
Chris@42 109 && XM(num_blocks)(p->sz->dims[0].n,
Chris@42 110 p->sz->dims[0].b[OB]) == 1)
Chris@42 111 ||
Chris@42 112 (p->flags == TRANSPOSED_IN && p->kind == HC2R
Chris@42 113 && XM(is_local_after)(1, p->sz, OB)
Chris@42 114 && XM(is_local_after)(2, p->sz, IB)
Chris@42 115 && XM(num_blocks)(p->sz->dims[0].n,
Chris@42 116 p->sz->dims[0].b[IB]) == 1))
Chris@42 117 && (!NO_SLOWP(plnr) /* slow if rdft2-serial is applicable */
Chris@42 118 || !XM(rdft2_serial_applicable)(p))
Chris@42 119 );
Chris@42 120 }
Chris@42 121
Chris@42 122 static void awake(plan *ego_, enum wakefulness wakefulness)
Chris@42 123 {
Chris@42 124 P *ego = (P *) ego_;
Chris@42 125 X(plan_awake)(ego->cld1, wakefulness);
Chris@42 126 X(plan_awake)(ego->cldt, wakefulness);
Chris@42 127 X(plan_awake)(ego->cld2, wakefulness);
Chris@42 128 }
Chris@42 129
Chris@42 130 static void destroy(plan *ego_)
Chris@42 131 {
Chris@42 132 P *ego = (P *) ego_;
Chris@42 133 X(plan_destroy_internal)(ego->cld2);
Chris@42 134 X(plan_destroy_internal)(ego->cldt);
Chris@42 135 X(plan_destroy_internal)(ego->cld1);
Chris@42 136 }
Chris@42 137
Chris@42 138 static void print(const plan *ego_, printer *p)
Chris@42 139 {
Chris@42 140 const P *ego = (const P *) ego_;
Chris@42 141 p->print(p, "(mpi-rdft2-rank-geq2-transposed%s%(%p%)%(%p%)%(%p%))",
Chris@42 142 ego->preserve_input==2 ?"/p":"",
Chris@42 143 ego->cld1, ego->cldt, ego->cld2);
Chris@42 144 }
Chris@42 145
Chris@42 146 static plan *mkplan(const solver *ego_, const problem *p_, planner *plnr)
Chris@42 147 {
Chris@42 148 const S *ego = (const S *) ego_;
Chris@42 149 const problem_mpi_rdft2 *p;
Chris@42 150 P *pln;
Chris@42 151 plan *cld1 = 0, *cldt = 0, *cld2 = 0;
Chris@42 152 R *r0, *r1, *cr, *ci, *ri, *ii, *ro, *io, *I, *O;
Chris@42 153 tensor *sz;
Chris@42 154 int i, my_pe, n_pes;
Chris@42 155 INT nrest, n1, b1;
Chris@42 156 static const plan_adt padt = {
Chris@42 157 XM(rdft2_solve), awake, print, destroy
Chris@42 158 };
Chris@42 159 block_kind k1, k2;
Chris@42 160
Chris@42 161 UNUSED(ego);
Chris@42 162
Chris@42 163 if (!applicable(ego, p_, plnr))
Chris@42 164 return (plan *) 0;
Chris@42 165
Chris@42 166 p = (const problem_mpi_rdft2 *) p_;
Chris@42 167
Chris@42 168 I = p->I; O = p->O;
Chris@42 169 if (p->kind == R2HC) {
Chris@42 170 k1 = IB; k2 = OB;
Chris@42 171 r1 = (r0 = I) + p->vn;
Chris@42 172 if (ego->preserve_input || NO_DESTROY_INPUTP(plnr)) {
Chris@42 173 ci = (cr = O) + 1;
Chris@42 174 I = O;
Chris@42 175 }
Chris@42 176 else
Chris@42 177 ci = (cr = I) + 1;
Chris@42 178 io = ii = (ro = ri = O) + 1;
Chris@42 179 }
Chris@42 180 else {
Chris@42 181 k1 = OB; k2 = IB;
Chris@42 182 r1 = (r0 = O) + p->vn;
Chris@42 183 ci = (cr = O) + 1;
Chris@42 184 if (ego->preserve_input || NO_DESTROY_INPUTP(plnr)) {
Chris@42 185 ri = (ii = I) + 1;
Chris@42 186 ro = (io = O) + 1;
Chris@42 187 I = O;
Chris@42 188 }
Chris@42 189 else
Chris@42 190 ro = ri = (io = ii = I) + 1;
Chris@42 191 }
Chris@42 192
Chris@42 193 MPI_Comm_rank(p->comm, &my_pe);
Chris@42 194 MPI_Comm_size(p->comm, &n_pes);
Chris@42 195
Chris@42 196 sz = X(mktensor)(p->sz->rnk - 1); /* tensor of last rnk-1 dimensions */
Chris@42 197 i = p->sz->rnk - 2; A(i >= 0);
Chris@42 198 sz->dims[i].n = p->sz->dims[i+1].n / 2 + 1;
Chris@42 199 sz->dims[i].is = sz->dims[i].os = 2 * p->vn;
Chris@42 200 for (--i; i >= 0; --i) {
Chris@42 201 sz->dims[i].n = p->sz->dims[i+1].n;
Chris@42 202 sz->dims[i].is = sz->dims[i].os = sz->dims[i+1].n * sz->dims[i+1].is;
Chris@42 203 }
Chris@42 204 nrest = 1; for (i = 1; i < sz->rnk; ++i) nrest *= sz->dims[i].n;
Chris@42 205 {
Chris@42 206 INT ivs = 1 + (p->kind == HC2R), ovs = 1 + (p->kind == R2HC);
Chris@42 207 INT is = sz->dims[0].n * sz->dims[0].is;
Chris@42 208 INT b = XM(block)(p->sz->dims[0].n, p->sz->dims[0].b[k1], my_pe);
Chris@42 209 sz->dims[p->sz->rnk - 2].n = p->sz->dims[p->sz->rnk - 1].n;
Chris@42 210 cld1 = X(mkplan_d)(plnr,
Chris@42 211 X(mkproblem_rdft2_d)(sz,
Chris@42 212 X(mktensor_2d)(b, is, is,
Chris@42 213 p->vn,ivs,ovs),
Chris@42 214 r0, r1, cr, ci, p->kind));
Chris@42 215 if (XM(any_true)(!cld1, p->comm)) goto nada;
Chris@42 216 }
Chris@42 217
Chris@42 218 nrest *= p->vn;
Chris@42 219 n1 = p->sz->dims[1].n;
Chris@42 220 b1 = p->sz->dims[1].b[k2];
Chris@42 221 if (p->sz->rnk == 2) { /* n1 dimension is cut in ~half */
Chris@42 222 n1 = n1 / 2 + 1;
Chris@42 223 b1 = b1 == p->sz->dims[1].n ? n1 : b1;
Chris@42 224 }
Chris@42 225
Chris@42 226 if (p->kind == R2HC)
Chris@42 227 cldt = X(mkplan_d)(plnr,
Chris@42 228 XM(mkproblem_transpose)(
Chris@42 229 p->sz->dims[0].n, n1, nrest * 2,
Chris@42 230 I, O,
Chris@42 231 p->sz->dims[0].b[IB], b1,
Chris@42 232 p->comm, 0));
Chris@42 233 else
Chris@42 234 cldt = X(mkplan_d)(plnr,
Chris@42 235 XM(mkproblem_transpose)(
Chris@42 236 n1, p->sz->dims[0].n, nrest * 2,
Chris@42 237 I, O,
Chris@42 238 b1, p->sz->dims[0].b[OB],
Chris@42 239 p->comm, 0));
Chris@42 240 if (XM(any_true)(!cldt, p->comm)) goto nada;
Chris@42 241
Chris@42 242 {
Chris@42 243 INT is = p->sz->dims[0].n * nrest * 2;
Chris@42 244 INT b = XM(block)(n1, b1, my_pe);
Chris@42 245 cld2 = X(mkplan_d)(plnr,
Chris@42 246 X(mkproblem_dft_d)(X(mktensor_1d)(
Chris@42 247 p->sz->dims[0].n,
Chris@42 248 nrest * 2, nrest * 2),
Chris@42 249 X(mktensor_2d)(b, is, is,
Chris@42 250 nrest, 2, 2),
Chris@42 251 ri, ii, ro, io));
Chris@42 252 if (XM(any_true)(!cld2, p->comm)) goto nada;
Chris@42 253 }
Chris@42 254
Chris@42 255 pln = MKPLAN_MPI_RDFT2(P, &padt, p->kind == R2HC ? apply_r2c : apply_c2r);
Chris@42 256 pln->cld1 = cld1;
Chris@42 257 pln->cldt = cldt;
Chris@42 258 pln->cld2 = cld2;
Chris@42 259 pln->preserve_input = ego->preserve_input ? 2 : NO_DESTROY_INPUTP(plnr);
Chris@42 260 pln->vn = p->vn;
Chris@42 261
Chris@42 262 X(ops_add)(&cld1->ops, &cld2->ops, &pln->super.super.ops);
Chris@42 263 X(ops_add2)(&cldt->ops, &pln->super.super.ops);
Chris@42 264
Chris@42 265 return &(pln->super.super);
Chris@42 266
Chris@42 267 nada:
Chris@42 268 X(plan_destroy_internal)(cld2);
Chris@42 269 X(plan_destroy_internal)(cldt);
Chris@42 270 X(plan_destroy_internal)(cld1);
Chris@42 271 return (plan *) 0;
Chris@42 272 }
Chris@42 273
Chris@42 274 static solver *mksolver(int preserve_input)
Chris@42 275 {
Chris@42 276 static const solver_adt sadt = { PROBLEM_MPI_RDFT2, mkplan, 0 };
Chris@42 277 S *slv = MKSOLVER(S, &sadt);
Chris@42 278 slv->preserve_input = preserve_input;
Chris@42 279 return &(slv->super);
Chris@42 280 }
Chris@42 281
Chris@42 282 void XM(rdft2_rank_geq2_transposed_register)(planner *p)
Chris@42 283 {
Chris@42 284 int preserve_input;
Chris@42 285 for (preserve_input = 0; preserve_input <= 1; ++preserve_input)
Chris@42 286 REGISTER_SOLVER(p, mksolver(preserve_input));
Chris@42 287 }