Chris@82: /* Chris@82: * Copyright (c) 2003, 2007-14 Matteo Frigo Chris@82: * Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology Chris@82: * Chris@82: * This program is free software; you can redistribute it and/or modify Chris@82: * it under the terms of the GNU General Public License as published by Chris@82: * the Free Software Foundation; either version 2 of the License, or Chris@82: * (at your option) any later version. Chris@82: * Chris@82: * This program is distributed in the hope that it will be useful, Chris@82: * but WITHOUT ANY WARRANTY; without even the implied warranty of Chris@82: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the Chris@82: * GNU General Public License for more details. Chris@82: * Chris@82: * You should have received a copy of the GNU General Public License Chris@82: * along with this program; if not, write to the Free Software Chris@82: * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA Chris@82: * Chris@82: */ Chris@82: Chris@82: #include "api/api.h" Chris@82: #include "fftw3-mpi.h" Chris@82: #include "ifftw-mpi.h" Chris@82: #include "mpi-transpose.h" Chris@82: #include "mpi-dft.h" Chris@82: #include "mpi-rdft.h" Chris@82: #include "mpi-rdft2.h" Chris@82: Chris@82: /* Convert API flags to internal MPI flags. */ Chris@82: #define MPI_FLAGS(f) ((f) >> 27) Chris@82: Chris@82: /*************************************************************************/ Chris@82: Chris@82: static int mpi_inited = 0; Chris@82: Chris@82: static MPI_Comm problem_comm(const problem *p) { Chris@82: switch (p->adt->problem_kind) { Chris@82: case PROBLEM_MPI_DFT: Chris@82: return ((const problem_mpi_dft *) p)->comm; Chris@82: case PROBLEM_MPI_RDFT: Chris@82: return ((const problem_mpi_rdft *) p)->comm; Chris@82: case PROBLEM_MPI_RDFT2: Chris@82: return ((const problem_mpi_rdft2 *) p)->comm; Chris@82: case PROBLEM_MPI_TRANSPOSE: Chris@82: return ((const problem_mpi_transpose *) p)->comm; Chris@82: default: Chris@82: return MPI_COMM_NULL; Chris@82: } Chris@82: } Chris@82: Chris@82: /* used to synchronize cost measurements (timing or estimation) Chris@82: across all processes for an MPI problem, which is critical to Chris@82: ensure that all processes decide to use the same MPI plans Chris@82: (whereas serial plans need not be syncronized). */ Chris@82: static double cost_hook(const problem *p, double t, cost_kind k) Chris@82: { Chris@82: MPI_Comm comm = problem_comm(p); Chris@82: double tsum; Chris@82: if (comm == MPI_COMM_NULL) return t; Chris@82: MPI_Allreduce(&t, &tsum, 1, MPI_DOUBLE, Chris@82: k == COST_SUM ? MPI_SUM : MPI_MAX, comm); Chris@82: return tsum; Chris@82: } Chris@82: Chris@82: /* Used to reject wisdom that is not in sync across all processes Chris@82: for an MPI problem, which is critical to ensure that all processes Chris@82: decide to use the same MPI plans. (Even though costs are synchronized, Chris@82: above, out-of-sync wisdom may result from plans being produced Chris@82: by communicators that do not span all processes, either from a Chris@82: user-specified communicator or e.g. from transpose-recurse. */ Chris@82: static int wisdom_ok_hook(const problem *p, flags_t flags) Chris@82: { Chris@82: MPI_Comm comm = problem_comm(p); Chris@82: int eq_me, eq_all; Chris@82: /* unpack flags bitfield, since MPI communications may involve Chris@82: byte-order changes and MPI cannot do this for bit fields */ Chris@82: #if SIZEOF_UNSIGNED_INT >= 4 /* must be big enough to hold 20-bit fields */ Chris@82: unsigned int f[5]; Chris@82: #else Chris@82: unsigned long f[5]; /* at least 32 bits as per C standard */ Chris@82: #endif Chris@82: Chris@82: if (comm == MPI_COMM_NULL) return 1; /* non-MPI wisdom is always ok */ Chris@82: Chris@82: if (XM(any_true)(0, comm)) return 0; /* some process had nowisdom_hook */ Chris@82: Chris@82: /* otherwise, check that the flags and solver index are identical Chris@82: on all processes in this problem's communicator. Chris@82: Chris@82: TO DO: possibly we can relax strict equality, but it is Chris@82: critical to ensure that any flags which affect what plan is Chris@82: created (and whether the solver is applicable) are the same, Chris@82: e.g. DESTROY_INPUT, NO_UGLY, etcetera. (If the MPI algorithm Chris@82: differs between processes, deadlocks/crashes generally result.) */ Chris@82: f[0] = flags.l; Chris@82: f[1] = flags.hash_info; Chris@82: f[2] = flags.timelimit_impatience; Chris@82: f[3] = flags.u; Chris@82: f[4] = flags.slvndx; Chris@82: MPI_Bcast(f, 5, Chris@82: SIZEOF_UNSIGNED_INT >= 4 ? MPI_UNSIGNED : MPI_UNSIGNED_LONG, Chris@82: 0, comm); Chris@82: eq_me = f[0] == flags.l && f[1] == flags.hash_info Chris@82: && f[2] == flags.timelimit_impatience Chris@82: && f[3] == flags.u && f[4] == flags.slvndx; Chris@82: MPI_Allreduce(&eq_me, &eq_all, 1, MPI_INT, MPI_LAND, comm); Chris@82: return eq_all; Chris@82: } Chris@82: Chris@82: /* This hook is called when wisdom is not found. The any_true here Chris@82: matches up with the any_true in wisdom_ok_hook, in order to handle Chris@82: the case where some processes had wisdom (and called wisdom_ok_hook) Chris@82: and some processes didn't have wisdom (and called nowisdom_hook). */ Chris@82: static void nowisdom_hook(const problem *p) Chris@82: { Chris@82: MPI_Comm comm = problem_comm(p); Chris@82: if (comm == MPI_COMM_NULL) return; /* nothing to do for non-MPI p */ Chris@82: XM(any_true)(1, comm); /* signal nowisdom to any wisdom_ok_hook */ Chris@82: } Chris@82: Chris@82: /* needed to synchronize planner bogosity flag, in case non-MPI problems Chris@82: on a subset of processes encountered bogus wisdom */ Chris@82: static wisdom_state_t bogosity_hook(wisdom_state_t state, const problem *p) Chris@82: { Chris@82: MPI_Comm comm = problem_comm(p); Chris@82: if (comm != MPI_COMM_NULL /* an MPI problem */ Chris@82: && XM(any_true)(state == WISDOM_IS_BOGUS, comm)) /* bogus somewhere */ Chris@82: return WISDOM_IS_BOGUS; Chris@82: return state; Chris@82: } Chris@82: Chris@82: void XM(init)(void) Chris@82: { Chris@82: if (!mpi_inited) { Chris@82: planner *plnr = X(the_planner)(); Chris@82: plnr->cost_hook = cost_hook; Chris@82: plnr->wisdom_ok_hook = wisdom_ok_hook; Chris@82: plnr->nowisdom_hook = nowisdom_hook; Chris@82: plnr->bogosity_hook = bogosity_hook; Chris@82: XM(conf_standard)(plnr); Chris@82: mpi_inited = 1; Chris@82: } Chris@82: } Chris@82: Chris@82: void XM(cleanup)(void) Chris@82: { Chris@82: X(cleanup)(); Chris@82: mpi_inited = 0; Chris@82: } Chris@82: Chris@82: /*************************************************************************/ Chris@82: Chris@82: static dtensor *mkdtensor_api(int rnk, const XM(ddim) *dims0) Chris@82: { Chris@82: dtensor *x = XM(mkdtensor)(rnk); Chris@82: int i; Chris@82: for (i = 0; i < rnk; ++i) { Chris@82: x->dims[i].n = dims0[i].n; Chris@82: x->dims[i].b[IB] = dims0[i].ib; Chris@82: x->dims[i].b[OB] = dims0[i].ob; Chris@82: } Chris@82: return x; Chris@82: } Chris@82: Chris@82: static dtensor *default_sz(int rnk, const XM(ddim) *dims0, int n_pes, Chris@82: int rdft2) Chris@82: { Chris@82: dtensor *sz = XM(mkdtensor)(rnk); Chris@82: dtensor *sz0 = mkdtensor_api(rnk, dims0); Chris@82: block_kind k; Chris@82: int i; Chris@82: Chris@82: for (i = 0; i < rnk; ++i) Chris@82: sz->dims[i].n = dims0[i].n; Chris@82: Chris@82: if (rdft2) sz->dims[rnk-1].n = dims0[rnk-1].n / 2 + 1; Chris@82: Chris@82: for (i = 0; i < rnk; ++i) { Chris@82: sz->dims[i].b[IB] = dims0[i].ib ? dims0[i].ib : sz->dims[i].n; Chris@82: sz->dims[i].b[OB] = dims0[i].ob ? dims0[i].ob : sz->dims[i].n; Chris@82: } Chris@82: Chris@82: /* If we haven't used all of the processes yet, and some of the Chris@82: block sizes weren't specified (i.e. 0), then set the Chris@82: unspecified blocks so as to use as many processes as Chris@82: possible with as few distributed dimensions as possible. */ Chris@82: FORALL_BLOCK_KIND(k) { Chris@82: INT nb = XM(num_blocks_total)(sz, k); Chris@82: INT np = n_pes / nb; Chris@82: for (i = 0; i < rnk && np > 1; ++i) Chris@82: if (!sz0->dims[i].b[k]) { Chris@82: sz->dims[i].b[k] = XM(default_block)(sz->dims[i].n, np); Chris@82: nb *= XM(num_blocks)(sz->dims[i].n, sz->dims[i].b[k]); Chris@82: np = n_pes / nb; Chris@82: } Chris@82: } Chris@82: Chris@82: if (rdft2) sz->dims[rnk-1].n = dims0[rnk-1].n; Chris@82: Chris@82: /* punt for 1d prime */ Chris@82: if (rnk == 1 && X(is_prime)(sz->dims[0].n)) Chris@82: sz->dims[0].b[IB] = sz->dims[0].b[OB] = sz->dims[0].n; Chris@82: Chris@82: XM(dtensor_destroy)(sz0); Chris@82: sz0 = XM(dtensor_canonical)(sz, 0); Chris@82: XM(dtensor_destroy)(sz); Chris@82: return sz0; Chris@82: } Chris@82: Chris@82: /* allocate simple local (serial) dims array corresponding to n[rnk] */ Chris@82: static XM(ddim) *simple_dims(int rnk, const ptrdiff_t *n) Chris@82: { Chris@82: XM(ddim) *dims = (XM(ddim) *) MALLOC(sizeof(XM(ddim)) * rnk, Chris@82: TENSORS); Chris@82: int i; Chris@82: for (i = 0; i < rnk; ++i) Chris@82: dims[i].n = dims[i].ib = dims[i].ob = n[i]; Chris@82: return dims; Chris@82: } Chris@82: Chris@82: /*************************************************************************/ Chris@82: Chris@82: static void local_size(int my_pe, const dtensor *sz, block_kind k, Chris@82: ptrdiff_t *local_n, ptrdiff_t *local_start) Chris@82: { Chris@82: int i; Chris@82: if (my_pe >= XM(num_blocks_total)(sz, k)) Chris@82: for (i = 0; i < sz->rnk; ++i) Chris@82: local_n[i] = local_start[i] = 0; Chris@82: else { Chris@82: XM(block_coords)(sz, k, my_pe, local_start); Chris@82: for (i = 0; i < sz->rnk; ++i) { Chris@82: local_n[i] = XM(block)(sz->dims[i].n, sz->dims[i].b[k], Chris@82: local_start[i]); Chris@82: local_start[i] *= sz->dims[i].b[k]; Chris@82: } Chris@82: } Chris@82: } Chris@82: Chris@82: static INT prod(int rnk, const ptrdiff_t *local_n) Chris@82: { Chris@82: int i; Chris@82: INT N = 1; Chris@82: for (i = 0; i < rnk; ++i) N *= local_n[i]; Chris@82: return N; Chris@82: } Chris@82: Chris@82: ptrdiff_t XM(local_size_guru)(int rnk, const XM(ddim) *dims0, Chris@82: ptrdiff_t howmany, MPI_Comm comm, Chris@82: ptrdiff_t *local_n_in, Chris@82: ptrdiff_t *local_start_in, Chris@82: ptrdiff_t *local_n_out, Chris@82: ptrdiff_t *local_start_out, Chris@82: int sign, unsigned flags) Chris@82: { Chris@82: INT N; Chris@82: int my_pe, n_pes, i; Chris@82: dtensor *sz; Chris@82: Chris@82: if (rnk == 0) Chris@82: return howmany; Chris@82: Chris@82: MPI_Comm_rank(comm, &my_pe); Chris@82: MPI_Comm_size(comm, &n_pes); Chris@82: sz = default_sz(rnk, dims0, n_pes, 0); Chris@82: Chris@82: /* Now, we must figure out how much local space the user should Chris@82: allocate (or at least an upper bound). This depends strongly Chris@82: on the exact algorithms we employ...ugh! FIXME: get this info Chris@82: from the solvers somehow? */ Chris@82: N = 1; /* never return zero allocation size */ Chris@82: if (rnk > 1 && XM(is_block1d)(sz, IB) && XM(is_block1d)(sz, OB)) { Chris@82: INT Nafter; Chris@82: ddim odims[2]; Chris@82: Chris@82: /* dft-rank-geq2-transposed */ Chris@82: odims[0] = sz->dims[0]; odims[1] = sz->dims[1]; /* save */ Chris@82: /* we may need extra space for transposed intermediate data */ Chris@82: for (i = 0; i < 2; ++i) Chris@82: if (XM(num_blocks)(sz->dims[i].n, sz->dims[i].b[IB]) == 1 && Chris@82: XM(num_blocks)(sz->dims[i].n, sz->dims[i].b[OB]) == 1) { Chris@82: sz->dims[i].b[IB] Chris@82: = XM(default_block)(sz->dims[i].n, n_pes); Chris@82: sz->dims[1-i].b[IB] = sz->dims[1-i].n; Chris@82: local_size(my_pe, sz, IB, local_n_in, local_start_in); Chris@82: N = X(imax)(N, prod(rnk, local_n_in)); Chris@82: sz->dims[i] = odims[i]; Chris@82: sz->dims[1-i] = odims[1-i]; Chris@82: break; Chris@82: } Chris@82: Chris@82: /* dft-rank-geq2 */ Chris@82: Nafter = howmany; Chris@82: for (i = 1; i < sz->rnk; ++i) Nafter *= sz->dims[i].n; Chris@82: N = X(imax)(N, (sz->dims[0].n Chris@82: * XM(block)(Nafter, XM(default_block)(Nafter, n_pes), Chris@82: my_pe) + howmany - 1) / howmany); Chris@82: Chris@82: /* dft-rank-geq2 with dimensions swapped */ Chris@82: Nafter = howmany * sz->dims[0].n; Chris@82: for (i = 2; i < sz->rnk; ++i) Nafter *= sz->dims[i].n; Chris@82: N = X(imax)(N, (sz->dims[1].n Chris@82: * XM(block)(Nafter, XM(default_block)(Nafter, n_pes), Chris@82: my_pe) + howmany - 1) / howmany); Chris@82: } Chris@82: else if (rnk == 1) { Chris@82: if (howmany >= n_pes && !MPI_FLAGS(flags)) { /* dft-rank1-bigvec */ Chris@82: ptrdiff_t n[2], start[2]; Chris@82: dtensor *sz2 = XM(mkdtensor)(2); Chris@82: sz2->dims[0] = sz->dims[0]; Chris@82: sz2->dims[0].b[IB] = sz->dims[0].n; Chris@82: sz2->dims[1].n = sz2->dims[1].b[OB] = howmany; Chris@82: sz2->dims[1].b[IB] = XM(default_block)(howmany, n_pes); Chris@82: local_size(my_pe, sz2, IB, n, start); Chris@82: XM(dtensor_destroy)(sz2); Chris@82: N = X(imax)(N, (prod(2, n) + howmany - 1) / howmany); Chris@82: } Chris@82: else { /* dft-rank1 */ Chris@82: INT r, m, rblock[2], mblock[2]; Chris@82: Chris@82: /* Since the 1d transforms are so different, we require Chris@82: the user to call local_size_1d for this case. Ugh. */ Chris@82: CK(sign == FFTW_FORWARD || sign == FFTW_BACKWARD); Chris@82: Chris@82: if ((r = XM(choose_radix)(sz->dims[0], n_pes, flags, sign, Chris@82: rblock, mblock))) { Chris@82: m = sz->dims[0].n / r; Chris@82: if (flags & FFTW_MPI_SCRAMBLED_IN) Chris@82: sz->dims[0].b[IB] = rblock[IB] * m; Chris@82: else { /* !SCRAMBLED_IN */ Chris@82: sz->dims[0].b[IB] = r * mblock[IB]; Chris@82: N = X(imax)(N, rblock[IB] * m); Chris@82: } Chris@82: if (flags & FFTW_MPI_SCRAMBLED_OUT) Chris@82: sz->dims[0].b[OB] = r * mblock[OB]; Chris@82: else { /* !SCRAMBLED_OUT */ Chris@82: N = X(imax)(N, r * mblock[OB]); Chris@82: sz->dims[0].b[OB] = rblock[OB] * m; Chris@82: } Chris@82: } Chris@82: } Chris@82: } Chris@82: Chris@82: local_size(my_pe, sz, IB, local_n_in, local_start_in); Chris@82: local_size(my_pe, sz, OB, local_n_out, local_start_out); Chris@82: Chris@82: /* at least, make sure we have enough space to store input & output */ Chris@82: N = X(imax)(N, X(imax)(prod(rnk, local_n_in), prod(rnk, local_n_out))); Chris@82: Chris@82: XM(dtensor_destroy)(sz); Chris@82: return N * howmany; Chris@82: } Chris@82: Chris@82: ptrdiff_t XM(local_size_many_transposed)(int rnk, const ptrdiff_t *n, Chris@82: ptrdiff_t howmany, Chris@82: ptrdiff_t xblock, ptrdiff_t yblock, Chris@82: MPI_Comm comm, Chris@82: ptrdiff_t *local_nx, Chris@82: ptrdiff_t *local_x_start, Chris@82: ptrdiff_t *local_ny, Chris@82: ptrdiff_t *local_y_start) Chris@82: { Chris@82: ptrdiff_t N; Chris@82: XM(ddim) *dims; Chris@82: ptrdiff_t *local; Chris@82: Chris@82: if (rnk == 0) { Chris@82: *local_nx = *local_ny = 1; Chris@82: *local_x_start = *local_y_start = 0; Chris@82: return howmany; Chris@82: } Chris@82: Chris@82: dims = simple_dims(rnk, n); Chris@82: local = (ptrdiff_t *) MALLOC(sizeof(ptrdiff_t) * rnk * 4, TENSORS); Chris@82: Chris@82: /* default 1d block distribution, with transposed output Chris@82: if yblock < n[1] */ Chris@82: dims[0].ib = xblock; Chris@82: if (rnk > 1) { Chris@82: if (yblock < n[1]) Chris@82: dims[1].ob = yblock; Chris@82: else Chris@82: dims[0].ob = xblock; Chris@82: } Chris@82: else Chris@82: dims[0].ob = xblock; /* FIXME: 1d not really supported here Chris@82: since we don't have flags/sign */ Chris@82: Chris@82: N = XM(local_size_guru)(rnk, dims, howmany, comm, Chris@82: local, local + rnk, Chris@82: local + 2*rnk, local + 3*rnk, Chris@82: 0, 0); Chris@82: *local_nx = local[0]; Chris@82: *local_x_start = local[rnk]; Chris@82: if (rnk > 1) { Chris@82: *local_ny = local[2*rnk + 1]; Chris@82: *local_y_start = local[3*rnk + 1]; Chris@82: } Chris@82: else { Chris@82: *local_ny = *local_nx; Chris@82: *local_y_start = *local_x_start; Chris@82: } Chris@82: X(ifree)(local); Chris@82: X(ifree)(dims); Chris@82: return N; Chris@82: } Chris@82: Chris@82: ptrdiff_t XM(local_size_many)(int rnk, const ptrdiff_t *n, Chris@82: ptrdiff_t howmany, Chris@82: ptrdiff_t xblock, Chris@82: MPI_Comm comm, Chris@82: ptrdiff_t *local_nx, Chris@82: ptrdiff_t *local_x_start) Chris@82: { Chris@82: ptrdiff_t local_ny, local_y_start; Chris@82: return XM(local_size_many_transposed)(rnk, n, howmany, Chris@82: xblock, rnk > 1 Chris@82: ? n[1] : FFTW_MPI_DEFAULT_BLOCK, Chris@82: comm, Chris@82: local_nx, local_x_start, Chris@82: &local_ny, &local_y_start); Chris@82: } Chris@82: Chris@82: Chris@82: ptrdiff_t XM(local_size_transposed)(int rnk, const ptrdiff_t *n, Chris@82: MPI_Comm comm, Chris@82: ptrdiff_t *local_nx, Chris@82: ptrdiff_t *local_x_start, Chris@82: ptrdiff_t *local_ny, Chris@82: ptrdiff_t *local_y_start) Chris@82: { Chris@82: return XM(local_size_many_transposed)(rnk, n, 1, Chris@82: FFTW_MPI_DEFAULT_BLOCK, Chris@82: FFTW_MPI_DEFAULT_BLOCK, Chris@82: comm, Chris@82: local_nx, local_x_start, Chris@82: local_ny, local_y_start); Chris@82: } Chris@82: Chris@82: ptrdiff_t XM(local_size)(int rnk, const ptrdiff_t *n, Chris@82: MPI_Comm comm, Chris@82: ptrdiff_t *local_nx, Chris@82: ptrdiff_t *local_x_start) Chris@82: { Chris@82: return XM(local_size_many)(rnk, n, 1, FFTW_MPI_DEFAULT_BLOCK, comm, Chris@82: local_nx, local_x_start); Chris@82: } Chris@82: Chris@82: ptrdiff_t XM(local_size_many_1d)(ptrdiff_t nx, ptrdiff_t howmany, Chris@82: MPI_Comm comm, int sign, unsigned flags, Chris@82: ptrdiff_t *local_nx, ptrdiff_t *local_x_start, Chris@82: ptrdiff_t *local_ny, ptrdiff_t *local_y_start) Chris@82: { Chris@82: XM(ddim) d; Chris@82: d.n = nx; Chris@82: d.ib = d.ob = FFTW_MPI_DEFAULT_BLOCK; Chris@82: return XM(local_size_guru)(1, &d, howmany, comm, Chris@82: local_nx, local_x_start, Chris@82: local_ny, local_y_start, sign, flags); Chris@82: } Chris@82: Chris@82: ptrdiff_t XM(local_size_1d)(ptrdiff_t nx, Chris@82: MPI_Comm comm, int sign, unsigned flags, Chris@82: ptrdiff_t *local_nx, ptrdiff_t *local_x_start, Chris@82: ptrdiff_t *local_ny, ptrdiff_t *local_y_start) Chris@82: { Chris@82: return XM(local_size_many_1d)(nx, 1, comm, sign, flags, Chris@82: local_nx, local_x_start, Chris@82: local_ny, local_y_start); Chris@82: } Chris@82: Chris@82: ptrdiff_t XM(local_size_2d_transposed)(ptrdiff_t nx, ptrdiff_t ny, Chris@82: MPI_Comm comm, Chris@82: ptrdiff_t *local_nx, Chris@82: ptrdiff_t *local_x_start, Chris@82: ptrdiff_t *local_ny, Chris@82: ptrdiff_t *local_y_start) Chris@82: { Chris@82: ptrdiff_t n[2]; Chris@82: n[0] = nx; n[1] = ny; Chris@82: return XM(local_size_transposed)(2, n, comm, Chris@82: local_nx, local_x_start, Chris@82: local_ny, local_y_start); Chris@82: } Chris@82: Chris@82: ptrdiff_t XM(local_size_2d)(ptrdiff_t nx, ptrdiff_t ny, MPI_Comm comm, Chris@82: ptrdiff_t *local_nx, ptrdiff_t *local_x_start) Chris@82: { Chris@82: ptrdiff_t n[2]; Chris@82: n[0] = nx; n[1] = ny; Chris@82: return XM(local_size)(2, n, comm, local_nx, local_x_start); Chris@82: } Chris@82: Chris@82: ptrdiff_t XM(local_size_3d_transposed)(ptrdiff_t nx, ptrdiff_t ny, Chris@82: ptrdiff_t nz, Chris@82: MPI_Comm comm, Chris@82: ptrdiff_t *local_nx, Chris@82: ptrdiff_t *local_x_start, Chris@82: ptrdiff_t *local_ny, Chris@82: ptrdiff_t *local_y_start) Chris@82: { Chris@82: ptrdiff_t n[3]; Chris@82: n[0] = nx; n[1] = ny; n[2] = nz; Chris@82: return XM(local_size_transposed)(3, n, comm, Chris@82: local_nx, local_x_start, Chris@82: local_ny, local_y_start); Chris@82: } Chris@82: Chris@82: ptrdiff_t XM(local_size_3d)(ptrdiff_t nx, ptrdiff_t ny, ptrdiff_t nz, Chris@82: MPI_Comm comm, Chris@82: ptrdiff_t *local_nx, ptrdiff_t *local_x_start) Chris@82: { Chris@82: ptrdiff_t n[3]; Chris@82: n[0] = nx; n[1] = ny; n[2] = nz; Chris@82: return XM(local_size)(3, n, comm, local_nx, local_x_start); Chris@82: } Chris@82: Chris@82: /*************************************************************************/ Chris@82: /* Transpose API */ Chris@82: Chris@82: X(plan) XM(plan_many_transpose)(ptrdiff_t nx, ptrdiff_t ny, Chris@82: ptrdiff_t howmany, Chris@82: ptrdiff_t xblock, ptrdiff_t yblock, Chris@82: R *in, R *out, Chris@82: MPI_Comm comm, unsigned flags) Chris@82: { Chris@82: int n_pes; Chris@82: XM(init)(); Chris@82: Chris@82: if (howmany < 0 || xblock < 0 || yblock < 0 || Chris@82: nx <= 0 || ny <= 0) return 0; Chris@82: Chris@82: MPI_Comm_size(comm, &n_pes); Chris@82: if (!xblock) xblock = XM(default_block)(nx, n_pes); Chris@82: if (!yblock) yblock = XM(default_block)(ny, n_pes); Chris@82: if (n_pes < XM(num_blocks)(nx, xblock) Chris@82: || n_pes < XM(num_blocks)(ny, yblock)) Chris@82: return 0; Chris@82: Chris@82: return Chris@82: X(mkapiplan)(FFTW_FORWARD, flags, Chris@82: XM(mkproblem_transpose)(nx, ny, howmany, Chris@82: in, out, xblock, yblock, Chris@82: comm, MPI_FLAGS(flags))); Chris@82: } Chris@82: Chris@82: X(plan) XM(plan_transpose)(ptrdiff_t nx, ptrdiff_t ny, R *in, R *out, Chris@82: MPI_Comm comm, unsigned flags) Chris@82: Chris@82: { Chris@82: return XM(plan_many_transpose)(nx, ny, 1, Chris@82: FFTW_MPI_DEFAULT_BLOCK, Chris@82: FFTW_MPI_DEFAULT_BLOCK, Chris@82: in, out, comm, flags); Chris@82: } Chris@82: Chris@82: /*************************************************************************/ Chris@82: /* Complex DFT API */ Chris@82: Chris@82: X(plan) XM(plan_guru_dft)(int rnk, const XM(ddim) *dims0, Chris@82: ptrdiff_t howmany, Chris@82: C *in, C *out, Chris@82: MPI_Comm comm, int sign, unsigned flags) Chris@82: { Chris@82: int n_pes, i; Chris@82: dtensor *sz; Chris@82: Chris@82: XM(init)(); Chris@82: Chris@82: if (howmany < 0 || rnk < 1) return 0; Chris@82: for (i = 0; i < rnk; ++i) Chris@82: if (dims0[i].n < 1 || dims0[i].ib < 0 || dims0[i].ob < 0) Chris@82: return 0; Chris@82: Chris@82: MPI_Comm_size(comm, &n_pes); Chris@82: sz = default_sz(rnk, dims0, n_pes, 0); Chris@82: Chris@82: if (XM(num_blocks_total)(sz, IB) > n_pes Chris@82: || XM(num_blocks_total)(sz, OB) > n_pes) { Chris@82: XM(dtensor_destroy)(sz); Chris@82: return 0; Chris@82: } Chris@82: Chris@82: return Chris@82: X(mkapiplan)(sign, flags, Chris@82: XM(mkproblem_dft_d)(sz, howmany, Chris@82: (R *) in, (R *) out, Chris@82: comm, sign, Chris@82: MPI_FLAGS(flags))); Chris@82: } Chris@82: Chris@82: X(plan) XM(plan_many_dft)(int rnk, const ptrdiff_t *n, Chris@82: ptrdiff_t howmany, Chris@82: ptrdiff_t iblock, ptrdiff_t oblock, Chris@82: C *in, C *out, Chris@82: MPI_Comm comm, int sign, unsigned flags) Chris@82: { Chris@82: XM(ddim) *dims = simple_dims(rnk, n); Chris@82: X(plan) pln; Chris@82: Chris@82: if (rnk == 1) { Chris@82: dims[0].ib = iblock; Chris@82: dims[0].ob = oblock; Chris@82: } Chris@82: else if (rnk > 1) { Chris@82: dims[0 != (flags & FFTW_MPI_TRANSPOSED_IN)].ib = iblock; Chris@82: dims[0 != (flags & FFTW_MPI_TRANSPOSED_OUT)].ob = oblock; Chris@82: } Chris@82: Chris@82: pln = XM(plan_guru_dft)(rnk,dims,howmany, in,out, comm, sign, flags); Chris@82: X(ifree)(dims); Chris@82: return pln; Chris@82: } Chris@82: Chris@82: X(plan) XM(plan_dft)(int rnk, const ptrdiff_t *n, C *in, C *out, Chris@82: MPI_Comm comm, int sign, unsigned flags) Chris@82: { Chris@82: return XM(plan_many_dft)(rnk, n, 1, Chris@82: FFTW_MPI_DEFAULT_BLOCK, Chris@82: FFTW_MPI_DEFAULT_BLOCK, Chris@82: in, out, comm, sign, flags); Chris@82: } Chris@82: Chris@82: X(plan) XM(plan_dft_1d)(ptrdiff_t nx, C *in, C *out, Chris@82: MPI_Comm comm, int sign, unsigned flags) Chris@82: { Chris@82: return XM(plan_dft)(1, &nx, in, out, comm, sign, flags); Chris@82: } Chris@82: Chris@82: X(plan) XM(plan_dft_2d)(ptrdiff_t nx, ptrdiff_t ny, C *in, C *out, Chris@82: MPI_Comm comm, int sign, unsigned flags) Chris@82: { Chris@82: ptrdiff_t n[2]; Chris@82: n[0] = nx; n[1] = ny; Chris@82: return XM(plan_dft)(2, n, in, out, comm, sign, flags); Chris@82: } Chris@82: Chris@82: X(plan) XM(plan_dft_3d)(ptrdiff_t nx, ptrdiff_t ny, ptrdiff_t nz, Chris@82: C *in, C *out, Chris@82: MPI_Comm comm, int sign, unsigned flags) Chris@82: { Chris@82: ptrdiff_t n[3]; Chris@82: n[0] = nx; n[1] = ny; n[2] = nz; Chris@82: return XM(plan_dft)(3, n, in, out, comm, sign, flags); Chris@82: } Chris@82: Chris@82: /*************************************************************************/ Chris@82: /* R2R API */ Chris@82: Chris@82: X(plan) XM(plan_guru_r2r)(int rnk, const XM(ddim) *dims0, Chris@82: ptrdiff_t howmany, Chris@82: R *in, R *out, Chris@82: MPI_Comm comm, const X(r2r_kind) *kind, Chris@82: unsigned flags) Chris@82: { Chris@82: int n_pes, i; Chris@82: dtensor *sz; Chris@82: rdft_kind *k; Chris@82: X(plan) pln; Chris@82: Chris@82: XM(init)(); Chris@82: Chris@82: if (howmany < 0 || rnk < 1) return 0; Chris@82: for (i = 0; i < rnk; ++i) Chris@82: if (dims0[i].n < 1 || dims0[i].ib < 0 || dims0[i].ob < 0) Chris@82: return 0; Chris@82: Chris@82: k = X(map_r2r_kind)(rnk, kind); Chris@82: Chris@82: MPI_Comm_size(comm, &n_pes); Chris@82: sz = default_sz(rnk, dims0, n_pes, 0); Chris@82: Chris@82: if (XM(num_blocks_total)(sz, IB) > n_pes Chris@82: || XM(num_blocks_total)(sz, OB) > n_pes) { Chris@82: XM(dtensor_destroy)(sz); Chris@82: return 0; Chris@82: } Chris@82: Chris@82: pln = X(mkapiplan)(0, flags, Chris@82: XM(mkproblem_rdft_d)(sz, howmany, Chris@82: in, out, Chris@82: comm, k, MPI_FLAGS(flags))); Chris@82: X(ifree0)(k); Chris@82: return pln; Chris@82: } Chris@82: Chris@82: X(plan) XM(plan_many_r2r)(int rnk, const ptrdiff_t *n, Chris@82: ptrdiff_t howmany, Chris@82: ptrdiff_t iblock, ptrdiff_t oblock, Chris@82: R *in, R *out, Chris@82: MPI_Comm comm, const X(r2r_kind) *kind, Chris@82: unsigned flags) Chris@82: { Chris@82: XM(ddim) *dims = simple_dims(rnk, n); Chris@82: X(plan) pln; Chris@82: Chris@82: if (rnk == 1) { Chris@82: dims[0].ib = iblock; Chris@82: dims[0].ob = oblock; Chris@82: } Chris@82: else if (rnk > 1) { Chris@82: dims[0 != (flags & FFTW_MPI_TRANSPOSED_IN)].ib = iblock; Chris@82: dims[0 != (flags & FFTW_MPI_TRANSPOSED_OUT)].ob = oblock; Chris@82: } Chris@82: Chris@82: pln = XM(plan_guru_r2r)(rnk,dims,howmany, in,out, comm, kind, flags); Chris@82: X(ifree)(dims); Chris@82: return pln; Chris@82: } Chris@82: Chris@82: X(plan) XM(plan_r2r)(int rnk, const ptrdiff_t *n, R *in, R *out, Chris@82: MPI_Comm comm, Chris@82: const X(r2r_kind) *kind, Chris@82: unsigned flags) Chris@82: { Chris@82: return XM(plan_many_r2r)(rnk, n, 1, Chris@82: FFTW_MPI_DEFAULT_BLOCK, Chris@82: FFTW_MPI_DEFAULT_BLOCK, Chris@82: in, out, comm, kind, flags); Chris@82: } Chris@82: Chris@82: X(plan) XM(plan_r2r_2d)(ptrdiff_t nx, ptrdiff_t ny, R *in, R *out, Chris@82: MPI_Comm comm, Chris@82: X(r2r_kind) kindx, X(r2r_kind) kindy, Chris@82: unsigned flags) Chris@82: { Chris@82: ptrdiff_t n[2]; Chris@82: X(r2r_kind) kind[2]; Chris@82: n[0] = nx; n[1] = ny; Chris@82: kind[0] = kindx; kind[1] = kindy; Chris@82: return XM(plan_r2r)(2, n, in, out, comm, kind, flags); Chris@82: } Chris@82: Chris@82: X(plan) XM(plan_r2r_3d)(ptrdiff_t nx, ptrdiff_t ny, ptrdiff_t nz, Chris@82: R *in, R *out, Chris@82: MPI_Comm comm, Chris@82: X(r2r_kind) kindx, X(r2r_kind) kindy, Chris@82: X(r2r_kind) kindz, Chris@82: unsigned flags) Chris@82: { Chris@82: ptrdiff_t n[3]; Chris@82: X(r2r_kind) kind[3]; Chris@82: n[0] = nx; n[1] = ny; n[2] = nz; Chris@82: kind[0] = kindx; kind[1] = kindy; kind[2] = kindz; Chris@82: return XM(plan_r2r)(3, n, in, out, comm, kind, flags); Chris@82: } Chris@82: Chris@82: /*************************************************************************/ Chris@82: /* R2C/C2R API */ Chris@82: Chris@82: static X(plan) plan_guru_rdft2(int rnk, const XM(ddim) *dims0, Chris@82: ptrdiff_t howmany, Chris@82: R *r, C *c, Chris@82: MPI_Comm comm, rdft_kind kind, unsigned flags) Chris@82: { Chris@82: int n_pes, i; Chris@82: dtensor *sz; Chris@82: R *cr = (R *) c; Chris@82: Chris@82: XM(init)(); Chris@82: Chris@82: if (howmany < 0 || rnk < 2) return 0; Chris@82: for (i = 0; i < rnk; ++i) Chris@82: if (dims0[i].n < 1 || dims0[i].ib < 0 || dims0[i].ob < 0) Chris@82: return 0; Chris@82: Chris@82: MPI_Comm_size(comm, &n_pes); Chris@82: sz = default_sz(rnk, dims0, n_pes, 1); Chris@82: Chris@82: sz->dims[rnk-1].n = dims0[rnk-1].n / 2 + 1; Chris@82: if (XM(num_blocks_total)(sz, IB) > n_pes Chris@82: || XM(num_blocks_total)(sz, OB) > n_pes) { Chris@82: XM(dtensor_destroy)(sz); Chris@82: return 0; Chris@82: } Chris@82: sz->dims[rnk-1].n = dims0[rnk-1].n; Chris@82: Chris@82: if (kind == R2HC) Chris@82: return X(mkapiplan)(0, flags, Chris@82: XM(mkproblem_rdft2_d)(sz, howmany, Chris@82: r, cr, comm, R2HC, Chris@82: MPI_FLAGS(flags))); Chris@82: else Chris@82: return X(mkapiplan)(0, flags, Chris@82: XM(mkproblem_rdft2_d)(sz, howmany, Chris@82: cr, r, comm, HC2R, Chris@82: MPI_FLAGS(flags))); Chris@82: } Chris@82: Chris@82: X(plan) XM(plan_many_dft_r2c)(int rnk, const ptrdiff_t *n, Chris@82: ptrdiff_t howmany, Chris@82: ptrdiff_t iblock, ptrdiff_t oblock, Chris@82: R *in, C *out, Chris@82: MPI_Comm comm, unsigned flags) Chris@82: { Chris@82: XM(ddim) *dims = simple_dims(rnk, n); Chris@82: X(plan) pln; Chris@82: Chris@82: if (rnk == 1) { Chris@82: dims[0].ib = iblock; Chris@82: dims[0].ob = oblock; Chris@82: } Chris@82: else if (rnk > 1) { Chris@82: dims[0 != (flags & FFTW_MPI_TRANSPOSED_IN)].ib = iblock; Chris@82: dims[0 != (flags & FFTW_MPI_TRANSPOSED_OUT)].ob = oblock; Chris@82: } Chris@82: Chris@82: pln = plan_guru_rdft2(rnk,dims,howmany, in,out, comm, R2HC, flags); Chris@82: X(ifree)(dims); Chris@82: return pln; Chris@82: } Chris@82: Chris@82: X(plan) XM(plan_many_dft_c2r)(int rnk, const ptrdiff_t *n, Chris@82: ptrdiff_t howmany, Chris@82: ptrdiff_t iblock, ptrdiff_t oblock, Chris@82: C *in, R *out, Chris@82: MPI_Comm comm, unsigned flags) Chris@82: { Chris@82: XM(ddim) *dims = simple_dims(rnk, n); Chris@82: X(plan) pln; Chris@82: Chris@82: if (rnk == 1) { Chris@82: dims[0].ib = iblock; Chris@82: dims[0].ob = oblock; Chris@82: } Chris@82: else if (rnk > 1) { Chris@82: dims[0 != (flags & FFTW_MPI_TRANSPOSED_IN)].ib = iblock; Chris@82: dims[0 != (flags & FFTW_MPI_TRANSPOSED_OUT)].ob = oblock; Chris@82: } Chris@82: Chris@82: pln = plan_guru_rdft2(rnk,dims,howmany, out,in, comm, HC2R, flags); Chris@82: X(ifree)(dims); Chris@82: return pln; Chris@82: } Chris@82: Chris@82: X(plan) XM(plan_dft_r2c)(int rnk, const ptrdiff_t *n, R *in, C *out, Chris@82: MPI_Comm comm, unsigned flags) Chris@82: { Chris@82: return XM(plan_many_dft_r2c)(rnk, n, 1, Chris@82: FFTW_MPI_DEFAULT_BLOCK, Chris@82: FFTW_MPI_DEFAULT_BLOCK, Chris@82: in, out, comm, flags); Chris@82: } Chris@82: Chris@82: X(plan) XM(plan_dft_r2c_2d)(ptrdiff_t nx, ptrdiff_t ny, R *in, C *out, Chris@82: MPI_Comm comm, unsigned flags) Chris@82: { Chris@82: ptrdiff_t n[2]; Chris@82: n[0] = nx; n[1] = ny; Chris@82: return XM(plan_dft_r2c)(2, n, in, out, comm, flags); Chris@82: } Chris@82: Chris@82: X(plan) XM(plan_dft_r2c_3d)(ptrdiff_t nx, ptrdiff_t ny, ptrdiff_t nz, Chris@82: R *in, C *out, MPI_Comm comm, unsigned flags) Chris@82: { Chris@82: ptrdiff_t n[3]; Chris@82: n[0] = nx; n[1] = ny; n[2] = nz; Chris@82: return XM(plan_dft_r2c)(3, n, in, out, comm, flags); Chris@82: } Chris@82: Chris@82: X(plan) XM(plan_dft_c2r)(int rnk, const ptrdiff_t *n, C *in, R *out, Chris@82: MPI_Comm comm, unsigned flags) Chris@82: { Chris@82: return XM(plan_many_dft_c2r)(rnk, n, 1, Chris@82: FFTW_MPI_DEFAULT_BLOCK, Chris@82: FFTW_MPI_DEFAULT_BLOCK, Chris@82: in, out, comm, flags); Chris@82: } Chris@82: Chris@82: X(plan) XM(plan_dft_c2r_2d)(ptrdiff_t nx, ptrdiff_t ny, C *in, R *out, Chris@82: MPI_Comm comm, unsigned flags) Chris@82: { Chris@82: ptrdiff_t n[2]; Chris@82: n[0] = nx; n[1] = ny; Chris@82: return XM(plan_dft_c2r)(2, n, in, out, comm, flags); Chris@82: } Chris@82: Chris@82: X(plan) XM(plan_dft_c2r_3d)(ptrdiff_t nx, ptrdiff_t ny, ptrdiff_t nz, Chris@82: C *in, R *out, MPI_Comm comm, unsigned flags) Chris@82: { Chris@82: ptrdiff_t n[3]; Chris@82: n[0] = nx; n[1] = ny; n[2] = nz; Chris@82: return XM(plan_dft_c2r)(3, n, in, out, comm, flags); Chris@82: } Chris@82: Chris@82: /*************************************************************************/ Chris@82: /* New-array execute functions */ Chris@82: Chris@82: void XM(execute_dft)(const X(plan) p, C *in, C *out) { Chris@82: /* internally, MPI plans are just rdft plans */ Chris@82: X(execute_r2r)(p, (R*) in, (R*) out); Chris@82: } Chris@82: Chris@82: void XM(execute_dft_r2c)(const X(plan) p, R *in, C *out) { Chris@82: /* internally, MPI plans are just rdft plans */ Chris@82: X(execute_r2r)(p, in, (R*) out); Chris@82: } Chris@82: Chris@82: void XM(execute_dft_c2r)(const X(plan) p, C *in, R *out) { Chris@82: /* internally, MPI plans are just rdft plans */ Chris@82: X(execute_r2r)(p, (R*) in, out); Chris@82: } Chris@82: Chris@82: void XM(execute_r2r)(const X(plan) p, R *in, R *out) { Chris@82: /* internally, MPI plans are just rdft plans */ Chris@82: X(execute_r2r)(p, in, out); Chris@82: }