cannam@95: /*
cannam@95:  * Copyright (c) 2003, 2007-11 Matteo Frigo
cannam@95:  * Copyright (c) 2003, 2007-11 Massachusetts Institute of Technology
cannam@95:  *
cannam@95:  * This program is free software; you can redistribute it and/or modify
cannam@95:  * it under the terms of the GNU General Public License as published by
cannam@95:  * the Free Software Foundation; either version 2 of the License, or
cannam@95:  * (at your option) any later version.
cannam@95:  *
cannam@95:  * This program is distributed in the hope that it will be useful,
cannam@95:  * but WITHOUT ANY WARRANTY; without even the implied warranty of
cannam@95:  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
cannam@95:  * GNU General Public License for more details.
cannam@95:  *
cannam@95:  * You should have received a copy of the GNU General Public License
cannam@95:  * along with this program; if not, write to the Free Software
cannam@95:  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
cannam@95:  *
cannam@95:  */
cannam@95: 
cannam@95: #include "ifftw-mpi.h"
cannam@95: 
cannam@95: /* Return the radix r for a 1d MPI transform of a distributed dimension d,
cannam@95:    with the given flags and transform size.   That is, decomposes d.n
cannam@95:    as r * m, Cooley-Tukey style.  Also computes the block sizes rblock
cannam@95:    and mblock.  Returns 0 if such a decomposition is not feasible.
cannam@95:    This is unfortunately somewhat complicated.
cannam@95: 
cannam@95:    A distributed Cooley-Tukey algorithm works as follows (see dft-rank1.c):
cannam@95: 
cannam@95:    d.n is initially distributed as an m x r array with block size mblock[IB].
cannam@95:    Then it is internally transposed to an r x m array with block size
cannam@95:    rblock[IB].  Then it is internally transposed to m x r again with block
cannam@95:    size mblock[OB].  Finally, it is transposed to r x m with block size
cannam@95:    rblock[IB].
cannam@95: 
cannam@95:    If flags & SCRAMBLED_IN, then the first transpose is skipped (the array
cannam@95:    starts out as r x m).  If flags & SCRAMBLED_OUT, then the last transpose
cannam@95:    is skipped (the array ends up as m x r).  To make sure the forward
cannam@95:    and backward transforms use the same "scrambling" format, we swap r
cannam@95:    and m when sign != FFT_SIGN.
cannam@95: 
cannam@95:    There are some downsides to this, especially in the case where
cannam@95:    either m or r is not divisible by n_pes.  For one thing, it means
cannam@95:    that in general we can't use the same block size for the input and
cannam@95:    output.  For another thing, it means that we can't in general honor
cannam@95:    a user's "requested" block sizes in d.b[].  Therefore, for simplicity,
cannam@95:    we simply ignore d.b[] for now.
cannam@95: */
cannam@95: INT XM(choose_radix)(ddim d, int n_pes, unsigned flags, int sign,
cannam@95: 		     INT rblock[2], INT mblock[2])
cannam@95: {
cannam@95:      INT r, m;
cannam@95: 
cannam@95:      UNUSED(flags); /* we would need this if we paid attention to d.b[*] */
cannam@95: 
cannam@95:      /* If n_pes is a factor of d.n, then choose r to be d.n / n_pes.
cannam@95:         This not only ensures that the input (the m dimension) is
cannam@95:         equally distributed if possible, and at the r dimension is
cannam@95:         maximally equally distributed (if d.n/n_pes >= n_pes), it also
cannam@95:         makes one of the local transpositions in the algorithm
cannam@95:         trivial. */
cannam@95:      if (d.n % n_pes == 0 /* it's good if n_pes divides d.n ...*/
cannam@95: 	 && d.n / n_pes >= n_pes /* .. unless we can't use n_pes processes */)
cannam@95: 	  r = d.n / n_pes;
cannam@95:      else {  /* n_pes does not divide d.n, pick a factor close to sqrt(d.n) */
cannam@95: 	  for (r = X(isqrt)(d.n); d.n % r != 0; ++r)
cannam@95: 	       ;
cannam@95:      }
cannam@95:      if (r == 1 || r == d.n) return 0; /* punt if we can't reduce size */
cannam@95: 
cannam@95:      if (sign != FFT_SIGN) { /* swap {m,r} so that scrambling is reversible */
cannam@95: 	  m = r;
cannam@95: 	  r = d.n / m;
cannam@95:      }
cannam@95:      else
cannam@95: 	  m = d.n / r;
cannam@95: 
cannam@95:      rblock[IB] = rblock[OB] = XM(default_block)(r, n_pes);
cannam@95:      mblock[IB] = mblock[OB] = XM(default_block)(m, n_pes);
cannam@95: 
cannam@95:      return r;
cannam@95: }