cannam@167: /**************************************************************************/
cannam@167: /* NOTE to users: this is the FFTW-MPI self-test and benchmark program.
cannam@167:    It is probably NOT a good place to learn FFTW usage, since it has a
cannam@167:    lot of added complexity in order to exercise and test the full API,
cannam@167:    etcetera.  We suggest reading the manual. */
cannam@167: /**************************************************************************/
cannam@167: 
cannam@167: #include <math.h>
cannam@167: #include <stdio.h>
cannam@167: #include <string.h>
cannam@167: #include "fftw3-mpi.h"
cannam@167: #include "tests/fftw-bench.h"
cannam@167: 
cannam@167: #if defined(BENCHFFT_SINGLE)
cannam@167: #  define BENCH_MPI_TYPE MPI_FLOAT
cannam@167: #elif defined(BENCHFFT_LDOUBLE)
cannam@167: #  define BENCH_MPI_TYPE MPI_LONG_DOUBLE
cannam@167: #elif defined(BENCHFFT_QUAD)
cannam@167: #  error MPI quad-precision type is unknown
cannam@167: #else
cannam@167: #  define BENCH_MPI_TYPE MPI_DOUBLE
cannam@167: #endif
cannam@167: 
cannam@167: #if SIZEOF_PTRDIFF_T == SIZEOF_INT
cannam@167: #  define FFTW_MPI_PTRDIFF_T MPI_INT
cannam@167: #elif SIZEOF_PTRDIFF_T == SIZEOF_LONG
cannam@167: #  define FFTW_MPI_PTRDIFF_T MPI_LONG
cannam@167: #elif SIZEOF_PTRDIFF_T == SIZEOF_LONG_LONG
cannam@167: #  define FFTW_MPI_PTRDIFF_T MPI_LONG_LONG
cannam@167: #else
cannam@167: #  error MPI type for ptrdiff_t is unknown
cannam@167: #  define FFTW_MPI_PTRDIFF_T MPI_LONG
cannam@167: #endif
cannam@167: 
cannam@167: static const char *mkversion(void) { return FFTW(version); }
cannam@167: static const char *mkcc(void) { return FFTW(cc); }
cannam@167: static const char *mkcodelet_optim(void) { return FFTW(codelet_optim); }
cannam@167: static const char *mknproc(void) {
cannam@167:      static char buf[32];
cannam@167:      int ncpus;
cannam@167:      MPI_Comm_size(MPI_COMM_WORLD, &ncpus);
cannam@167: #ifdef HAVE_SNPRINTF
cannam@167:      snprintf(buf, 32, "%d", ncpus);
cannam@167: #else
cannam@167:      sprintf(buf, "%d", ncpus);
cannam@167: #endif
cannam@167:      return buf;
cannam@167: }
cannam@167: 
cannam@167: BEGIN_BENCH_DOC
cannam@167: BENCH_DOC("name", "fftw3_mpi")
cannam@167: BENCH_DOCF("version", mkversion)
cannam@167: BENCH_DOCF("cc", mkcc)
cannam@167: BENCH_DOCF("codelet-optim", mkcodelet_optim)
cannam@167: BENCH_DOCF("nproc", mknproc)
cannam@167: END_BENCH_DOC 
cannam@167: 
cannam@167: static int n_pes = 1, my_pe = 0;
cannam@167: 
cannam@167: /* global variables describing the shape of the data and its distribution */
cannam@167: static int rnk;
cannam@167: static ptrdiff_t vn, iNtot, oNtot;
cannam@167: static ptrdiff_t *local_ni=0, *local_starti=0;
cannam@167: static ptrdiff_t *local_no=0, *local_starto=0;
cannam@167: static ptrdiff_t *all_local_ni=0, *all_local_starti=0; /* n_pes x rnk arrays */
cannam@167: static ptrdiff_t *all_local_no=0, *all_local_starto=0; /* n_pes x rnk arrays */
cannam@167: static ptrdiff_t *istrides = 0, *ostrides = 0;
cannam@167: static ptrdiff_t *total_ni=0, *total_no=0;
cannam@167: static int *isend_cnt = 0, *isend_off = 0; /* for MPI_Scatterv */
cannam@167: static int *orecv_cnt = 0, *orecv_off = 0; /* for MPI_Gatherv */
cannam@167: 
cannam@167: static bench_real *local_in = 0, *local_out = 0;
cannam@167: static bench_real *all_local_in = 0, *all_local_out = 0;
cannam@167: static int all_local_in_alloc = 0, all_local_out_alloc = 0;
cannam@167: static FFTW(plan) plan_scramble_in = 0, plan_unscramble_out = 0;
cannam@167: 
cannam@167: static void alloc_rnk(int rnk_) {
cannam@167:      rnk = rnk_;
cannam@167:      bench_free(local_ni);
cannam@167:      if (rnk == 0)
cannam@167: 	  local_ni = 0;
cannam@167:      else
cannam@167: 	  local_ni = (ptrdiff_t *) bench_malloc(sizeof(ptrdiff_t) * rnk
cannam@167: 						* (8 + n_pes * 4));
cannam@167: 
cannam@167:      local_starti = local_ni + rnk;
cannam@167:      local_no = local_ni + 2 * rnk;
cannam@167:      local_starto = local_ni + 3 * rnk;
cannam@167:      istrides = local_ni + 4 * rnk;
cannam@167:      ostrides = local_ni + 5 * rnk;
cannam@167:      total_ni = local_ni + 6 * rnk;
cannam@167:      total_no = local_ni + 7 * rnk;
cannam@167:      all_local_ni = local_ni + 8 * rnk;
cannam@167:      all_local_starti = local_ni + (8 + n_pes) * rnk;
cannam@167:      all_local_no = local_ni + (8 + 2 * n_pes) * rnk;
cannam@167:      all_local_starto = local_ni + (8 + 3 * n_pes) * rnk;
cannam@167: }
cannam@167: 
cannam@167: static void setup_gather_scatter(void)
cannam@167: {
cannam@167:      int i, j;
cannam@167:      ptrdiff_t off;
cannam@167: 
cannam@167:      MPI_Gather(local_ni, rnk, FFTW_MPI_PTRDIFF_T,
cannam@167: 		all_local_ni, rnk, FFTW_MPI_PTRDIFF_T,
cannam@167: 		0, MPI_COMM_WORLD);
cannam@167:      MPI_Bcast(all_local_ni, rnk*n_pes, FFTW_MPI_PTRDIFF_T, 0, MPI_COMM_WORLD);
cannam@167:      MPI_Gather(local_starti, rnk, FFTW_MPI_PTRDIFF_T,
cannam@167: 		all_local_starti, rnk, FFTW_MPI_PTRDIFF_T,
cannam@167: 		0, MPI_COMM_WORLD);
cannam@167:      MPI_Bcast(all_local_starti, rnk*n_pes, FFTW_MPI_PTRDIFF_T, 0, MPI_COMM_WORLD);
cannam@167: 
cannam@167:      MPI_Gather(local_no, rnk, FFTW_MPI_PTRDIFF_T,
cannam@167: 		all_local_no, rnk, FFTW_MPI_PTRDIFF_T,
cannam@167: 		0, MPI_COMM_WORLD);
cannam@167:      MPI_Bcast(all_local_no, rnk*n_pes, FFTW_MPI_PTRDIFF_T, 0, MPI_COMM_WORLD);
cannam@167:      MPI_Gather(local_starto, rnk, FFTW_MPI_PTRDIFF_T,
cannam@167: 		all_local_starto, rnk, FFTW_MPI_PTRDIFF_T,
cannam@167: 		0, MPI_COMM_WORLD);
cannam@167:      MPI_Bcast(all_local_starto, rnk*n_pes, FFTW_MPI_PTRDIFF_T, 0, MPI_COMM_WORLD);
cannam@167: 
cannam@167:      off = 0;
cannam@167:      for (i = 0; i < n_pes; ++i) {
cannam@167: 	  ptrdiff_t N = vn;
cannam@167: 	  for (j = 0; j < rnk; ++j)
cannam@167: 	       N *= all_local_ni[i * rnk + j];
cannam@167: 	  isend_cnt[i] = N;
cannam@167: 	  isend_off[i] = off;
cannam@167: 	  off += N;
cannam@167:      }
cannam@167:      iNtot = off;
cannam@167:      all_local_in_alloc = 1;
cannam@167: 
cannam@167:      istrides[rnk - 1] = vn;
cannam@167:      for (j = rnk - 2; j >= 0; --j)
cannam@167: 	  istrides[j] = total_ni[j + 1] * istrides[j + 1];
cannam@167: 
cannam@167:      off = 0;
cannam@167:      for (i = 0; i < n_pes; ++i) {
cannam@167: 	  ptrdiff_t N = vn;
cannam@167: 	  for (j = 0; j < rnk; ++j)
cannam@167: 	       N *= all_local_no[i * rnk + j];
cannam@167: 	  orecv_cnt[i] = N;
cannam@167: 	  orecv_off[i] = off;
cannam@167: 	  off += N;
cannam@167:      }
cannam@167:      oNtot = off;
cannam@167:      all_local_out_alloc = 1;
cannam@167: 
cannam@167:      ostrides[rnk - 1] = vn;
cannam@167:      for (j = rnk - 2; j >= 0; --j)
cannam@167: 	  ostrides[j] = total_no[j + 1] * ostrides[j + 1];
cannam@167: }
cannam@167: 
cannam@167: static void copy_block_out(const bench_real *in,
cannam@167: 			   int rnk, ptrdiff_t *n, ptrdiff_t *start, 
cannam@167: 			   ptrdiff_t is, ptrdiff_t *os, ptrdiff_t vn,
cannam@167: 			   bench_real *out)
cannam@167: {
cannam@167:      ptrdiff_t i;
cannam@167:      if (rnk == 0) { 
cannam@167: 	  for (i = 0; i < vn; ++i)
cannam@167: 	       out[i] = in[i];
cannam@167:      }
cannam@167:      else if (rnk == 1) { /* this case is just an optimization */
cannam@167: 	  ptrdiff_t j;
cannam@167: 	  out += start[0] * os[0];
cannam@167: 	  for (j = 0; j < n[0]; ++j) {
cannam@167: 	       for (i = 0; i < vn; ++i)
cannam@167: 		    out[i] = in[i];
cannam@167: 	       in += is;
cannam@167: 	       out += os[0];
cannam@167: 	  }
cannam@167:      }
cannam@167:      else {
cannam@167: 	  /* we should do n[0] for locality, but this way is simpler to code */
cannam@167: 	  for (i = 0; i < n[rnk - 1]; ++i) 
cannam@167: 	       copy_block_out(in + i * is,
cannam@167: 			      rnk - 1, n, start, is * n[rnk - 1], os, vn,
cannam@167: 			      out + (start[rnk - 1] + i) * os[rnk - 1]);
cannam@167:      }
cannam@167: }
cannam@167: 
cannam@167: static void copy_block_in(bench_real *in,
cannam@167: 			  int rnk, ptrdiff_t *n, ptrdiff_t *start, 
cannam@167: 			  ptrdiff_t is, ptrdiff_t *os, ptrdiff_t vn,
cannam@167: 			  const bench_real *out)
cannam@167: {
cannam@167:      ptrdiff_t i;
cannam@167:      if (rnk == 0) { 
cannam@167: 	  for (i = 0; i < vn; ++i)
cannam@167: 	       in[i] = out[i];
cannam@167:      }
cannam@167:      else if (rnk == 1) { /* this case is just an optimization */
cannam@167: 	  ptrdiff_t j;
cannam@167: 	  out += start[0] * os[0];
cannam@167: 	  for (j = 0; j < n[0]; ++j) {
cannam@167: 	       for (i = 0; i < vn; ++i)
cannam@167: 		    in[i] = out[i];
cannam@167: 	       in += is;
cannam@167: 	       out += os[0];
cannam@167: 	  }
cannam@167:      }
cannam@167:      else {
cannam@167: 	  /* we should do n[0] for locality, but this way is simpler to code */
cannam@167: 	  for (i = 0; i < n[rnk - 1]; ++i) 
cannam@167: 	       copy_block_in(in + i * is,
cannam@167: 			     rnk - 1, n, start, is * n[rnk - 1], os, vn,
cannam@167: 			     out + (start[rnk - 1] + i) * os[rnk - 1]);
cannam@167:      }
cannam@167: }
cannam@167: 
cannam@167: static void do_scatter_in(bench_real *in)
cannam@167: {
cannam@167:      bench_real *ali;
cannam@167:      int i;
cannam@167:      if (all_local_in_alloc) {
cannam@167:           bench_free(all_local_in);
cannam@167: 	  all_local_in = (bench_real*) bench_malloc(iNtot*sizeof(bench_real));
cannam@167: 	  all_local_in_alloc = 0;
cannam@167:      }
cannam@167:      ali = all_local_in;
cannam@167:      for (i = 0; i < n_pes; ++i) {
cannam@167: 	  copy_block_in(ali,
cannam@167: 			rnk, all_local_ni + i * rnk, 
cannam@167: 			all_local_starti + i * rnk,
cannam@167: 			vn, istrides, vn,
cannam@167: 			in);
cannam@167: 	  ali += isend_cnt[i];
cannam@167:      }
cannam@167:      MPI_Scatterv(all_local_in, isend_cnt, isend_off, BENCH_MPI_TYPE,
cannam@167: 		  local_in, isend_cnt[my_pe], BENCH_MPI_TYPE,
cannam@167: 		  0, MPI_COMM_WORLD);
cannam@167: }
cannam@167: 
cannam@167: static void do_gather_out(bench_real *out)
cannam@167: {
cannam@167:      bench_real *alo;
cannam@167:      int i;
cannam@167: 
cannam@167:      if (all_local_out_alloc) {
cannam@167:           bench_free(all_local_out);
cannam@167: 	  all_local_out = (bench_real*) bench_malloc(oNtot*sizeof(bench_real));
cannam@167: 	  all_local_out_alloc = 0;
cannam@167:      }
cannam@167:      MPI_Gatherv(local_out, orecv_cnt[my_pe], BENCH_MPI_TYPE,
cannam@167: 		 all_local_out, orecv_cnt, orecv_off, BENCH_MPI_TYPE,
cannam@167: 		 0, MPI_COMM_WORLD);
cannam@167:      MPI_Bcast(all_local_out, oNtot, BENCH_MPI_TYPE, 0, MPI_COMM_WORLD);
cannam@167:      alo = all_local_out;
cannam@167:      for (i = 0; i < n_pes; ++i) {
cannam@167: 	  copy_block_out(alo,
cannam@167: 			 rnk, all_local_no + i * rnk, 
cannam@167: 			 all_local_starto + i * rnk,
cannam@167: 			 vn, ostrides, vn,
cannam@167: 			 out);
cannam@167: 	  alo += orecv_cnt[i];
cannam@167:      }
cannam@167: }
cannam@167: 
cannam@167: static void alloc_local(ptrdiff_t nreal, int inplace)
cannam@167: {
cannam@167:      bench_free(local_in);
cannam@167:      if (local_out != local_in) bench_free(local_out);
cannam@167:      local_in = local_out = 0;
cannam@167:      if (nreal > 0) {
cannam@167: 	  ptrdiff_t i;
cannam@167: 	  local_in = (bench_real*) bench_malloc(nreal * sizeof(bench_real));
cannam@167: 	  if (inplace)
cannam@167: 	       local_out = local_in;
cannam@167: 	  else
cannam@167: 	       local_out = (bench_real*) bench_malloc(nreal * sizeof(bench_real));
cannam@167: 	  for (i = 0; i < nreal; ++i) local_in[i] = local_out[i] = 0.0;
cannam@167:      }
cannam@167: }
cannam@167: 
cannam@167: void after_problem_rcopy_from(bench_problem *p, bench_real *ri)
cannam@167: {
cannam@167:      UNUSED(p);
cannam@167:      do_scatter_in(ri);
cannam@167:      if (plan_scramble_in) FFTW(execute)(plan_scramble_in);
cannam@167: }
cannam@167: 
cannam@167: void after_problem_rcopy_to(bench_problem *p, bench_real *ro)
cannam@167: {
cannam@167:      UNUSED(p);
cannam@167:      if (plan_unscramble_out) FFTW(execute)(plan_unscramble_out);
cannam@167:      do_gather_out(ro);
cannam@167: }
cannam@167: 
cannam@167: void after_problem_ccopy_from(bench_problem *p, bench_real *ri, bench_real *ii)
cannam@167: {
cannam@167:      UNUSED(ii);
cannam@167:      after_problem_rcopy_from(p, ri);
cannam@167: }
cannam@167: 
cannam@167: void after_problem_ccopy_to(bench_problem *p, bench_real *ro, bench_real *io)
cannam@167: {
cannam@167:      UNUSED(io);
cannam@167:      after_problem_rcopy_to(p, ro);
cannam@167: }
cannam@167: 
cannam@167: void after_problem_hccopy_from(bench_problem *p, bench_real *ri, bench_real *ii)
cannam@167: {
cannam@167:      UNUSED(ii);
cannam@167:      after_problem_rcopy_from(p, ri);
cannam@167: }
cannam@167: 
cannam@167: void after_problem_hccopy_to(bench_problem *p, bench_real *ro, bench_real *io)
cannam@167: {
cannam@167:      UNUSED(io);
cannam@167:      after_problem_rcopy_to(p, ro);
cannam@167: }
cannam@167: 
cannam@167: static FFTW(plan) mkplan_transpose_local(ptrdiff_t nx, ptrdiff_t ny, 
cannam@167: 					 ptrdiff_t vn, 
cannam@167: 					 bench_real *in, bench_real *out)
cannam@167: {
cannam@167:      FFTW(iodim64) hdims[3];
cannam@167:      FFTW(r2r_kind) k[3];
cannam@167:      FFTW(plan) pln;
cannam@167: 
cannam@167:      hdims[0].n = nx;
cannam@167:      hdims[0].is = ny * vn;
cannam@167:      hdims[0].os = vn;
cannam@167:      hdims[1].n = ny;
cannam@167:      hdims[1].is = vn;
cannam@167:      hdims[1].os = nx * vn;
cannam@167:      hdims[2].n = vn;
cannam@167:      hdims[2].is = 1;
cannam@167:      hdims[2].os = 1;
cannam@167:      k[0] = k[1] = k[2] = FFTW_R2HC;
cannam@167:      pln = FFTW(plan_guru64_r2r)(0, 0, 3, hdims, in, out, k, FFTW_ESTIMATE);
cannam@167:      BENCH_ASSERT(pln != 0);
cannam@167:      return pln;
cannam@167: }
cannam@167: 
cannam@167: static int tensor_rowmajor_transposedp(bench_tensor *t)
cannam@167: {
cannam@167:      bench_iodim *d;
cannam@167:      int i;
cannam@167: 
cannam@167:      BENCH_ASSERT(BENCH_FINITE_RNK(t->rnk));
cannam@167:      if (t->rnk < 2)
cannam@167: 	  return 0;
cannam@167: 
cannam@167:      d = t->dims;
cannam@167:      if (d[0].is != d[1].is * d[1].n
cannam@167: 	 || d[0].os != d[1].is
cannam@167: 	 || d[1].os != d[0].os * d[0].n)
cannam@167: 	  return 0;
cannam@167:      if (t->rnk > 2 && d[1].is != d[2].is * d[2].n)
cannam@167: 	  return 0;
cannam@167:      for (i = 2; i + 1 < t->rnk; ++i) {
cannam@167:           d = t->dims + i;
cannam@167:           if (d[0].is != d[1].is * d[1].n
cannam@167: 	      || d[0].os != d[1].os * d[1].n)
cannam@167:                return 0;
cannam@167:      }
cannam@167: 
cannam@167:      if (t->rnk > 2 && t->dims[t->rnk-1].is != t->dims[t->rnk-1].os)
cannam@167: 	  return 0;
cannam@167:      return 1;
cannam@167: }
cannam@167: 
cannam@167: static int tensor_contiguousp(bench_tensor *t, int s)
cannam@167: {
cannam@167:      return (t->dims[t->rnk-1].is == s
cannam@167: 	     && ((tensor_rowmajorp(t) && 
cannam@167: 		  t->dims[t->rnk-1].is == t->dims[t->rnk-1].os)
cannam@167: 		 || tensor_rowmajor_transposedp(t)));
cannam@167: }
cannam@167: 
cannam@167: static FFTW(plan) mkplan_complex(bench_problem *p, unsigned flags)
cannam@167: {
cannam@167:      FFTW(plan) pln = 0;
cannam@167:      int i; 
cannam@167:      ptrdiff_t ntot;
cannam@167: 
cannam@167:      vn = p->vecsz->rnk == 1 ? p->vecsz->dims[0].n : 1;
cannam@167: 
cannam@167:      if (p->sz->rnk < 1
cannam@167: 	 || p->split
cannam@167: 	 || !tensor_contiguousp(p->sz, vn)
cannam@167: 	 || tensor_rowmajor_transposedp(p->sz)
cannam@167: 	 || p->vecsz->rnk > 1
cannam@167: 	 || (p->vecsz->rnk == 1 && (p->vecsz->dims[0].is != 1
cannam@167: 				    || p->vecsz->dims[0].os != 1)))
cannam@167: 	  return 0;
cannam@167: 
cannam@167:      alloc_rnk(p->sz->rnk);
cannam@167:      for (i = 0; i < rnk; ++i) {
cannam@167: 	  total_ni[i] = total_no[i] = p->sz->dims[i].n;
cannam@167: 	  local_ni[i] = local_no[i] = total_ni[i];
cannam@167: 	  local_starti[i] = local_starto[i] = 0;
cannam@167:      }
cannam@167:      if (rnk > 1) {
cannam@167: 	  ptrdiff_t n, start, nT, startT;
cannam@167: 	  ntot = FFTW(mpi_local_size_many_transposed)
cannam@167: 	       (p->sz->rnk, total_ni, vn,
cannam@167: 		FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK,
cannam@167: 		MPI_COMM_WORLD,
cannam@167: 		&n, &start, &nT, &startT);
cannam@167: 	  if  (flags & FFTW_MPI_TRANSPOSED_IN) {
cannam@167: 	       local_ni[1] = nT;
cannam@167: 	       local_starti[1] = startT;
cannam@167: 	  }
cannam@167: 	  else {
cannam@167: 	       local_ni[0] = n;
cannam@167: 	       local_starti[0] = start;
cannam@167: 	  }
cannam@167: 	  if  (flags & FFTW_MPI_TRANSPOSED_OUT) {
cannam@167: 	       local_no[1] = nT;
cannam@167: 	       local_starto[1] = startT;
cannam@167: 	  }
cannam@167: 	  else {
cannam@167: 	       local_no[0] = n;
cannam@167: 	       local_starto[0] = start;
cannam@167: 	  }
cannam@167:      }
cannam@167:      else if (rnk == 1) {
cannam@167: 	  ntot = FFTW(mpi_local_size_many_1d)
cannam@167: 	       (total_ni[0], vn, MPI_COMM_WORLD, p->sign, flags,
cannam@167: 		local_ni, local_starti, local_no, local_starto);
cannam@167:      }
cannam@167:      alloc_local(ntot * 2, p->in == p->out);
cannam@167: 
cannam@167:      pln = FFTW(mpi_plan_many_dft)(p->sz->rnk, total_ni, vn, 
cannam@167: 				   FFTW_MPI_DEFAULT_BLOCK,
cannam@167: 				   FFTW_MPI_DEFAULT_BLOCK,
cannam@167: 				   (FFTW(complex) *) local_in, 
cannam@167: 				   (FFTW(complex) *) local_out,
cannam@167: 				   MPI_COMM_WORLD, p->sign, flags);
cannam@167: 
cannam@167:      vn *= 2;
cannam@167: 
cannam@167:      if (rnk > 1) {
cannam@167: 	  ptrdiff_t nrest = 1;
cannam@167: 	  for (i = 2; i < rnk; ++i) nrest *= p->sz->dims[i].n;
cannam@167: 	  if (flags & FFTW_MPI_TRANSPOSED_IN)
cannam@167: 	       plan_scramble_in = mkplan_transpose_local(
cannam@167: 		    p->sz->dims[0].n, local_ni[1], vn * nrest,
cannam@167: 		    local_in, local_in);
cannam@167: 	  if (flags & FFTW_MPI_TRANSPOSED_OUT)
cannam@167: 	       plan_unscramble_out = mkplan_transpose_local(
cannam@167: 		    local_no[1], p->sz->dims[0].n, vn * nrest,
cannam@167: 		    local_out, local_out);
cannam@167:      }
cannam@167:      
cannam@167:      return pln;
cannam@167: }
cannam@167: 
cannam@167: static int tensor_real_contiguousp(bench_tensor *t, int sign, int s)
cannam@167: {
cannam@167:      return (t->dims[t->rnk-1].is == s
cannam@167: 	     && ((tensor_real_rowmajorp(t, sign, 1) && 
cannam@167: 		  t->dims[t->rnk-1].is == t->dims[t->rnk-1].os)));
cannam@167: }
cannam@167: 
cannam@167: static FFTW(plan) mkplan_real(bench_problem *p, unsigned flags)
cannam@167: {
cannam@167:      FFTW(plan) pln = 0;
cannam@167:      int i; 
cannam@167:      ptrdiff_t ntot;
cannam@167: 
cannam@167:      vn = p->vecsz->rnk == 1 ? p->vecsz->dims[0].n : 1;
cannam@167: 
cannam@167:      if (p->sz->rnk < 2
cannam@167: 	 || p->split
cannam@167: 	 || !tensor_real_contiguousp(p->sz, p->sign, vn)
cannam@167: 	 || tensor_rowmajor_transposedp(p->sz)
cannam@167: 	 || p->vecsz->rnk > 1
cannam@167: 	 || (p->vecsz->rnk == 1 && (p->vecsz->dims[0].is != 1
cannam@167: 				    || p->vecsz->dims[0].os != 1)))
cannam@167: 	  return 0;
cannam@167: 
cannam@167:      alloc_rnk(p->sz->rnk);
cannam@167:      for (i = 0; i < rnk; ++i) {
cannam@167: 	  total_ni[i] = total_no[i] = p->sz->dims[i].n;
cannam@167: 	  local_ni[i] = local_no[i] = total_ni[i];
cannam@167: 	  local_starti[i] = local_starto[i] = 0;
cannam@167:      }
cannam@167:      local_ni[rnk-1] = local_no[rnk-1] = total_ni[rnk-1] = total_no[rnk-1] 
cannam@167: 	  = p->sz->dims[rnk-1].n / 2 + 1;
cannam@167:      {
cannam@167: 	  ptrdiff_t n, start, nT, startT;
cannam@167: 	  ntot = FFTW(mpi_local_size_many_transposed)
cannam@167: 	       (p->sz->rnk, total_ni, vn,
cannam@167: 		FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK,
cannam@167: 		MPI_COMM_WORLD,
cannam@167: 		&n, &start, &nT, &startT);
cannam@167: 	  if  (flags & FFTW_MPI_TRANSPOSED_IN) {
cannam@167: 	       local_ni[1] = nT;
cannam@167: 	       local_starti[1] = startT;
cannam@167: 	  }
cannam@167: 	  else {
cannam@167: 	       local_ni[0] = n;
cannam@167: 	       local_starti[0] = start;
cannam@167: 	  }
cannam@167: 	  if  (flags & FFTW_MPI_TRANSPOSED_OUT) {
cannam@167: 	       local_no[1] = nT;
cannam@167: 	       local_starto[1] = startT;
cannam@167: 	  }
cannam@167: 	  else {
cannam@167: 	       local_no[0] = n;
cannam@167: 	       local_starto[0] = start;
cannam@167: 	  }
cannam@167:      }
cannam@167:      alloc_local(ntot * 2, p->in == p->out);
cannam@167: 
cannam@167:      total_ni[rnk - 1] = p->sz->dims[rnk - 1].n;
cannam@167:      if (p->sign < 0)
cannam@167: 	  pln = FFTW(mpi_plan_many_dft_r2c)(p->sz->rnk, total_ni, vn, 
cannam@167: 					    FFTW_MPI_DEFAULT_BLOCK,
cannam@167: 					    FFTW_MPI_DEFAULT_BLOCK,
cannam@167: 					    local_in, 
cannam@167: 					    (FFTW(complex) *) local_out,
cannam@167: 					    MPI_COMM_WORLD, flags);
cannam@167:      else
cannam@167: 	  pln = FFTW(mpi_plan_many_dft_c2r)(p->sz->rnk, total_ni, vn, 
cannam@167: 					    FFTW_MPI_DEFAULT_BLOCK,
cannam@167: 					    FFTW_MPI_DEFAULT_BLOCK,
cannam@167: 					    (FFTW(complex) *) local_in, 
cannam@167: 					    local_out,
cannam@167: 					    MPI_COMM_WORLD, flags);
cannam@167: 
cannam@167:      total_ni[rnk - 1] = p->sz->dims[rnk - 1].n / 2 + 1;
cannam@167:      vn *= 2;
cannam@167: 
cannam@167:      {
cannam@167: 	  ptrdiff_t nrest = 1;
cannam@167: 	  for (i = 2; i < rnk; ++i) nrest *= total_ni[i];
cannam@167: 	  if (flags & FFTW_MPI_TRANSPOSED_IN)
cannam@167: 	       plan_scramble_in = mkplan_transpose_local(
cannam@167: 		    total_ni[0], local_ni[1], vn * nrest,
cannam@167: 		    local_in, local_in);
cannam@167: 	  if (flags & FFTW_MPI_TRANSPOSED_OUT)
cannam@167: 	       plan_unscramble_out = mkplan_transpose_local(
cannam@167: 		    local_no[1], total_ni[0], vn * nrest,
cannam@167: 		    local_out, local_out);
cannam@167:      }
cannam@167:      
cannam@167:      return pln;
cannam@167: }
cannam@167: 
cannam@167: static FFTW(plan) mkplan_transpose(bench_problem *p, unsigned flags)
cannam@167: {
cannam@167:      ptrdiff_t ntot, nx, ny;
cannam@167:      int ix=0, iy=1, i;
cannam@167:      const bench_iodim *d = p->vecsz->dims;
cannam@167:      FFTW(plan) pln;
cannam@167: 
cannam@167:      if (p->vecsz->rnk == 3) {
cannam@167: 	  for (i = 0; i < 3; ++i)
cannam@167: 	       if (d[i].is == 1 && d[i].os == 1) {
cannam@167: 		    vn = d[i].n;
cannam@167: 		    ix = (i + 1) % 3;
cannam@167: 		    iy = (i + 2) % 3;
cannam@167: 		    break;
cannam@167: 	       }
cannam@167: 	  if (i == 3) return 0;
cannam@167:      }
cannam@167:      else {
cannam@167: 	  vn = 1;
cannam@167: 	  ix = 0;
cannam@167: 	  iy = 1;
cannam@167:      }
cannam@167: 
cannam@167:      if (d[ix].is == d[iy].n * vn && d[ix].os == vn
cannam@167: 	 && d[iy].os == d[ix].n * vn && d[iy].is == vn) {
cannam@167: 	  nx = d[ix].n;
cannam@167: 	  ny = d[iy].n;
cannam@167:      }
cannam@167:      else if (d[iy].is == d[ix].n * vn && d[iy].os == vn
cannam@167: 	      && d[ix].os == d[iy].n * vn && d[ix].is == vn) {
cannam@167: 	  nx = d[iy].n;
cannam@167: 	  ny = d[ix].n;
cannam@167:      }
cannam@167:      else
cannam@167: 	  return 0;
cannam@167: 
cannam@167:      alloc_rnk(2);
cannam@167:      ntot = vn * FFTW(mpi_local_size_2d_transposed)(nx, ny, MPI_COMM_WORLD,
cannam@167: 						    &local_ni[0], 
cannam@167: 						    &local_starti[0],
cannam@167: 						    &local_no[0], 
cannam@167: 						    &local_starto[0]);
cannam@167:      local_ni[1] = ny;
cannam@167:      local_starti[1] = 0;
cannam@167:      local_no[1] = nx;
cannam@167:      local_starto[1] = 0;
cannam@167:      total_ni[0] = nx; total_ni[1] = ny;
cannam@167:      total_no[1] = nx; total_no[0] = ny;
cannam@167:      alloc_local(ntot, p->in == p->out);
cannam@167: 
cannam@167:      pln = FFTW(mpi_plan_many_transpose)(nx, ny, vn,
cannam@167: 					 FFTW_MPI_DEFAULT_BLOCK,
cannam@167: 					 FFTW_MPI_DEFAULT_BLOCK,
cannam@167: 					 local_in, local_out,
cannam@167: 					 MPI_COMM_WORLD, flags);
cannam@167:      
cannam@167:      if (flags & FFTW_MPI_TRANSPOSED_IN)
cannam@167: 	  plan_scramble_in = mkplan_transpose_local(local_ni[0], ny, vn,
cannam@167: 						    local_in, local_in);
cannam@167:      if (flags & FFTW_MPI_TRANSPOSED_OUT)
cannam@167: 	  plan_unscramble_out = mkplan_transpose_local
cannam@167: 	       (nx, local_no[0], vn, local_out, local_out);
cannam@167:      
cannam@167: #if 0
cannam@167:      if (pln && vn == 1) {
cannam@167: 	  int i, j;
cannam@167: 	  bench_real *ri = (bench_real *) p->in;
cannam@167: 	  bench_real *ro = (bench_real *) p->out;
cannam@167: 	  if (!ri || !ro) return pln;
cannam@167: 	  setup_gather_scatter();
cannam@167: 	  for (i = 0; i < nx * ny; ++i)
cannam@167: 	       ri[i] = i;
cannam@167: 	  after_problem_rcopy_from(p, ri);
cannam@167: 	  FFTW(execute)(pln);
cannam@167: 	  after_problem_rcopy_to(p, ro);
cannam@167: 	  if (my_pe == 0) {
cannam@167: 	       for (i = 0; i < nx; ++i) {
cannam@167: 		    for (j = 0; j < ny; ++j)
cannam@167: 			 printf("  %3g", ro[j * nx + i]);
cannam@167: 		    printf("\n");
cannam@167: 	       }
cannam@167: 	  }
cannam@167:      }
cannam@167: #endif
cannam@167: 
cannam@167:      return pln;
cannam@167: }
cannam@167: 
cannam@167: static FFTW(plan) mkplan_r2r(bench_problem *p, unsigned flags)
cannam@167: {
cannam@167:      FFTW(plan) pln = 0;
cannam@167:      int i; 
cannam@167:      ptrdiff_t ntot;
cannam@167:      FFTW(r2r_kind) *k;
cannam@167: 
cannam@167:      if ((p->sz->rnk == 0 || (p->sz->rnk == 1 && p->sz->dims[0].n == 1))
cannam@167: 	 && p->vecsz->rnk >= 2 && p->vecsz->rnk <= 3)
cannam@167: 	  return mkplan_transpose(p, flags);
cannam@167: 
cannam@167:      vn = p->vecsz->rnk == 1 ? p->vecsz->dims[0].n : 1;
cannam@167: 
cannam@167:      if (p->sz->rnk < 1
cannam@167: 	 || p->split
cannam@167: 	 || !tensor_contiguousp(p->sz, vn)
cannam@167: 	 || tensor_rowmajor_transposedp(p->sz)
cannam@167: 	 || p->vecsz->rnk > 1
cannam@167: 	 || (p->vecsz->rnk == 1 && (p->vecsz->dims[0].is != 1
cannam@167: 				    || p->vecsz->dims[0].os != 1)))
cannam@167: 	  return 0;
cannam@167: 
cannam@167:      alloc_rnk(p->sz->rnk);
cannam@167:      for (i = 0; i < rnk; ++i) {
cannam@167: 	  total_ni[i] = total_no[i] = p->sz->dims[i].n;
cannam@167: 	  local_ni[i] = local_no[i] = total_ni[i];
cannam@167: 	  local_starti[i] = local_starto[i] = 0;
cannam@167:      }
cannam@167:      if (rnk > 1) {
cannam@167: 	  ptrdiff_t n, start, nT, startT;
cannam@167: 	  ntot = FFTW(mpi_local_size_many_transposed)
cannam@167: 	       (p->sz->rnk, total_ni, vn,
cannam@167: 		FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK,
cannam@167: 		MPI_COMM_WORLD,
cannam@167: 		&n, &start, &nT, &startT);
cannam@167: 	  if  (flags & FFTW_MPI_TRANSPOSED_IN) {
cannam@167: 	       local_ni[1] = nT;
cannam@167: 	       local_starti[1] = startT;
cannam@167: 	  }
cannam@167: 	  else {
cannam@167: 	       local_ni[0] = n;
cannam@167: 	       local_starti[0] = start;
cannam@167: 	  }
cannam@167: 	  if  (flags & FFTW_MPI_TRANSPOSED_OUT) {
cannam@167: 	       local_no[1] = nT;
cannam@167: 	       local_starto[1] = startT;
cannam@167: 	  }
cannam@167: 	  else {
cannam@167: 	       local_no[0] = n;
cannam@167: 	       local_starto[0] = start;
cannam@167: 	  }
cannam@167:      }
cannam@167:      else if (rnk == 1) {
cannam@167: 	  ntot = FFTW(mpi_local_size_many_1d)
cannam@167: 	       (total_ni[0], vn, MPI_COMM_WORLD, p->sign, flags,
cannam@167: 		local_ni, local_starti, local_no, local_starto);
cannam@167:      }
cannam@167:      alloc_local(ntot, p->in == p->out);
cannam@167: 
cannam@167:      k = (FFTW(r2r_kind) *) bench_malloc(sizeof(FFTW(r2r_kind)) * p->sz->rnk);
cannam@167:      for (i = 0; i < p->sz->rnk; ++i)
cannam@167: 	  switch (p->k[i]) {
cannam@167: 	      case R2R_R2HC: k[i] = FFTW_R2HC; break;
cannam@167: 	      case R2R_HC2R: k[i] = FFTW_HC2R; break;
cannam@167: 	      case R2R_DHT: k[i] = FFTW_DHT; break;
cannam@167: 	      case R2R_REDFT00: k[i] = FFTW_REDFT00; break;
cannam@167: 	      case R2R_REDFT01: k[i] = FFTW_REDFT01; break;
cannam@167: 	      case R2R_REDFT10: k[i] = FFTW_REDFT10; break;
cannam@167: 	      case R2R_REDFT11: k[i] = FFTW_REDFT11; break;
cannam@167: 	      case R2R_RODFT00: k[i] = FFTW_RODFT00; break;
cannam@167: 	      case R2R_RODFT01: k[i] = FFTW_RODFT01; break;
cannam@167: 	      case R2R_RODFT10: k[i] = FFTW_RODFT10; break;
cannam@167: 	      case R2R_RODFT11: k[i] = FFTW_RODFT11; break;
cannam@167: 	      default: BENCH_ASSERT(0);
cannam@167: 	  }
cannam@167: 
cannam@167:      pln = FFTW(mpi_plan_many_r2r)(p->sz->rnk, total_ni, vn, 
cannam@167: 				   FFTW_MPI_DEFAULT_BLOCK,
cannam@167: 				   FFTW_MPI_DEFAULT_BLOCK,
cannam@167: 				   local_in, local_out,
cannam@167: 				   MPI_COMM_WORLD, k, flags);
cannam@167:      bench_free(k);
cannam@167: 
cannam@167:      if (rnk > 1) {
cannam@167: 	  ptrdiff_t nrest = 1;
cannam@167: 	  for (i = 2; i < rnk; ++i) nrest *= p->sz->dims[i].n;
cannam@167: 	  if (flags & FFTW_MPI_TRANSPOSED_IN)
cannam@167: 	       plan_scramble_in = mkplan_transpose_local(
cannam@167: 		    p->sz->dims[0].n, local_ni[1], vn * nrest,
cannam@167: 		    local_in, local_in);
cannam@167: 	  if (flags & FFTW_MPI_TRANSPOSED_OUT)
cannam@167: 	       plan_unscramble_out = mkplan_transpose_local(
cannam@167: 		    local_no[1], p->sz->dims[0].n, vn * nrest,
cannam@167: 		    local_out, local_out);
cannam@167:      }
cannam@167:      
cannam@167:      return pln;
cannam@167: }
cannam@167: 
cannam@167: FFTW(plan) mkplan(bench_problem *p, unsigned flags)
cannam@167: {
cannam@167:      FFTW(plan) pln = 0;
cannam@167:      FFTW(destroy_plan)(plan_scramble_in); plan_scramble_in = 0;
cannam@167:      FFTW(destroy_plan)(plan_unscramble_out); plan_unscramble_out = 0;
cannam@167:      if (p->scrambled_in) {
cannam@167: 	  if (p->sz->rnk == 1 && p->sz->dims[0].n != 1) 
cannam@167: 	       flags |= FFTW_MPI_SCRAMBLED_IN;
cannam@167: 	  else
cannam@167: 	       flags |= FFTW_MPI_TRANSPOSED_IN;
cannam@167:      }
cannam@167:      if (p->scrambled_out) {
cannam@167: 	  if (p->sz->rnk == 1 && p->sz->dims[0].n != 1) 
cannam@167: 	       flags |= FFTW_MPI_SCRAMBLED_OUT;
cannam@167: 	  else
cannam@167: 	       flags |= FFTW_MPI_TRANSPOSED_OUT;
cannam@167:      }
cannam@167:      switch (p->kind) {
cannam@167:          case PROBLEM_COMPLEX: 
cannam@167: 	      pln =mkplan_complex(p, flags);
cannam@167: 	      break;
cannam@167:          case PROBLEM_REAL: 
cannam@167: 	      pln = mkplan_real(p, flags);
cannam@167: 	      break;
cannam@167:          case PROBLEM_R2R:
cannam@167: 	      pln = mkplan_r2r(p, flags);
cannam@167: 	      break;
cannam@167:          default: BENCH_ASSERT(0);
cannam@167:      }
cannam@167:      if (pln) setup_gather_scatter();
cannam@167:      return pln;
cannam@167: }
cannam@167: 
cannam@167: void main_init(int *argc, char ***argv)
cannam@167: {
cannam@167: #ifdef HAVE_SMP
cannam@167: # if MPI_VERSION >= 2 /* for MPI_Init_thread */
cannam@167:      int provided;
cannam@167:      MPI_Init_thread(argc, argv, MPI_THREAD_FUNNELED, &provided);
cannam@167:      threads_ok = provided >= MPI_THREAD_FUNNELED;
cannam@167: # else
cannam@167:      MPI_Init(argc, argv);
cannam@167:      threads_ok = 0;
cannam@167: # endif
cannam@167: #else
cannam@167:      MPI_Init(argc, argv);
cannam@167: #endif
cannam@167:      MPI_Comm_rank(MPI_COMM_WORLD, &my_pe);
cannam@167:      MPI_Comm_size(MPI_COMM_WORLD, &n_pes);
cannam@167:      if (my_pe != 0) verbose = -999;
cannam@167:      no_speed_allocation = 1; /* so we can benchmark transforms > memory */
cannam@167:      always_pad_real = 1; /* out-of-place real transforms are padded */
cannam@167:      isend_cnt = (int *) bench_malloc(sizeof(int) * n_pes);
cannam@167:      isend_off = (int *) bench_malloc(sizeof(int) * n_pes);
cannam@167:      orecv_cnt = (int *) bench_malloc(sizeof(int) * n_pes);
cannam@167:      orecv_off = (int *) bench_malloc(sizeof(int) * n_pes);
cannam@167: 
cannam@167:      /* init_threads must be called before any other FFTW function,
cannam@167: 	including mpi_init, because it has to register the threads hooks
cannam@167: 	before the planner is initalized */
cannam@167: #ifdef HAVE_SMP
cannam@167:      if (threads_ok) { BENCH_ASSERT(FFTW(init_threads)()); }
cannam@167: #endif
cannam@167:      FFTW(mpi_init)();
cannam@167: }
cannam@167: 
cannam@167: void initial_cleanup(void)
cannam@167: {
cannam@167:      alloc_rnk(0);
cannam@167:      alloc_local(0, 0);
cannam@167:      bench_free(all_local_in); all_local_in = 0;
cannam@167:      bench_free(all_local_out); all_local_out = 0;
cannam@167:      bench_free(isend_off); isend_off = 0;
cannam@167:      bench_free(isend_cnt); isend_cnt = 0;
cannam@167:      bench_free(orecv_off); orecv_off = 0;
cannam@167:      bench_free(orecv_cnt); orecv_cnt = 0;
cannam@167:      FFTW(destroy_plan)(plan_scramble_in); plan_scramble_in = 0;
cannam@167:      FFTW(destroy_plan)(plan_unscramble_out); plan_unscramble_out = 0;
cannam@167: }
cannam@167: 
cannam@167: void final_cleanup(void)
cannam@167: {
cannam@167:      MPI_Finalize();
cannam@167: }
cannam@167: 
cannam@167: void bench_exit(int status)
cannam@167: {
cannam@167:      MPI_Abort(MPI_COMM_WORLD, status);
cannam@167: }
cannam@167: 
cannam@167: double bench_cost_postprocess(double cost)
cannam@167: {
cannam@167:      double cost_max;
cannam@167:      MPI_Allreduce(&cost, &cost_max, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
cannam@167:      return cost_max;
cannam@167: }
cannam@167: 
cannam@167: 
cannam@167: int import_wisdom(FILE *f)
cannam@167: {
cannam@167:      int success = 1, sall;
cannam@167:      if (my_pe == 0) success = FFTW(import_wisdom_from_file)(f);
cannam@167:      FFTW(mpi_broadcast_wisdom)(MPI_COMM_WORLD);
cannam@167:      MPI_Allreduce(&success, &sall, 1, MPI_INT, MPI_LAND, MPI_COMM_WORLD);
cannam@167:      return sall;
cannam@167: }
cannam@167: 
cannam@167: void export_wisdom(FILE *f)
cannam@167: {
cannam@167:      FFTW(mpi_gather_wisdom)(MPI_COMM_WORLD);
cannam@167:      if (my_pe == 0) FFTW(export_wisdom_to_file)(f);
cannam@167: }