src/gromacs/fft/fft5d.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2009,2010,2012,2013,2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #include "config.h"
  36
  37 #include <algorithm>
  38
  39 #include <stdio.h>
  40 #include <stdlib.h>
  41 #include <string.h>
  42
  43 #ifdef NOGMX
  44 #define GMX_PARALLEL_ENV_INITIALIZED 1
  45 #else
  46 #ifdef GMX_MPI
  47 #define GMX_PARALLEL_ENV_INITIALIZED 1
  48 #else
  49 #define GMX_PARALLEL_ENV_INITIALIZED 0
  50 #endif
  51 #endif
  52
  53 #include "gromacs/utility/gmxmpi.h"
  54
  55 #ifdef GMX_OPENMP
  56 /* TODO: Do we still need this? Are we still planning ot use fftw + OpenMP? */
  57 #define FFT5D_THREADS
  58 /* requires fftw compiled with openmp */
  59 /* #define FFT5D_FFTW_THREADS (now set by cmake) */
  60 #endif
  61
  62 #include "fft5d.h"
  63 #include <float.h>
  64 #include <math.h>
  65 #include <assert.h>
  66 #include "gromacs/utility/smalloc.h"
  67
  68 #ifndef __FLT_EPSILON__
  69 #define __FLT_EPSILON__ FLT_EPSILON
  70 #define __DBL_EPSILON__ DBL_EPSILON
  71 #endif
  72
  73 #ifdef NOGMX
  74 FILE* debug = 0;
  75 #endif
  76
  77 #include "gromacs/utility/fatalerror.h"
  78
  79
  80 #ifdef GMX_FFT_FFTW3
  81 #include "thread_mpi/mutex.h"
  82 #include "gromacs/utility/exceptions.h"
  83 /* none of the fftw3 calls, except execute(), are thread-safe, so
  84    we need to serialize them with this mutex. */
  85 static tMPI::mutex big_fftw_mutex;
  86 #define FFTW_LOCK try { big_fftw_mutex.lock(); } GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
  87 #define FFTW_UNLOCK try { big_fftw_mutex.unlock(); } GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
  88 #endif /* GMX_FFT_FFTW3 */
  89
  90 /* largest factor smaller than sqrt */
  91 static int lfactor(int z)
  92 {
  93     int i;
  94     for (i = static_cast<int>(sqrt(static_cast<double>(z)));; i--)
  95     {
  96         if (z%i == 0)
  97         {
  98             return i;
  99         }
 100     }
 101     return 1;
 102 }
 103
 104 /* largest factor */
 105 static int l2factor(int z)
 106 {
 107     int i;
 108     if (z == 1)
 109     {
 110         return 1;
 111     }
 112     for (i = z/2;; i--)
 113     {
 114         if (z%i == 0)
 115         {
 116             return i;
 117         }
 118     }
 119     return 1;
 120 }
 121
 122 /* largest prime factor: WARNING: slow recursion, only use for small numbers */
 123 static int lpfactor(int z)
 124 {
 125     int f = l2factor(z);
 126     if (f == 1)
 127     {
 128         return z;
 129     }
 130     return std::max(lpfactor(f), lpfactor(z/f));
 131 }
 132
 133 #ifndef GMX_MPI
 134 #ifdef HAVE_GETTIMEOFDAY
 135 #include <sys/time.h>
 136 double MPI_Wtime()
 137 {
 138     struct timeval tv;
 139     gettimeofday(&tv, 0);
 140     return tv.tv_sec+tv.tv_usec*1e-6;
 141 }
 142 #else
 143 double MPI_Wtime()
 144 {
 145     return 0.0;
 146 }
 147 #endif
 148 #endif
 149
 150 static int vmax(int* a, int s)
 151 {
 152     int i, max = 0;
 153     for (i = 0; i < s; i++)
 154     {
 155         if (a[i] > max)
 156         {
 157             max = a[i];
 158         }
 159     }
 160     return max;
 161 }
 162
 163
 164 /* NxMxK the size of the data
 165  * comm communicator to use for fft5d
 166  * P0 number of processor in 1st axes (can be null for automatic)
 167  * lin is allocated by fft5d because size of array is only known after planning phase
 168  * rlout2 is only used as intermediate buffer - only returned after allocation to reuse for back transform - should not be used by caller
 169  */
 170 fft5d_plan fft5d_plan_3d(int NG, int MG, int KG, MPI_Comm comm[2], int flags, t_complex** rlin, t_complex** rlout, t_complex** rlout2, t_complex** rlout3, int nthreads)
 171 {
 172
 173     int        P[2], bMaster, prank[2], i, t;
 174     int        rNG, rMG, rKG;
 175     int       *N0 = 0, *N1 = 0, *M0 = 0, *M1 = 0, *K0 = 0, *K1 = 0, *oN0 = 0, *oN1 = 0, *oM0 = 0, *oM1 = 0, *oK0 = 0, *oK1 = 0;
 176     int        N[3], M[3], K[3], pN[3], pM[3], pK[3], oM[3], oK[3], *iNin[3] = {0}, *oNin[3] = {0}, *iNout[3] = {0}, *oNout[3] = {0};
 177     int        C[3], rC[3], nP[2];
 178     int        lsize;
 179     t_complex *lin = 0, *lout = 0, *lout2 = 0, *lout3 = 0;
 180     fft5d_plan plan;
 181     int        s;
 182
 183     /* comm, prank and P are in the order of the decomposition (plan->cart is in the order of transposes) */
 184 #ifdef GMX_MPI
 185     if (GMX_PARALLEL_ENV_INITIALIZED && comm[0] != MPI_COMM_NULL)
 186     {
 187         MPI_Comm_size(comm[0], &P[0]);
 188         MPI_Comm_rank(comm[0], &prank[0]);
 189     }
 190     else
 191 #endif
 192     {
 193         P[0]     = 1;
 194         prank[0] = 0;
 195     }
 196 #ifdef GMX_MPI
 197     if (GMX_PARALLEL_ENV_INITIALIZED && comm[1] != MPI_COMM_NULL)
 198     {
 199         MPI_Comm_size(comm[1], &P[1]);
 200         MPI_Comm_rank(comm[1], &prank[1]);
 201     }
 202     else
 203 #endif
 204     {
 205         P[1]     = 1;
 206         prank[1] = 0;
 207     }
 208
 209     bMaster = (prank[0] == 0 && prank[1] == 0);
 210
 211
 212     if (debug)
 213     {
 214         fprintf(debug, "FFT5D: Using %dx%d rank grid, rank %d,%d\n",
 215                 P[0], P[1], prank[0], prank[1]);
 216     }
 217
 218     if (bMaster)
 219     {
 220         if (debug)
 221         {
 222             fprintf(debug, "FFT5D: N: %d, M: %d, K: %d, P: %dx%d, real2complex: %d, backward: %d, order yz: %d, debug %d\n",
 223                     NG, MG, KG, P[0], P[1], (flags&FFT5D_REALCOMPLEX) > 0, (flags&FFT5D_BACKWARD) > 0, (flags&FFT5D_ORDER_YZ) > 0, (flags&FFT5D_DEBUG) > 0);
 224         }
 225         /* The check below is not correct, one prime factor 11 or 13 is ok.
 226            if (fft5d_fmax(fft5d_fmax(lpfactor(NG),lpfactor(MG)),lpfactor(KG))>7) {
 227             printf("WARNING: FFT very slow with prime factors larger 7\n");
 228             printf("Change FFT size or in case you cannot change it look at\n");
 229             printf("http://www.fftw.org/fftw3_doc/Generating-your-own-code.html\n");
 230            }
 231          */
 232     }
 233
 234     if (NG == 0 || MG == 0 || KG == 0)
 235     {
 236         if (bMaster)
 237         {
 238             printf("FFT5D: FATAL: Datasize cannot be zero in any dimension\n");
 239         }
 240         return 0;
 241     }
 242
 243     rNG = NG; rMG = MG; rKG = KG;
 244
 245     if (flags&FFT5D_REALCOMPLEX)
 246     {
 247         if (!(flags&FFT5D_BACKWARD))
 248         {
 249             NG = NG/2+1;
 250         }
 251         else
 252         {
 253             if (!(flags&FFT5D_ORDER_YZ))
 254             {
 255                 MG = MG/2+1;
 256             }
 257             else
 258             {
 259                 KG = KG/2+1;
 260             }
 261         }
 262     }
 263
 264
 265     /*for transpose we need to know the size for each processor not only our own size*/
 266
 267     N0  = (int*)malloc(P[0]*sizeof(int)); N1 = (int*)malloc(P[1]*sizeof(int));
 268     M0  = (int*)malloc(P[0]*sizeof(int)); M1 = (int*)malloc(P[1]*sizeof(int));
 269     K0  = (int*)malloc(P[0]*sizeof(int)); K1 = (int*)malloc(P[1]*sizeof(int));
 270     oN0 = (int*)malloc(P[0]*sizeof(int)); oN1 = (int*)malloc(P[1]*sizeof(int));
 271     oM0 = (int*)malloc(P[0]*sizeof(int)); oM1 = (int*)malloc(P[1]*sizeof(int));
 272     oK0 = (int*)malloc(P[0]*sizeof(int)); oK1 = (int*)malloc(P[1]*sizeof(int));
 273
 274     for (i = 0; i < P[0]; i++)
 275     {
 276         #define EVENDIST
 277         #ifndef EVENDIST
 278         oN0[i] = i*ceil((double)NG/P[0]);
 279         oM0[i] = i*ceil((double)MG/P[0]);
 280         oK0[i] = i*ceil((double)KG/P[0]);
 281         #else
 282         oN0[i] = (NG*i)/P[0];
 283         oM0[i] = (MG*i)/P[0];
 284         oK0[i] = (KG*i)/P[0];
 285         #endif
 286     }
 287     for (i = 0; i < P[1]; i++)
 288     {
 289         #ifndef EVENDIST
 290         oN1[i] = i*ceil((double)NG/P[1]);
 291         oM1[i] = i*ceil((double)MG/P[1]);
 292         oK1[i] = i*ceil((double)KG/P[1]);
 293         #else
 294         oN1[i] = (NG*i)/P[1];
 295         oM1[i] = (MG*i)/P[1];
 296         oK1[i] = (KG*i)/P[1];
 297         #endif
 298     }
 299     for (i = 0; i < P[0]-1; i++)
 300     {
 301         N0[i] = oN0[i+1]-oN0[i];
 302         M0[i] = oM0[i+1]-oM0[i];
 303         K0[i] = oK0[i+1]-oK0[i];
 304     }
 305     N0[P[0]-1] = NG-oN0[P[0]-1];
 306     M0[P[0]-1] = MG-oM0[P[0]-1];
 307     K0[P[0]-1] = KG-oK0[P[0]-1];
 308     for (i = 0; i < P[1]-1; i++)
 309     {
 310         N1[i] = oN1[i+1]-oN1[i];
 311         M1[i] = oM1[i+1]-oM1[i];
 312         K1[i] = oK1[i+1]-oK1[i];
 313     }
 314     N1[P[1]-1] = NG-oN1[P[1]-1];
 315     M1[P[1]-1] = MG-oM1[P[1]-1];
 316     K1[P[1]-1] = KG-oK1[P[1]-1];
 317
 318     /* for step 1-3 the local N,M,K sizes of the transposed system
 319        C: contiguous dimension, and nP: number of processor in subcommunicator
 320        for that step */
 321
 322
 323     pM[0] = M0[prank[0]];
 324     oM[0] = oM0[prank[0]];
 325     pK[0] = K1[prank[1]];
 326     oK[0] = oK1[prank[1]];
 327     C[0]  = NG;
 328     rC[0] = rNG;
 329     if (!(flags&FFT5D_ORDER_YZ))
 330     {
 331         N[0]     = vmax(N1, P[1]);
 332         M[0]     = M0[prank[0]];
 333         K[0]     = vmax(K1, P[1]);
 334         pN[0]    = N1[prank[1]];
 335         iNout[0] = N1;
 336         oNout[0] = oN1;
 337         nP[0]    = P[1];
 338         C[1]     = KG;
 339         rC[1]    = rKG;
 340         N[1]     = vmax(K0, P[0]);
 341         pN[1]    = K0[prank[0]];
 342         iNin[1]  = K1;
 343         oNin[1]  = oK1;
 344         iNout[1] = K0;
 345         oNout[1] = oK0;
 346         M[1]     = vmax(M0, P[0]);
 347         pM[1]    = M0[prank[0]];
 348         oM[1]    = oM0[prank[0]];
 349         K[1]     = N1[prank[1]];
 350         pK[1]    = N1[prank[1]];
 351         oK[1]    = oN1[prank[1]];
 352         nP[1]    = P[0];
 353         C[2]     = MG;
 354         rC[2]    = rMG;
 355         iNin[2]  = M0;
 356         oNin[2]  = oM0;
 357         M[2]     = vmax(K0, P[0]);
 358         pM[2]    = K0[prank[0]];
 359         oM[2]    = oK0[prank[0]];
 360         K[2]     = vmax(N1, P[1]);
 361         pK[2]    = N1[prank[1]];
 362         oK[2]    = oN1[prank[1]];
 363         free(N0); free(oN0); /*these are not used for this order*/
 364         free(M1); free(oM1); /*the rest is freed in destroy*/
 365     }
 366     else
 367     {
 368         N[0]     = vmax(N0, P[0]);
 369         M[0]     = vmax(M0, P[0]);
 370         K[0]     = K1[prank[1]];
 371         pN[0]    = N0[prank[0]];
 372         iNout[0] = N0;
 373         oNout[0] = oN0;
 374         nP[0]    = P[0];
 375         C[1]     = MG;
 376         rC[1]    = rMG;
 377         N[1]     = vmax(M1, P[1]);
 378         pN[1]    = M1[prank[1]];
 379         iNin[1]  = M0;
 380         oNin[1]  = oM0;
 381         iNout[1] = M1;
 382         oNout[1] = oM1;
 383         M[1]     = N0[prank[0]];
 384         pM[1]    = N0[prank[0]];
 385         oM[1]    = oN0[prank[0]];
 386         K[1]     = vmax(K1, P[1]);
 387         pK[1]    = K1[prank[1]];
 388         oK[1]    = oK1[prank[1]];
 389         nP[1]    = P[1];
 390         C[2]     = KG;
 391         rC[2]    = rKG;
 392         iNin[2]  = K1;
 393         oNin[2]  = oK1;
 394         M[2]     = vmax(N0, P[0]);
 395         pM[2]    = N0[prank[0]];
 396         oM[2]    = oN0[prank[0]];
 397         K[2]     = vmax(M1, P[1]);
 398         pK[2]    = M1[prank[1]];
 399         oK[2]    = oM1[prank[1]];
 400         free(N1); free(oN1); /*these are not used for this order*/
 401         free(K0); free(oK0); /*the rest is freed in destroy*/
 402     }
 403     N[2] = pN[2] = -1;       /*not used*/
 404
 405     /*
 406        Difference between x-y-z regarding 2d decomposition is whether they are
 407        distributed along axis 1, 2 or both
 408      */
 409
 410     /* int lsize = fmax(N[0]*M[0]*K[0]*nP[0],N[1]*M[1]*K[1]*nP[1]); */
 411     lsize = std::max(N[0]*M[0]*K[0]*nP[0], std::max(N[1]*M[1]*K[1]*nP[1], C[2]*M[2]*K[2]));
 412     /* int lsize = fmax(C[0]*M[0]*K[0],fmax(C[1]*M[1]*K[1],C[2]*M[2]*K[2])); */
 413     if (!(flags&FFT5D_NOMALLOC))
 414     {
 415         snew_aligned(lin, lsize, 32);
 416         snew_aligned(lout, lsize, 32);
 417         if (nthreads > 1)
 418         {
 419             /* We need extra transpose buffers to avoid OpenMP barriers */
 420             snew_aligned(lout2, lsize, 32);
 421             snew_aligned(lout3, lsize, 32);
 422         }
 423         else
 424         {
 425             /* We can reuse the buffers to avoid cache misses */
 426             lout2 = lin;
 427             lout3 = lout;
 428         }
 429     }
 430     else
 431     {
 432         lin  = *rlin;
 433         lout = *rlout;
 434         if (nthreads > 1)
 435         {
 436             lout2 = *rlout2;
 437             lout3 = *rlout3;
 438         }
 439         else
 440         {
 441             lout2 = lin;
 442             lout3 = lout;
 443         }
 444     }
 445
 446     plan = (fft5d_plan)calloc(1, sizeof(struct fft5d_plan_t));
 447
 448
 449     if (debug)
 450     {
 451         fprintf(debug, "Running on %d threads\n", nthreads);
 452     }
 453
 454 #ifdef GMX_FFT_FFTW3                                                            /*if not FFTW - then we don't do a 3d plan but instead use only 1D plans */
 455     /* It is possible to use the 3d plan with OMP threads - but in that case it is not allowed to be called from
 456      * within a parallel region. For now deactivated. If it should be supported it has to made sure that
 457      * that the execute of the 3d plan is in a master/serial block (since it contains it own parallel region)
 458      * and that the 3d plan is faster than the 1d plan.
 459      */
 460     if ((!(flags&FFT5D_INPLACE)) && (!(P[0] > 1 || P[1] > 1)) && nthreads == 1) /*don't do 3d plan in parallel or if in_place requested */
 461     {
 462         int fftwflags = FFTW_DESTROY_INPUT;
 463         FFTW(iodim) dims[3];
 464         int inNG = NG, outMG = MG, outKG = KG;
 465
 466         FFTW_LOCK;
 467         if (!(flags&FFT5D_NOMEASURE))
 468         {
 469             fftwflags |= FFTW_MEASURE;
 470         }
 471         if (flags&FFT5D_REALCOMPLEX)
 472         {
 473             if (!(flags&FFT5D_BACKWARD))        /*input pointer is not complex*/
 474             {
 475                 inNG *= 2;
 476             }
 477             else                                /*output pointer is not complex*/
 478             {
 479                 if (!(flags&FFT5D_ORDER_YZ))
 480                 {
 481                     outMG *= 2;
 482                 }
 483                 else
 484                 {
 485                     outKG *= 2;
 486                 }
 487             }
 488         }
 489
 490         if (!(flags&FFT5D_BACKWARD))
 491         {
 492             dims[0].n  = KG;
 493             dims[1].n  = MG;
 494             dims[2].n  = rNG;
 495
 496             dims[0].is = inNG*MG;         /*N M K*/
 497             dims[1].is = inNG;
 498             dims[2].is = 1;
 499             if (!(flags&FFT5D_ORDER_YZ))
 500             {
 501                 dims[0].os = MG;           /*M K N*/
 502                 dims[1].os = 1;
 503                 dims[2].os = MG*KG;
 504             }
 505             else
 506             {
 507                 dims[0].os = 1;           /*K N M*/
 508                 dims[1].os = KG*NG;
 509                 dims[2].os = KG;
 510             }
 511         }
 512         else
 513         {
 514             if (!(flags&FFT5D_ORDER_YZ))
 515             {
 516                 dims[0].n  = NG;
 517                 dims[1].n  = KG;
 518                 dims[2].n  = rMG;
 519
 520                 dims[0].is = 1;
 521                 dims[1].is = NG*MG;
 522                 dims[2].is = NG;
 523
 524                 dims[0].os = outMG*KG;
 525                 dims[1].os = outMG;
 526                 dims[2].os = 1;
 527             }
 528             else
 529             {
 530                 dims[0].n  = MG;
 531                 dims[1].n  = NG;
 532                 dims[2].n  = rKG;
 533
 534                 dims[0].is = NG;
 535                 dims[1].is = 1;
 536                 dims[2].is = NG*MG;
 537
 538                 dims[0].os = outKG*NG;
 539                 dims[1].os = outKG;
 540                 dims[2].os = 1;
 541             }
 542         }
 543 #ifdef FFT5D_THREADS
 544 #ifdef FFT5D_FFTW_THREADS
 545         FFTW(plan_with_nthreads)(nthreads);
 546 #endif
 547 #endif
 548         if ((flags&FFT5D_REALCOMPLEX) && !(flags&FFT5D_BACKWARD))
 549         {
 550             plan->p3d = FFTW(plan_guru_dft_r2c)(/*rank*/ 3, dims,
 551                                                          /*howmany*/ 0, /*howmany_dims*/ 0,
 552                                                          (real*)lin, (FFTW(complex) *) lout,
 553                                                          /*flags*/ fftwflags);
 554         }
 555         else if ((flags&FFT5D_REALCOMPLEX) && (flags&FFT5D_BACKWARD))
 556         {
 557             plan->p3d = FFTW(plan_guru_dft_c2r)(/*rank*/ 3, dims,
 558                                                          /*howmany*/ 0, /*howmany_dims*/ 0,
 559                                                          (FFTW(complex) *) lin, (real*)lout,
 560                                                          /*flags*/ fftwflags);
 561         }
 562         else
 563         {
 564             plan->p3d = FFTW(plan_guru_dft)(/*rank*/ 3, dims,
 565                                                      /*howmany*/ 0, /*howmany_dims*/ 0,
 566                                                      (FFTW(complex) *) lin, (FFTW(complex) *) lout,
 567                                                      /*sign*/ (flags&FFT5D_BACKWARD) ? 1 : -1, /*flags*/ fftwflags);
 568         }
 569 #ifdef FFT5D_THREADS
 570 #ifdef FFT5D_FFTW_THREADS
 571         FFTW(plan_with_nthreads)(1);
 572 #endif
 573 #endif
 574         FFTW_UNLOCK;
 575     }
 576     if (!plan->p3d) /* for decomposition and if 3d plan did not work */
 577     {
 578 #endif              /* GMX_FFT_FFTW3 */
 579     for (s = 0; s < 3; s++)
 580     {
 581         if (debug)
 582         {
 583             fprintf(debug, "FFT5D: Plan s %d rC %d M %d pK %d C %d lsize %d\n",
 584                     s, rC[s], M[s], pK[s], C[s], lsize);
 585         }
 586         plan->p1d[s] = (gmx_fft_t*)malloc(sizeof(gmx_fft_t)*nthreads);
 587
 588         /* Make sure that the init routines are only called by one thread at a time and in order
 589            (later is only important to not confuse valgrind)
 590          */
 591 #pragma omp parallel for num_threads(nthreads) schedule(static) ordered
 592         for (t = 0; t < nthreads; t++)
 593         {
 594 #pragma omp ordered
 595             {
 596                 int tsize = ((t+1)*pM[s]*pK[s]/nthreads)-(t*pM[s]*pK[s]/nthreads);
 597
 598                 if ((flags&FFT5D_REALCOMPLEX) && ((!(flags&FFT5D_BACKWARD) && s == 0) || ((flags&FFT5D_BACKWARD) && s == 2)))
 599                 {
 600                     gmx_fft_init_many_1d_real( &plan->p1d[s][t], rC[s], tsize, (flags&FFT5D_NOMEASURE) ? GMX_FFT_FLAG_CONSERVATIVE : 0 );
 601                 }
 602                 else
 603                 {
 604                     gmx_fft_init_many_1d     ( &plan->p1d[s][t],  C[s], tsize, (flags&FFT5D_NOMEASURE) ? GMX_FFT_FLAG_CONSERVATIVE : 0 );
 605                 }
 606             }
 607         }
 608     }
 609
 610 #ifdef GMX_FFT_FFTW3
 611 }
 612 #endif
 613     if ((flags&FFT5D_ORDER_YZ))   /*plan->cart is in the order of transposes */
 614     {
 615         plan->cart[0] = comm[0]; plan->cart[1] = comm[1];
 616     }
 617     else
 618     {
 619         plan->cart[1] = comm[0]; plan->cart[0] = comm[1];
 620     }
 621 #ifdef FFT5D_MPI_TRANSPOSE
 622     FFTW_LOCK;
 623     for (s = 0; s < 2; s++)
 624     {
 625         if ((s == 0 && !(flags&FFT5D_ORDER_YZ)) || (s == 1 && (flags&FFT5D_ORDER_YZ)))
 626         {
 627             plan->mpip[s] = FFTW(mpi_plan_many_transpose)(nP[s], nP[s], N[s]*K[s]*pM[s]*2, 1, 1, (real*)lout2, (real*)lout3, plan->cart[s], FFTW_PATIENT);
 628         }
 629         else
 630         {
 631             plan->mpip[s] = FFTW(mpi_plan_many_transpose)(nP[s], nP[s], N[s]*pK[s]*M[s]*2, 1, 1, (real*)lout2, (real*)lout3, plan->cart[s], FFTW_PATIENT);
 632         }
 633     }
 634     FFTW_UNLOCK;
 635 #endif
 636
 637
 638     plan->lin   = lin;
 639     plan->lout  = lout;
 640     plan->lout2 = lout2;
 641     plan->lout3 = lout3;
 642
 643     plan->NG = NG; plan->MG = MG; plan->KG = KG;
 644
 645     for (s = 0; s < 3; s++)
 646     {
 647         plan->N[s]    = N[s]; plan->M[s] = M[s]; plan->K[s] = K[s]; plan->pN[s] = pN[s]; plan->pM[s] = pM[s]; plan->pK[s] = pK[s];
 648         plan->oM[s]   = oM[s]; plan->oK[s] = oK[s];
 649         plan->C[s]    = C[s]; plan->rC[s] = rC[s];
 650         plan->iNin[s] = iNin[s]; plan->oNin[s] = oNin[s]; plan->iNout[s] = iNout[s]; plan->oNout[s] = oNout[s];
 651     }
 652     for (s = 0; s < 2; s++)
 653     {
 654         plan->P[s] = nP[s]; plan->coor[s] = prank[s];
 655     }
 656
 657 /*    plan->fftorder=fftorder;
 658     plan->direction=direction;
 659     plan->realcomplex=realcomplex;
 660  */
 661     plan->flags    = flags;
 662     plan->nthreads = nthreads;
 663     *rlin          = lin;
 664     *rlout         = lout;
 665     *rlout2        = lout2;
 666     *rlout3        = lout3;
 667     return plan;
 668 }
 669
 670
 671 enum order {
 672     XYZ,
 673     XZY,
 674     YXZ,
 675     YZX,
 676     ZXY,
 677     ZYX
 678 };
 679
 680
 681
 682 /*here x,y,z and N,M,K is in rotated coordinate system!!
 683    x (and N) is mayor (consecutive) dimension, y (M) middle and z (K) major
 684    maxN,maxM,maxK is max size of local data
 685    pN, pM, pK is local size specific to current processor (only different to max if not divisible)
 686    NG, MG, KG is size of global data*/
 687 static void splitaxes(t_complex* lout, const t_complex* lin,
 688                       int maxN, int maxM, int maxK, int pM,
 689                       int P, int NG, int *N, int* oN, int starty, int startz, int endy, int endz)
 690 {
 691     int x, y, z, i;
 692     int in_i, out_i, in_z, out_z, in_y, out_y;
 693     int s_y, e_y;
 694
 695     for (z = startz; z < endz+1; z++) /*3. z l*/
 696     {
 697         if (z == startz)
 698         {
 699             s_y = starty;
 700         }
 701         else
 702         {
 703             s_y = 0;
 704         }
 705         if (z == endz)
 706         {
 707             e_y = endy;
 708         }
 709         else
 710         {
 711             e_y = pM;
 712         }
 713         out_z  = z*maxN*maxM;
 714         in_z   = z*NG*pM;
 715
 716         for (i = 0; i < P; i++) /*index cube along long axis*/
 717         {
 718             out_i  = out_z  + i*maxN*maxM*maxK;
 719             in_i   = in_z + oN[i];
 720             for (y = s_y; y < e_y; y++)   /*2. y k*/
 721             {
 722                 out_y  = out_i  + y*maxN;
 723                 in_y   = in_i + y*NG;
 724                 for (x = 0; x < N[i]; x++)       /*1. x j*/
 725                 {
 726                     lout[out_y+x] = lin[in_y+x]; /*in=z*NG*pM+oN[i]+y*NG+x*/
 727                     /*after split important that each processor chunk i has size maxN*maxM*maxK and thus being the same size*/
 728                     /*before split data contiguos - thus if different processor get different amount oN is different*/
 729                 }
 730             }
 731         }
 732     }
 733 }
 734
 735 /*make axis contiguous again (after AllToAll) and also do local transpose*/
 736 /*transpose mayor and major dimension
 737    variables see above
 738    the major, middle, minor order is only correct for x,y,z (N,M,K) for the input
 739    N,M,K local dimensions
 740    KG global size*/
 741 static void joinAxesTrans13(t_complex* lout, const t_complex* lin,
 742                             int maxN, int maxM, int maxK, int pM,
 743                             int P, int KG, int* K, int* oK, int starty, int startx, int endy, int endx)
 744 {
 745     int i, x, y, z;
 746     int out_i, in_i, out_x, in_x, out_z, in_z;
 747     int s_y, e_y;
 748
 749     for (x = startx; x < endx+1; x++) /*1.j*/
 750     {
 751         if (x == startx)
 752         {
 753             s_y = starty;
 754         }
 755         else
 756         {
 757             s_y = 0;
 758         }
 759         if (x == endx)
 760         {
 761             e_y = endy;
 762         }
 763         else
 764         {
 765             e_y = pM;
 766         }
 767
 768         out_x  = x*KG*pM;
 769         in_x   = x;
 770
 771         for (i = 0; i < P; i++) /*index cube along long axis*/
 772         {
 773             out_i  = out_x  + oK[i];
 774             in_i   = in_x + i*maxM*maxN*maxK;
 775             for (z = 0; z < K[i]; z++) /*3.l*/
 776             {
 777                 out_z  = out_i  + z;
 778                 in_z   = in_i + z*maxM*maxN;
 779                 for (y = s_y; y < e_y; y++)              /*2.k*/
 780                 {
 781                     lout[out_z+y*KG] = lin[in_z+y*maxN]; /*out=x*KG*pM+oK[i]+z+y*KG*/
 782                 }
 783             }
 784         }
 785     }
 786 }
 787
 788 /*make axis contiguous again (after AllToAll) and also do local transpose
 789    tranpose mayor and middle dimension
 790    variables see above
 791    the minor, middle, major order is only correct for x,y,z (N,M,K) for the input
 792    N,M,K local size
 793    MG, global size*/
 794 static void joinAxesTrans12(t_complex* lout, const t_complex* lin, int maxN, int maxM, int maxK, int pN,
 795                             int P, int MG, int* M, int* oM, int startx, int startz, int endx, int endz)
 796 {
 797     int i, z, y, x;
 798     int out_i, in_i, out_z, in_z, out_x, in_x;
 799     int s_x, e_x;
 800
 801     for (z = startz; z < endz+1; z++)
 802     {
 803         if (z == startz)
 804         {
 805             s_x = startx;
 806         }
 807         else
 808         {
 809             s_x = 0;
 810         }
 811         if (z == endz)
 812         {
 813             e_x = endx;
 814         }
 815         else
 816         {
 817             e_x = pN;
 818         }
 819         out_z  = z*MG*pN;
 820         in_z   = z*maxM*maxN;
 821
 822         for (i = 0; i < P; i++) /*index cube along long axis*/
 823         {
 824             out_i  = out_z  + oM[i];
 825             in_i   = in_z + i*maxM*maxN*maxK;
 826             for (x = s_x; x < e_x; x++)
 827             {
 828                 out_x  = out_i  + x*MG;
 829                 in_x   = in_i + x;
 830                 for (y = 0; y < M[i]; y++)
 831                 {
 832                     lout[out_x+y] = lin[in_x+y*maxN]; /*out=z*MG*pN+oM[i]+x*MG+y*/
 833                 }
 834             }
 835         }
 836     }
 837 }
 838
 839
 840 static void rotate_offsets(int x[])
 841 {
 842     int t = x[0];
 843 /*    x[0]=x[2];
 844     x[2]=x[1];
 845     x[1]=t;*/
 846     x[0] = x[1];
 847     x[1] = x[2];
 848     x[2] = t;
 849 }
 850
 851 /*compute the offset to compare or print transposed local data in original input coordinates
 852    xs matrix dimension size, xl dimension length, xc decomposition offset
 853    s: step in computation = number of transposes*/
 854 static void compute_offsets(fft5d_plan plan, int xs[], int xl[], int xc[], int NG[], int s)
 855 {
 856 /*    int direction = plan->direction;
 857     int fftorder = plan->fftorder;*/
 858
 859     int  o = 0;
 860     int  pos[3], i;
 861     int *pM = plan->pM, *pK = plan->pK, *oM = plan->oM, *oK = plan->oK,
 862     *C      = plan->C, *rC = plan->rC;
 863
 864     NG[0] = plan->NG; NG[1] = plan->MG; NG[2] = plan->KG;
 865
 866     if (!(plan->flags&FFT5D_ORDER_YZ))
 867     {
 868         switch (s)
 869         {
 870             case 0: o = XYZ; break;
 871             case 1: o = ZYX; break;
 872             case 2: o = YZX; break;
 873             default: assert(0);
 874         }
 875     }
 876     else
 877     {
 878         switch (s)
 879         {
 880             case 0: o = XYZ; break;
 881             case 1: o = YXZ; break;
 882             case 2: o = ZXY; break;
 883             default: assert(0);
 884         }
 885     }
 886
 887     switch (o)
 888     {
 889         case XYZ: pos[0] = 1; pos[1] = 2; pos[2] = 3; break;
 890         case XZY: pos[0] = 1; pos[1] = 3; pos[2] = 2; break;
 891         case YXZ: pos[0] = 2; pos[1] = 1; pos[2] = 3; break;
 892         case YZX: pos[0] = 3; pos[1] = 1; pos[2] = 2; break;
 893         case ZXY: pos[0] = 2; pos[1] = 3; pos[2] = 1; break;
 894         case ZYX: pos[0] = 3; pos[1] = 2; pos[2] = 1; break;
 895     }
 896     /*if (debug) printf("pos: %d %d %d\n",pos[0],pos[1],pos[2]);*/
 897
 898     /*xs, xl give dimension size and data length in local transposed coordinate system
 899        for 0(/1/2): x(/y/z) in original coordinate system*/
 900     for (i = 0; i < 3; i++)
 901     {
 902         switch (pos[i])
 903         {
 904             case 1: xs[i] = 1;         xc[i] = 0;     xl[i] = C[s]; break;
 905             case 2: xs[i] = C[s];      xc[i] = oM[s]; xl[i] = pM[s]; break;
 906             case 3: xs[i] = C[s]*pM[s]; xc[i] = oK[s]; xl[i] = pK[s]; break;
 907         }
 908     }
 909     /*input order is different for test program to match FFTW order
 910        (important for complex to real)*/
 911     if (plan->flags&FFT5D_BACKWARD)
 912     {
 913         rotate_offsets(xs);
 914         rotate_offsets(xl);
 915         rotate_offsets(xc);
 916         rotate_offsets(NG);
 917         if (plan->flags&FFT5D_ORDER_YZ)
 918         {
 919             rotate_offsets(xs);
 920             rotate_offsets(xl);
 921             rotate_offsets(xc);
 922             rotate_offsets(NG);
 923         }
 924     }
 925     if ((plan->flags&FFT5D_REALCOMPLEX) && ((!(plan->flags&FFT5D_BACKWARD) && s == 0) || ((plan->flags&FFT5D_BACKWARD) && s == 2)))
 926     {
 927         xl[0] = rC[s];
 928     }
 929 }
 930
 931 static void print_localdata(const t_complex* lin, const char* txt, int s, fft5d_plan plan)
 932 {
 933     int  x, y, z, l;
 934     int *coor = plan->coor;
 935     int  xs[3], xl[3], xc[3], NG[3];
 936     int  ll = (plan->flags&FFT5D_REALCOMPLEX) ? 1 : 2;
 937     compute_offsets(plan, xs, xl, xc, NG, s);
 938     fprintf(debug, txt, coor[0], coor[1], s);
 939     /*printf("xs: %d %d %d, xl: %d %d %d\n",xs[0],xs[1],xs[2],xl[0],xl[1],xl[2]);*/
 940     for (z = 0; z < xl[2]; z++)
 941     {
 942         for (y = 0; y < xl[1]; y++)
 943         {
 944             fprintf(debug, "%d %d: ", coor[0], coor[1]);
 945             for (x = 0; x < xl[0]; x++)
 946             {
 947                 for (l = 0; l < ll; l++)
 948                 {
 949                     fprintf(debug, "%f ", ((real*)lin)[(z*xs[2]+y*xs[1])*2+(x*xs[0])*ll+l]);
 950                 }
 951                 fprintf(debug, ",");
 952             }
 953             fprintf(debug, "\n");
 954         }
 955     }
 956 }
 957
 958 void fft5d_execute(fft5d_plan plan, int thread, fft5d_time times)
 959 {
 960     t_complex  *lin   = plan->lin;
 961     t_complex  *lout  = plan->lout;
 962     t_complex  *lout2 = plan->lout2;
 963     t_complex  *lout3 = plan->lout3;
 964     t_complex  *fftout, *joinin;
 965
 966     gmx_fft_t **p1d = plan->p1d;
 967 #ifdef FFT5D_MPI_TRANSPOSE
 968     FFTW(plan) *mpip = plan->mpip;
 969 #endif
 970 #ifdef GMX_MPI
 971     MPI_Comm *cart = plan->cart;
 972 #endif
 973 #ifdef NOGMX
 974     double time_fft = 0, time_local = 0, time_mpi[2] = {0}, time = 0;
 975 #endif
 976     int   *N = plan->N, *M = plan->M, *K = plan->K, *pN = plan->pN, *pM = plan->pM, *pK = plan->pK,
 977     *C       = plan->C, *P = plan->P, **iNin = plan->iNin, **oNin = plan->oNin, **iNout = plan->iNout, **oNout = plan->oNout;
 978     int    s = 0, tstart, tend, bParallelDim;
 979
 980
 981 #ifdef GMX_FFT_FFTW3
 982     if (plan->p3d)
 983     {
 984         if (thread == 0)
 985         {
 986 #ifdef NOGMX
 987             if (times != 0)
 988             {
 989                 time = MPI_Wtime();
 990             }
 991 #endif
 992             FFTW(execute)(plan->p3d);
 993 #ifdef NOGMX
 994             if (times != 0)
 995             {
 996                 times->fft += MPI_Wtime()-time;
 997             }
 998 #endif
 999         }
1000         return;
1001     }
1002 #endif
1003
1004     s = 0;
1005
1006     /*lin: x,y,z*/
1007     if (plan->flags&FFT5D_DEBUG && thread == 0)
1008     {
1009         print_localdata(lin, "%d %d: copy in lin\n", s, plan);
1010     }
1011
1012     for (s = 0; s < 2; s++)  /*loop over first two FFT steps (corner rotations)*/
1013
1014     {
1015 #ifdef GMX_MPI
1016         if (GMX_PARALLEL_ENV_INITIALIZED && cart[s] != MPI_COMM_NULL && P[s] > 1)
1017         {
1018             bParallelDim = 1;
1019         }
1020         else
1021 #endif
1022         {
1023             bParallelDim = 0;
1024         }
1025
1026         /* ---------- START FFT ------------ */
1027 #ifdef NOGMX
1028         if (times != 0 && thread == 0)
1029         {
1030             time = MPI_Wtime();
1031         }
1032 #endif
1033
1034         if (bParallelDim || plan->nthreads == 1)
1035         {
1036             fftout = lout;
1037         }
1038         else
1039         {
1040             if (s == 0)
1041             {
1042                 fftout = lout3;
1043             }
1044             else
1045             {
1046                 fftout = lout2;
1047             }
1048         }
1049
1050         tstart = (thread*pM[s]*pK[s]/plan->nthreads)*C[s];
1051         if ((plan->flags&FFT5D_REALCOMPLEX) && !(plan->flags&FFT5D_BACKWARD) && s == 0)
1052         {
1053             gmx_fft_many_1d_real(p1d[s][thread], (plan->flags&FFT5D_BACKWARD) ? GMX_FFT_COMPLEX_TO_REAL : GMX_FFT_REAL_TO_COMPLEX, lin+tstart, fftout+tstart);
1054         }
1055         else
1056         {
1057             gmx_fft_many_1d(     p1d[s][thread], (plan->flags&FFT5D_BACKWARD) ? GMX_FFT_BACKWARD : GMX_FFT_FORWARD,               lin+tstart, fftout+tstart);
1058
1059         }
1060
1061 #ifdef NOGMX
1062         if (times != NULL && thread == 0)
1063         {
1064             time_fft += MPI_Wtime()-time;
1065         }
1066 #endif
1067         if (plan->flags&FFT5D_DEBUG && thread == 0)
1068         {
1069             print_localdata(lout, "%d %d: FFT %d\n", s, plan);
1070         }
1071         /* ---------- END FFT ------------ */
1072
1073         /* ---------- START SPLIT + TRANSPOSE------------ (if parallel in in this dimension)*/
1074         if (bParallelDim)
1075         {
1076 #ifdef NOGMX
1077             if (times != NULL && thread == 0)
1078             {
1079                 time = MPI_Wtime();
1080             }
1081 #endif
1082             /*prepare for A
1083                llToAll
1084                1. (most outer) axes (x) is split into P[s] parts of size N[s]
1085                for sending*/
1086             if (pM[s] > 0)
1087             {
1088                 tend    = ((thread+1)*pM[s]*pK[s]/plan->nthreads);
1089                 tstart /= C[s];
1090                 splitaxes(lout2, lout, N[s], M[s], K[s], pM[s], P[s], C[s], iNout[s], oNout[s], tstart%pM[s], tstart/pM[s], tend%pM[s], tend/pM[s]);
1091             }
1092 #pragma omp barrier /*barrier required before AllToAll (all input has to be their) - before timing to make timing more acurate*/
1093 #ifdef NOGMX
1094             if (times != NULL && thread == 0)
1095             {
1096                 time_local += MPI_Wtime()-time;
1097             }
1098 #endif
1099
1100             /* ---------- END SPLIT , START TRANSPOSE------------ */
1101
1102             if (thread == 0)
1103             {
1104 #ifdef NOGMX
1105                 if (times != 0)
1106                 {
1107                     time = MPI_Wtime();
1108                 }
1109 #else
1110                 wallcycle_start(times, ewcPME_FFTCOMM);
1111 #endif
1112 #ifdef FFT5D_MPI_TRANSPOSE
1113                 FFTW(execute)(mpip[s]);
1114 #else
1115 #ifdef GMX_MPI
1116                 if ((s == 0 && !(plan->flags&FFT5D_ORDER_YZ)) || (s == 1 && (plan->flags&FFT5D_ORDER_YZ)))
1117                 {
1118                     MPI_Alltoall((real *)lout2, N[s]*pM[s]*K[s]*sizeof(t_complex)/sizeof(real), GMX_MPI_REAL, (real *)lout3, N[s]*pM[s]*K[s]*sizeof(t_complex)/sizeof(real), GMX_MPI_REAL, cart[s]);
1119                 }
1120                 else
1121                 {
1122                     MPI_Alltoall((real *)lout2, N[s]*M[s]*pK[s]*sizeof(t_complex)/sizeof(real), GMX_MPI_REAL, (real *)lout3, N[s]*M[s]*pK[s]*sizeof(t_complex)/sizeof(real), GMX_MPI_REAL, cart[s]);
1123                 }
1124 #else
1125                 gmx_incons("fft5d MPI call without MPI configuration");
1126 #endif /*GMX_MPI*/
1127 #endif /*FFT5D_MPI_TRANSPOSE*/
1128 #ifdef NOGMX
1129                 if (times != 0)
1130                 {
1131                     time_mpi[s] = MPI_Wtime()-time;
1132                 }
1133 #else
1134                 wallcycle_stop(times, ewcPME_FFTCOMM);
1135 #endif
1136             } /*master*/
1137         }     /* bPrallelDim */
1138 #pragma omp barrier  /*both needed for parallel and non-parallel dimension (either have to wait on data from AlltoAll or from last FFT*/
1139
1140         /* ---------- END SPLIT + TRANSPOSE------------ */
1141
1142         /* ---------- START JOIN ------------ */
1143 #ifdef NOGMX
1144         if (times != NULL && thread == 0)
1145         {
1146             time = MPI_Wtime();
1147         }
1148 #endif
1149
1150         if (bParallelDim)
1151         {
1152             joinin = lout3;
1153         }
1154         else
1155         {
1156             joinin = fftout;
1157         }
1158         /*bring back in matrix form
1159            thus make  new 1. axes contiguos
1160            also local transpose 1 and 2/3
1161            runs on thread used for following FFT (thus needing a barrier before but not afterwards)
1162          */
1163         if ((s == 0 && !(plan->flags&FFT5D_ORDER_YZ)) || (s == 1 && (plan->flags&FFT5D_ORDER_YZ)))
1164         {
1165             if (pM[s] > 0)
1166             {
1167                 tstart = ( thread   *pM[s]*pN[s]/plan->nthreads);
1168                 tend   = ((thread+1)*pM[s]*pN[s]/plan->nthreads);
1169                 joinAxesTrans13(lin, joinin, N[s], pM[s], K[s], pM[s], P[s], C[s+1], iNin[s+1], oNin[s+1], tstart%pM[s], tstart/pM[s], tend%pM[s], tend/pM[s]);
1170             }
1171         }
1172         else
1173         {
1174             if (pN[s] > 0)
1175             {
1176                 tstart = ( thread   *pK[s]*pN[s]/plan->nthreads);
1177                 tend   = ((thread+1)*pK[s]*pN[s]/plan->nthreads);
1178                 joinAxesTrans12(lin, joinin, N[s], M[s], pK[s], pN[s], P[s], C[s+1], iNin[s+1], oNin[s+1], tstart%pN[s], tstart/pN[s], tend%pN[s], tend/pN[s]);
1179             }
1180         }
1181
1182 #ifdef NOGMX
1183         if (times != NULL && thread == 0)
1184         {
1185             time_local += MPI_Wtime()-time;
1186         }
1187 #endif
1188         if (plan->flags&FFT5D_DEBUG && thread == 0)
1189         {
1190             print_localdata(lin, "%d %d: tranposed %d\n", s+1, plan);
1191         }
1192         /* ---------- END JOIN ------------ */
1193
1194         /*if (debug) print_localdata(lin, "%d %d: transposed x-z\n", N1, M0, K, ZYX, coor);*/
1195     }  /* for(s=0;s<2;s++) */
1196 #ifdef NOGMX
1197     if (times != NULL && thread == 0)
1198     {
1199         time = MPI_Wtime();
1200     }
1201 #endif
1202
1203     if (plan->flags&FFT5D_INPLACE)
1204     {
1205         lout = lin;                          /*in place currently not supported*/
1206
1207     }
1208     /*  ----------- FFT ----------- */
1209     tstart = (thread*pM[s]*pK[s]/plan->nthreads)*C[s];
1210     if ((plan->flags&FFT5D_REALCOMPLEX) && (plan->flags&FFT5D_BACKWARD))
1211     {
1212         gmx_fft_many_1d_real(p1d[s][thread], (plan->flags&FFT5D_BACKWARD) ? GMX_FFT_COMPLEX_TO_REAL : GMX_FFT_REAL_TO_COMPLEX, lin+tstart, lout+tstart);
1213     }
1214     else
1215     {
1216         gmx_fft_many_1d(     p1d[s][thread], (plan->flags&FFT5D_BACKWARD) ? GMX_FFT_BACKWARD : GMX_FFT_FORWARD,               lin+tstart, lout+tstart);
1217     }
1218     /* ------------ END FFT ---------*/
1219
1220 #ifdef NOGMX
1221     if (times != NULL && thread == 0)
1222     {
1223         time_fft += MPI_Wtime()-time;
1224
1225         times->fft   += time_fft;
1226         times->local += time_local;
1227         times->mpi2  += time_mpi[1];
1228         times->mpi1  += time_mpi[0];
1229     }
1230 #endif
1231
1232     if (plan->flags&FFT5D_DEBUG && thread == 0)
1233     {
1234         print_localdata(lout, "%d %d: FFT %d\n", s, plan);
1235     }
1236 }
1237
1238 void fft5d_destroy(fft5d_plan plan)
1239 {
1240     int s, t;
1241
1242     for (s = 0; s < 3; s++)
1243     {
1244         if (plan->p1d[s])
1245         {
1246             for (t = 0; t < plan->nthreads; t++)
1247             {
1248                 gmx_many_fft_destroy(plan->p1d[s][t]);
1249             }
1250             free(plan->p1d[s]);
1251         }
1252         if (plan->iNin[s])
1253         {
1254             free(plan->iNin[s]);
1255             plan->iNin[s] = 0;
1256         }
1257         if (plan->oNin[s])
1258         {
1259             free(plan->oNin[s]);
1260             plan->oNin[s] = 0;
1261         }
1262         if (plan->iNout[s])
1263         {
1264             free(plan->iNout[s]);
1265             plan->iNout[s] = 0;
1266         }
1267         if (plan->oNout[s])
1268         {
1269             free(plan->oNout[s]);
1270             plan->oNout[s] = 0;
1271         }
1272     }
1273 #ifdef GMX_FFT_FFTW3
1274     FFTW_LOCK;
1275 #ifdef FFT5D_MPI_TRANSPOS
1276     for (s = 0; s < 2; s++)
1277     {
1278         FFTW(destroy_plan)(plan->mpip[s]);
1279     }
1280 #endif /* FFT5D_MPI_TRANSPOS */
1281     if (plan->p3d)
1282     {
1283         FFTW(destroy_plan)(plan->p3d);
1284     }
1285     FFTW_UNLOCK;
1286 #endif /* GMX_FFT_FFTW3 */
1287
1288     if (!(plan->flags&FFT5D_NOMALLOC))
1289     {
1290         sfree_aligned(plan->lin);
1291         sfree_aligned(plan->lout);
1292         if (plan->nthreads > 1)
1293         {
1294             sfree_aligned(plan->lout2);
1295             sfree_aligned(plan->lout3);
1296         }
1297     }
1298
1299 #ifdef FFT5D_THREADS
1300 #ifdef FFT5D_FFTW_THREADS
1301     /*FFTW(cleanup_threads)();*/
1302 #endif
1303 #endif
1304
1305     free(plan);
1306 }
1307
1308 /*Is this better than direct access of plan? enough data?
1309    here 0,1 reference divided by which processor grid dimension (not FFT step!)*/
1310 void fft5d_local_size(fft5d_plan plan, int* N1, int* M0, int* K0, int* K1, int** coor)
1311 {
1312     *N1 = plan->N[0];
1313     *M0 = plan->M[0];
1314     *K1 = plan->K[0];
1315     *K0 = plan->N[1];
1316
1317     *coor = plan->coor;
1318 }
1319
1320
1321 /*same as fft5d_plan_3d but with cartesian coordinator and automatic splitting
1322    of processor dimensions*/
1323 fft5d_plan fft5d_plan_3d_cart(int NG, int MG, int KG, MPI_Comm comm, int P0, int flags, t_complex** rlin, t_complex** rlout, t_complex** rlout2, t_complex** rlout3, int nthreads)
1324 {
1325     MPI_Comm cart[2] = {0};
1326 #ifdef GMX_MPI
1327     int      size = 1, prank = 0;
1328     int      P[2];
1329     int      coor[2];
1330     int      wrap[] = {0, 0};
1331     MPI_Comm gcart;
1332     int      rdim1[] = {0, 1}, rdim2[] = {1, 0};
1333
1334     MPI_Comm_size(comm, &size);
1335     MPI_Comm_rank(comm, &prank);
1336
1337     if (P0 == 0)
1338     {
1339         P0 = lfactor(size);
1340     }
1341     if (size%P0 != 0)
1342     {
1343         if (prank == 0)
1344         {
1345             printf("FFT5D: WARNING: Number of ranks %d not evenly divisible by %d\n", size, P0);
1346         }
1347         P0 = lfactor(size);
1348     }
1349
1350     P[0] = P0; P[1] = size/P0; /*number of processors in the two dimensions*/
1351
1352     /*Difference between x-y-z regarding 2d decomposition is whether they are
1353        distributed along axis 1, 2 or both*/
1354
1355     MPI_Cart_create(comm, 2, P, wrap, 1, &gcart); /*parameter 4: value 1: reorder*/
1356     MPI_Cart_get(gcart, 2, P, wrap, coor);
1357     MPI_Cart_sub(gcart, rdim1, &cart[0]);
1358     MPI_Cart_sub(gcart, rdim2, &cart[1]);
1359 #else
1360     (void)P0;
1361     (void)comm;
1362 #endif
1363     return fft5d_plan_3d(NG, MG, KG, cart, flags, rlin, rlout, rlout2, rlout3, nthreads);
1364 }
1365
1366
1367
1368 /*prints in original coordinate system of data (as the input to FFT)*/
1369 void fft5d_compare_data(const t_complex* lin, const t_complex* in, fft5d_plan plan, int bothLocal, int normalize)
1370 {
1371     int  xs[3], xl[3], xc[3], NG[3];
1372     int  x, y, z, l;
1373     int *coor = plan->coor;
1374     int  ll   = 2; /*compare ll values per element (has to be 2 for complex)*/
1375     if ((plan->flags&FFT5D_REALCOMPLEX) && (plan->flags&FFT5D_BACKWARD))
1376     {
1377         ll = 1;
1378     }
1379
1380     compute_offsets(plan, xs, xl, xc, NG, 2);
1381     if (plan->flags&FFT5D_DEBUG)
1382     {
1383         printf("Compare2\n");
1384     }
1385     for (z = 0; z < xl[2]; z++)
1386     {
1387         for (y = 0; y < xl[1]; y++)
1388         {
1389             if (plan->flags&FFT5D_DEBUG)
1390             {
1391                 printf("%d %d: ", coor[0], coor[1]);
1392             }
1393             for (x = 0; x < xl[0]; x++)
1394             {
1395                 for (l = 0; l < ll; l++)   /*loop over real/complex parts*/
1396                 {
1397                     real a, b;
1398                     a = ((real*)lin)[(z*xs[2]+y*xs[1])*2+x*xs[0]*ll+l];
1399                     if (normalize)
1400                     {
1401                         a /= plan->rC[0]*plan->rC[1]*plan->rC[2];
1402                     }
1403                     if (!bothLocal)
1404                     {
1405                         b = ((real*)in)[((z+xc[2])*NG[0]*NG[1]+(y+xc[1])*NG[0])*2+(x+xc[0])*ll+l];
1406                     }
1407                     else
1408                     {
1409                         b = ((real*)in)[(z*xs[2]+y*xs[1])*2+x*xs[0]*ll+l];
1410                     }
1411                     if (plan->flags&FFT5D_DEBUG)
1412                     {
1413                         printf("%f %f, ", a, b);
1414                     }
1415                     else
1416                     {
1417                         if (fabs(a-b) > 2*NG[0]*NG[1]*NG[2]*GMX_REAL_EPS)
1418                         {
1419                             printf("result incorrect on %d,%d at %d,%d,%d: FFT5D:%f reference:%f\n", coor[0], coor[1], x, y, z, a, b);
1420                         }
1421 /*                        assert(fabs(a-b)<2*NG[0]*NG[1]*NG[2]*GMX_REAL_EPS);*/
1422                     }
1423                 }
1424                 if (plan->flags&FFT5D_DEBUG)
1425                 {
1426                     printf(",");
1427                 }
1428             }
1429             if (plan->flags&FFT5D_DEBUG)
1430             {
1431                 printf("\n");
1432             }
1433         }
1434     }
1435
1436 }