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,2015, 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 "gmxpre.h"
  36
  37 #include "fft5d.h"
  38
  39 #include "config.h"
  40
  41 #include <assert.h>
  42 #include <float.h>
  43 #include <math.h>
  44 #include <stdio.h>
  45 #include <stdlib.h>
  46 #include <string.h>
  47
  48 #include <algorithm>
  49
  50 #include "gromacs/utility/fatalerror.h"
  51 #include "gromacs/utility/gmxmpi.h"
  52 #include "gromacs/utility/smalloc.h"
  53
  54 #ifdef NOGMX
  55 #define GMX_PARALLEL_ENV_INITIALIZED 1
  56 #else
  57 #ifdef GMX_MPI
  58 #define GMX_PARALLEL_ENV_INITIALIZED 1
  59 #else
  60 #define GMX_PARALLEL_ENV_INITIALIZED 0
  61 #endif
  62 #endif
  63
  64 #ifdef GMX_OPENMP
  65 /* TODO: Do we still need this? Are we still planning ot use fftw + OpenMP? */
  66 #define FFT5D_THREADS
  67 /* requires fftw compiled with openmp */
  68 /* #define FFT5D_FFTW_THREADS (now set by cmake) */
  69 #endif
  70
  71 #ifndef __FLT_EPSILON__
  72 #define __FLT_EPSILON__ FLT_EPSILON
  73 #define __DBL_EPSILON__ DBL_EPSILON
  74 #endif
  75
  76 #ifdef NOGMX
  77 FILE* debug = 0;
  78 #endif
  79
  80 #ifdef GMX_FFT_FFTW3
  81 #include "gromacs/utility/exceptions.h"
  82 #include "gromacs/utility/mutex.h"
  83 /* none of the fftw3 calls, except execute(), are thread-safe, so
  84    we need to serialize them with this mutex. */
  85 static gmx::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 #ifdef GMX_MPI
  91 /* largest factor smaller than sqrt */
  92 static int lfactor(int z)
  93 {
  94     int i = static_cast<int>(sqrt(static_cast<double>(z)));
  95     while (z%i != 0)
  96     {
  97         i--;
  98     }
  99     return i;
 100 }
 101 #endif
 102
 103 /* largest factor */
 104 static int l2factor(int z)
 105 {
 106     int i;
 107     if (z == 1)
 108     {
 109         i = 1;
 110     }
 111     else
 112     {
 113         i = z/2;
 114         while (z%i != 0)
 115         {
 116             i--;
 117         }
 118     }
 119     return i;
 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
 455     /* Don't add more stuff here! We have already had at least one bug because we are reimplementing
 456      * the low-level FFT interface instead of using the Gromacs FFT module. If we need more
 457      * generic functionality it is far better to extend the interface so we can use it for
 458      * all FFT libraries instead of writing FFTW-specific code here.
 459      */
 460
 461     /*if not FFTW - then we don't do a 3d plan but instead use only 1D plans */
 462     /* It is possible to use the 3d plan with OMP threads - but in that case it is not allowed to be called from
 463      * within a parallel region. For now deactivated. If it should be supported it has to made sure that
 464      * that the execute of the 3d plan is in a master/serial block (since it contains it own parallel region)
 465      * and that the 3d plan is faster than the 1d plan.
 466      */
 467     if ((!(flags&FFT5D_INPLACE)) && (!(P[0] > 1 || P[1] > 1)) && nthreads == 1) /*don't do 3d plan in parallel or if in_place requested */
 468     {
 469         int fftwflags = FFTW_DESTROY_INPUT;
 470         FFTW(iodim) dims[3];
 471         int inNG = NG, outMG = MG, outKG = KG;
 472
 473         FFTW_LOCK;
 474
 475         fftwflags |= (flags & FFT5D_NOMEASURE) ? FFTW_ESTIMATE : FFTW_MEASURE;
 476
 477         if (flags&FFT5D_REALCOMPLEX)
 478         {
 479             if (!(flags&FFT5D_BACKWARD))        /*input pointer is not complex*/
 480             {
 481                 inNG *= 2;
 482             }
 483             else                                /*output pointer is not complex*/
 484             {
 485                 if (!(flags&FFT5D_ORDER_YZ))
 486                 {
 487                     outMG *= 2;
 488                 }
 489                 else
 490                 {
 491                     outKG *= 2;
 492                 }
 493             }
 494         }
 495
 496         if (!(flags&FFT5D_BACKWARD))
 497         {
 498             dims[0].n  = KG;
 499             dims[1].n  = MG;
 500             dims[2].n  = rNG;
 501
 502             dims[0].is = inNG*MG;         /*N M K*/
 503             dims[1].is = inNG;
 504             dims[2].is = 1;
 505             if (!(flags&FFT5D_ORDER_YZ))
 506             {
 507                 dims[0].os = MG;           /*M K N*/
 508                 dims[1].os = 1;
 509                 dims[2].os = MG*KG;
 510             }
 511             else
 512             {
 513                 dims[0].os = 1;           /*K N M*/
 514                 dims[1].os = KG*NG;
 515                 dims[2].os = KG;
 516             }
 517         }
 518         else
 519         {
 520             if (!(flags&FFT5D_ORDER_YZ))
 521             {
 522                 dims[0].n  = NG;
 523                 dims[1].n  = KG;
 524                 dims[2].n  = rMG;
 525
 526                 dims[0].is = 1;
 527                 dims[1].is = NG*MG;
 528                 dims[2].is = NG;
 529
 530                 dims[0].os = outMG*KG;
 531                 dims[1].os = outMG;
 532                 dims[2].os = 1;
 533             }
 534             else
 535             {
 536                 dims[0].n  = MG;
 537                 dims[1].n  = NG;
 538                 dims[2].n  = rKG;
 539
 540                 dims[0].is = NG;
 541                 dims[1].is = 1;
 542                 dims[2].is = NG*MG;
 543
 544                 dims[0].os = outKG*NG;
 545                 dims[1].os = outKG;
 546                 dims[2].os = 1;
 547             }
 548         }
 549 #ifdef FFT5D_THREADS
 550 #ifdef FFT5D_FFTW_THREADS
 551         FFTW(plan_with_nthreads)(nthreads);
 552 #endif
 553 #endif
 554         if ((flags&FFT5D_REALCOMPLEX) && !(flags&FFT5D_BACKWARD))
 555         {
 556             plan->p3d = FFTW(plan_guru_dft_r2c)(/*rank*/ 3, dims,
 557                                                          /*howmany*/ 0, /*howmany_dims*/ 0,
 558                                                          (real*)lin, (FFTW(complex) *) lout,
 559                                                          /*flags*/ fftwflags);
 560         }
 561         else if ((flags&FFT5D_REALCOMPLEX) && (flags&FFT5D_BACKWARD))
 562         {
 563             plan->p3d = FFTW(plan_guru_dft_c2r)(/*rank*/ 3, dims,
 564                                                          /*howmany*/ 0, /*howmany_dims*/ 0,
 565                                                          (FFTW(complex) *) lin, (real*)lout,
 566                                                          /*flags*/ fftwflags);
 567         }
 568         else
 569         {
 570             plan->p3d = FFTW(plan_guru_dft)(/*rank*/ 3, dims,
 571                                                      /*howmany*/ 0, /*howmany_dims*/ 0,
 572                                                      (FFTW(complex) *) lin, (FFTW(complex) *) lout,
 573                                                      /*sign*/ (flags&FFT5D_BACKWARD) ? 1 : -1, /*flags*/ fftwflags);
 574         }
 575 #ifdef FFT5D_THREADS
 576 #ifdef FFT5D_FFTW_THREADS
 577         FFTW(plan_with_nthreads)(1);
 578 #endif
 579 #endif
 580         FFTW_UNLOCK;
 581     }
 582     if (!plan->p3d) /* for decomposition and if 3d plan did not work */
 583     {
 584 #endif              /* GMX_FFT_FFTW3 */
 585     for (s = 0; s < 3; s++)
 586     {
 587         if (debug)
 588         {
 589             fprintf(debug, "FFT5D: Plan s %d rC %d M %d pK %d C %d lsize %d\n",
 590                     s, rC[s], M[s], pK[s], C[s], lsize);
 591         }
 592         plan->p1d[s] = (gmx_fft_t*)malloc(sizeof(gmx_fft_t)*nthreads);
 593
 594         /* Make sure that the init routines are only called by one thread at a time and in order
 595            (later is only important to not confuse valgrind)
 596          */
 597 #pragma omp parallel for num_threads(nthreads) schedule(static) ordered
 598         for (t = 0; t < nthreads; t++)
 599         {
 600 #pragma omp ordered
 601             {
 602                 int tsize = ((t+1)*pM[s]*pK[s]/nthreads)-(t*pM[s]*pK[s]/nthreads);
 603
 604                 if ((flags&FFT5D_REALCOMPLEX) && ((!(flags&FFT5D_BACKWARD) && s == 0) || ((flags&FFT5D_BACKWARD) && s == 2)))
 605                 {
 606                     gmx_fft_init_many_1d_real( &plan->p1d[s][t], rC[s], tsize, (flags&FFT5D_NOMEASURE) ? GMX_FFT_FLAG_CONSERVATIVE : 0 );
 607                 }
 608                 else
 609                 {
 610                     gmx_fft_init_many_1d     ( &plan->p1d[s][t],  C[s], tsize, (flags&FFT5D_NOMEASURE) ? GMX_FFT_FLAG_CONSERVATIVE : 0 );
 611                 }
 612             }
 613         }
 614     }
 615
 616 #ifdef GMX_FFT_FFTW3
 617 }
 618 #endif
 619     if ((flags&FFT5D_ORDER_YZ))   /*plan->cart is in the order of transposes */
 620     {
 621         plan->cart[0] = comm[0]; plan->cart[1] = comm[1];
 622     }
 623     else
 624     {
 625         plan->cart[1] = comm[0]; plan->cart[0] = comm[1];
 626     }
 627 #ifdef FFT5D_MPI_TRANSPOSE
 628     FFTW_LOCK;
 629     for (s = 0; s < 2; s++)
 630     {
 631         if ((s == 0 && !(flags&FFT5D_ORDER_YZ)) || (s == 1 && (flags&FFT5D_ORDER_YZ)))
 632         {
 633             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);
 634         }
 635         else
 636         {
 637             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);
 638         }
 639     }
 640     FFTW_UNLOCK;
 641 #endif
 642
 643
 644     plan->lin   = lin;
 645     plan->lout  = lout;
 646     plan->lout2 = lout2;
 647     plan->lout3 = lout3;
 648
 649     plan->NG = NG; plan->MG = MG; plan->KG = KG;
 650
 651     for (s = 0; s < 3; s++)
 652     {
 653         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];
 654         plan->oM[s]   = oM[s]; plan->oK[s] = oK[s];
 655         plan->C[s]    = C[s]; plan->rC[s] = rC[s];
 656         plan->iNin[s] = iNin[s]; plan->oNin[s] = oNin[s]; plan->iNout[s] = iNout[s]; plan->oNout[s] = oNout[s];
 657     }
 658     for (s = 0; s < 2; s++)
 659     {
 660         plan->P[s] = nP[s]; plan->coor[s] = prank[s];
 661     }
 662
 663 /*    plan->fftorder=fftorder;
 664     plan->direction=direction;
 665     plan->realcomplex=realcomplex;
 666  */
 667     plan->flags    = flags;
 668     plan->nthreads = nthreads;
 669     *rlin          = lin;
 670     *rlout         = lout;
 671     *rlout2        = lout2;
 672     *rlout3        = lout3;
 673     return plan;
 674 }
 675
 676
 677 enum order {
 678     XYZ,
 679     XZY,
 680     YXZ,
 681     YZX,
 682     ZXY,
 683     ZYX
 684 };
 685
 686
 687
 688 /*here x,y,z and N,M,K is in rotated coordinate system!!
 689    x (and N) is mayor (consecutive) dimension, y (M) middle and z (K) major
 690    maxN,maxM,maxK is max size of local data
 691    pN, pM, pK is local size specific to current processor (only different to max if not divisible)
 692    NG, MG, KG is size of global data*/
 693 static void splitaxes(t_complex* lout, const t_complex* lin,
 694                       int maxN, int maxM, int maxK, int pM,
 695                       int P, int NG, int *N, int* oN, int starty, int startz, int endy, int endz)
 696 {
 697     int x, y, z, i;
 698     int in_i, out_i, in_z, out_z, in_y, out_y;
 699     int s_y, e_y;
 700
 701     for (z = startz; z < endz+1; z++) /*3. z l*/
 702     {
 703         if (z == startz)
 704         {
 705             s_y = starty;
 706         }
 707         else
 708         {
 709             s_y = 0;
 710         }
 711         if (z == endz)
 712         {
 713             e_y = endy;
 714         }
 715         else
 716         {
 717             e_y = pM;
 718         }
 719         out_z  = z*maxN*maxM;
 720         in_z   = z*NG*pM;
 721
 722         for (i = 0; i < P; i++) /*index cube along long axis*/
 723         {
 724             out_i  = out_z  + i*maxN*maxM*maxK;
 725             in_i   = in_z + oN[i];
 726             for (y = s_y; y < e_y; y++)   /*2. y k*/
 727             {
 728                 out_y  = out_i  + y*maxN;
 729                 in_y   = in_i + y*NG;
 730                 for (x = 0; x < N[i]; x++)       /*1. x j*/
 731                 {
 732                     lout[out_y+x] = lin[in_y+x]; /*in=z*NG*pM+oN[i]+y*NG+x*/
 733                     /*after split important that each processor chunk i has size maxN*maxM*maxK and thus being the same size*/
 734                     /*before split data contiguos - thus if different processor get different amount oN is different*/
 735                 }
 736             }
 737         }
 738     }
 739 }
 740
 741 /*make axis contiguous again (after AllToAll) and also do local transpose*/
 742 /*transpose mayor and major dimension
 743    variables see above
 744    the major, middle, minor order is only correct for x,y,z (N,M,K) for the input
 745    N,M,K local dimensions
 746    KG global size*/
 747 static void joinAxesTrans13(t_complex* lout, const t_complex* lin,
 748                             int maxN, int maxM, int maxK, int pM,
 749                             int P, int KG, int* K, int* oK, int starty, int startx, int endy, int endx)
 750 {
 751     int i, x, y, z;
 752     int out_i, in_i, out_x, in_x, out_z, in_z;
 753     int s_y, e_y;
 754
 755     for (x = startx; x < endx+1; x++) /*1.j*/
 756     {
 757         if (x == startx)
 758         {
 759             s_y = starty;
 760         }
 761         else
 762         {
 763             s_y = 0;
 764         }
 765         if (x == endx)
 766         {
 767             e_y = endy;
 768         }
 769         else
 770         {
 771             e_y = pM;
 772         }
 773
 774         out_x  = x*KG*pM;
 775         in_x   = x;
 776
 777         for (i = 0; i < P; i++) /*index cube along long axis*/
 778         {
 779             out_i  = out_x  + oK[i];
 780             in_i   = in_x + i*maxM*maxN*maxK;
 781             for (z = 0; z < K[i]; z++) /*3.l*/
 782             {
 783                 out_z  = out_i  + z;
 784                 in_z   = in_i + z*maxM*maxN;
 785                 for (y = s_y; y < e_y; y++)              /*2.k*/
 786                 {
 787                     lout[out_z+y*KG] = lin[in_z+y*maxN]; /*out=x*KG*pM+oK[i]+z+y*KG*/
 788                 }
 789             }
 790         }
 791     }
 792 }
 793
 794 /*make axis contiguous again (after AllToAll) and also do local transpose
 795    tranpose mayor and middle dimension
 796    variables see above
 797    the minor, middle, major order is only correct for x,y,z (N,M,K) for the input
 798    N,M,K local size
 799    MG, global size*/
 800 static void joinAxesTrans12(t_complex* lout, const t_complex* lin, int maxN, int maxM, int maxK, int pN,
 801                             int P, int MG, int* M, int* oM, int startx, int startz, int endx, int endz)
 802 {
 803     int i, z, y, x;
 804     int out_i, in_i, out_z, in_z, out_x, in_x;
 805     int s_x, e_x;
 806
 807     for (z = startz; z < endz+1; z++)
 808     {
 809         if (z == startz)
 810         {
 811             s_x = startx;
 812         }
 813         else
 814         {
 815             s_x = 0;
 816         }
 817         if (z == endz)
 818         {
 819             e_x = endx;
 820         }
 821         else
 822         {
 823             e_x = pN;
 824         }
 825         out_z  = z*MG*pN;
 826         in_z   = z*maxM*maxN;
 827
 828         for (i = 0; i < P; i++) /*index cube along long axis*/
 829         {
 830             out_i  = out_z  + oM[i];
 831             in_i   = in_z + i*maxM*maxN*maxK;
 832             for (x = s_x; x < e_x; x++)
 833             {
 834                 out_x  = out_i  + x*MG;
 835                 in_x   = in_i + x;
 836                 for (y = 0; y < M[i]; y++)
 837                 {
 838                     lout[out_x+y] = lin[in_x+y*maxN]; /*out=z*MG*pN+oM[i]+x*MG+y*/
 839                 }
 840             }
 841         }
 842     }
 843 }
 844
 845
 846 static void rotate_offsets(int x[])
 847 {
 848     int t = x[0];
 849 /*    x[0]=x[2];
 850     x[2]=x[1];
 851     x[1]=t;*/
 852     x[0] = x[1];
 853     x[1] = x[2];
 854     x[2] = t;
 855 }
 856
 857 /*compute the offset to compare or print transposed local data in original input coordinates
 858    xs matrix dimension size, xl dimension length, xc decomposition offset
 859    s: step in computation = number of transposes*/
 860 static void compute_offsets(fft5d_plan plan, int xs[], int xl[], int xc[], int NG[], int s)
 861 {
 862 /*    int direction = plan->direction;
 863     int fftorder = plan->fftorder;*/
 864
 865     int  o = 0;
 866     int  pos[3], i;
 867     int *pM = plan->pM, *pK = plan->pK, *oM = plan->oM, *oK = plan->oK,
 868     *C      = plan->C, *rC = plan->rC;
 869
 870     NG[0] = plan->NG; NG[1] = plan->MG; NG[2] = plan->KG;
 871
 872     if (!(plan->flags&FFT5D_ORDER_YZ))
 873     {
 874         switch (s)
 875         {
 876             case 0: o = XYZ; break;
 877             case 1: o = ZYX; break;
 878             case 2: o = YZX; break;
 879             default: assert(0);
 880         }
 881     }
 882     else
 883     {
 884         switch (s)
 885         {
 886             case 0: o = XYZ; break;
 887             case 1: o = YXZ; break;
 888             case 2: o = ZXY; break;
 889             default: assert(0);
 890         }
 891     }
 892
 893     switch (o)
 894     {
 895         case XYZ: pos[0] = 1; pos[1] = 2; pos[2] = 3; break;
 896         case XZY: pos[0] = 1; pos[1] = 3; pos[2] = 2; break;
 897         case YXZ: pos[0] = 2; pos[1] = 1; pos[2] = 3; break;
 898         case YZX: pos[0] = 3; pos[1] = 1; pos[2] = 2; break;
 899         case ZXY: pos[0] = 2; pos[1] = 3; pos[2] = 1; break;
 900         case ZYX: pos[0] = 3; pos[1] = 2; pos[2] = 1; break;
 901     }
 902     /*if (debug) printf("pos: %d %d %d\n",pos[0],pos[1],pos[2]);*/
 903
 904     /*xs, xl give dimension size and data length in local transposed coordinate system
 905        for 0(/1/2): x(/y/z) in original coordinate system*/
 906     for (i = 0; i < 3; i++)
 907     {
 908         switch (pos[i])
 909         {
 910             case 1: xs[i] = 1;         xc[i] = 0;     xl[i] = C[s]; break;
 911             case 2: xs[i] = C[s];      xc[i] = oM[s]; xl[i] = pM[s]; break;
 912             case 3: xs[i] = C[s]*pM[s]; xc[i] = oK[s]; xl[i] = pK[s]; break;
 913         }
 914     }
 915     /*input order is different for test program to match FFTW order
 916        (important for complex to real)*/
 917     if (plan->flags&FFT5D_BACKWARD)
 918     {
 919         rotate_offsets(xs);
 920         rotate_offsets(xl);
 921         rotate_offsets(xc);
 922         rotate_offsets(NG);
 923         if (plan->flags&FFT5D_ORDER_YZ)
 924         {
 925             rotate_offsets(xs);
 926             rotate_offsets(xl);
 927             rotate_offsets(xc);
 928             rotate_offsets(NG);
 929         }
 930     }
 931     if ((plan->flags&FFT5D_REALCOMPLEX) && ((!(plan->flags&FFT5D_BACKWARD) && s == 0) || ((plan->flags&FFT5D_BACKWARD) && s == 2)))
 932     {
 933         xl[0] = rC[s];
 934     }
 935 }
 936
 937 static void print_localdata(const t_complex* lin, const char* txt, int s, fft5d_plan plan)
 938 {
 939     int  x, y, z, l;
 940     int *coor = plan->coor;
 941     int  xs[3], xl[3], xc[3], NG[3];
 942     int  ll = (plan->flags&FFT5D_REALCOMPLEX) ? 1 : 2;
 943     compute_offsets(plan, xs, xl, xc, NG, s);
 944     fprintf(debug, txt, coor[0], coor[1], s);
 945     /*printf("xs: %d %d %d, xl: %d %d %d\n",xs[0],xs[1],xs[2],xl[0],xl[1],xl[2]);*/
 946     for (z = 0; z < xl[2]; z++)
 947     {
 948         for (y = 0; y < xl[1]; y++)
 949         {
 950             fprintf(debug, "%d %d: ", coor[0], coor[1]);
 951             for (x = 0; x < xl[0]; x++)
 952             {
 953                 for (l = 0; l < ll; l++)
 954                 {
 955                     fprintf(debug, "%f ", ((real*)lin)[(z*xs[2]+y*xs[1])*2+(x*xs[0])*ll+l]);
 956                 }
 957                 fprintf(debug, ",");
 958             }
 959             fprintf(debug, "\n");
 960         }
 961     }
 962 }
 963
 964 void fft5d_execute(fft5d_plan plan, int thread, fft5d_time times)
 965 {
 966     t_complex  *lin   = plan->lin;
 967     t_complex  *lout  = plan->lout;
 968     t_complex  *lout2 = plan->lout2;
 969     t_complex  *lout3 = plan->lout3;
 970     t_complex  *fftout, *joinin;
 971
 972     gmx_fft_t **p1d = plan->p1d;
 973 #ifdef FFT5D_MPI_TRANSPOSE
 974     FFTW(plan) *mpip = plan->mpip;
 975 #endif
 976 #ifdef GMX_MPI
 977     MPI_Comm *cart = plan->cart;
 978 #endif
 979 #ifdef NOGMX
 980     double time_fft = 0, time_local = 0, time_mpi[2] = {0}, time = 0;
 981 #endif
 982     int   *N = plan->N, *M = plan->M, *K = plan->K, *pN = plan->pN, *pM = plan->pM, *pK = plan->pK,
 983     *C       = plan->C, *P = plan->P, **iNin = plan->iNin, **oNin = plan->oNin, **iNout = plan->iNout, **oNout = plan->oNout;
 984     int    s = 0, tstart, tend, bParallelDim;
 985
 986
 987 #ifdef GMX_FFT_FFTW3
 988     if (plan->p3d)
 989     {
 990         if (thread == 0)
 991         {
 992 #ifdef NOGMX
 993             if (times != 0)
 994             {
 995                 time = MPI_Wtime();
 996             }
 997 #endif
 998             FFTW(execute)(plan->p3d);
 999 #ifdef NOGMX
1000             if (times != 0)
1001             {
1002                 times->fft += MPI_Wtime()-time;
1003             }
1004 #endif
1005         }
1006         return;
1007     }
1008 #endif
1009
1010     s = 0;
1011
1012     /*lin: x,y,z*/
1013     if ((plan->flags&FFT5D_DEBUG) && thread == 0)
1014     {
1015         print_localdata(lin, "%d %d: copy in lin\n", s, plan);
1016     }
1017
1018     for (s = 0; s < 2; s++)  /*loop over first two FFT steps (corner rotations)*/
1019
1020     {
1021 #ifdef GMX_MPI
1022         if (GMX_PARALLEL_ENV_INITIALIZED && cart[s] != MPI_COMM_NULL && P[s] > 1)
1023         {
1024             bParallelDim = 1;
1025         }
1026         else
1027 #endif
1028         {
1029             bParallelDim = 0;
1030         }
1031
1032         /* ---------- START FFT ------------ */
1033 #ifdef NOGMX
1034         if (times != 0 && thread == 0)
1035         {
1036             time = MPI_Wtime();
1037         }
1038 #endif
1039
1040         if (bParallelDim || plan->nthreads == 1)
1041         {
1042             fftout = lout;
1043         }
1044         else
1045         {
1046             if (s == 0)
1047             {
1048                 fftout = lout3;
1049             }
1050             else
1051             {
1052                 fftout = lout2;
1053             }
1054         }
1055
1056         tstart = (thread*pM[s]*pK[s]/plan->nthreads)*C[s];
1057         if ((plan->flags&FFT5D_REALCOMPLEX) && !(plan->flags&FFT5D_BACKWARD) && s == 0)
1058         {
1059             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);
1060         }
1061         else
1062         {
1063             gmx_fft_many_1d(     p1d[s][thread], (plan->flags&FFT5D_BACKWARD) ? GMX_FFT_BACKWARD : GMX_FFT_FORWARD,               lin+tstart, fftout+tstart);
1064
1065         }
1066
1067 #ifdef NOGMX
1068         if (times != NULL && thread == 0)
1069         {
1070             time_fft += MPI_Wtime()-time;
1071         }
1072 #endif
1073         if ((plan->flags&FFT5D_DEBUG) && thread == 0)
1074         {
1075             print_localdata(lout, "%d %d: FFT %d\n", s, plan);
1076         }
1077         /* ---------- END FFT ------------ */
1078
1079         /* ---------- START SPLIT + TRANSPOSE------------ (if parallel in in this dimension)*/
1080         if (bParallelDim)
1081         {
1082 #ifdef NOGMX
1083             if (times != NULL && thread == 0)
1084             {
1085                 time = MPI_Wtime();
1086             }
1087 #endif
1088             /*prepare for A
1089                llToAll
1090                1. (most outer) axes (x) is split into P[s] parts of size N[s]
1091                for sending*/
1092             if (pM[s] > 0)
1093             {
1094                 tend    = ((thread+1)*pM[s]*pK[s]/plan->nthreads);
1095                 tstart /= C[s];
1096                 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]);
1097             }
1098 #pragma omp barrier /*barrier required before AllToAll (all input has to be their) - before timing to make timing more acurate*/
1099 #ifdef NOGMX
1100             if (times != NULL && thread == 0)
1101             {
1102                 time_local += MPI_Wtime()-time;
1103             }
1104 #endif
1105
1106             /* ---------- END SPLIT , START TRANSPOSE------------ */
1107
1108             if (thread == 0)
1109             {
1110 #ifdef NOGMX
1111                 if (times != 0)
1112                 {
1113                     time = MPI_Wtime();
1114                 }
1115 #else
1116                 wallcycle_start(times, ewcPME_FFTCOMM);
1117 #endif
1118 #ifdef FFT5D_MPI_TRANSPOSE
1119                 FFTW(execute)(mpip[s]);
1120 #else
1121 #ifdef GMX_MPI
1122                 if ((s == 0 && !(plan->flags&FFT5D_ORDER_YZ)) || (s == 1 && (plan->flags&FFT5D_ORDER_YZ)))
1123                 {
1124                     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]);
1125                 }
1126                 else
1127                 {
1128                     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]);
1129                 }
1130 #else
1131                 gmx_incons("fft5d MPI call without MPI configuration");
1132 #endif /*GMX_MPI*/
1133 #endif /*FFT5D_MPI_TRANSPOSE*/
1134 #ifdef NOGMX
1135                 if (times != 0)
1136                 {
1137                     time_mpi[s] = MPI_Wtime()-time;
1138                 }
1139 #else
1140                 wallcycle_stop(times, ewcPME_FFTCOMM);
1141 #endif
1142             } /*master*/
1143         }     /* bPrallelDim */
1144 #pragma omp barrier  /*both needed for parallel and non-parallel dimension (either have to wait on data from AlltoAll or from last FFT*/
1145
1146         /* ---------- END SPLIT + TRANSPOSE------------ */
1147
1148         /* ---------- START JOIN ------------ */
1149 #ifdef NOGMX
1150         if (times != NULL && thread == 0)
1151         {
1152             time = MPI_Wtime();
1153         }
1154 #endif
1155
1156         if (bParallelDim)
1157         {
1158             joinin = lout3;
1159         }
1160         else
1161         {
1162             joinin = fftout;
1163         }
1164         /*bring back in matrix form
1165            thus make  new 1. axes contiguos
1166            also local transpose 1 and 2/3
1167            runs on thread used for following FFT (thus needing a barrier before but not afterwards)
1168          */
1169         if ((s == 0 && !(plan->flags&FFT5D_ORDER_YZ)) || (s == 1 && (plan->flags&FFT5D_ORDER_YZ)))
1170         {
1171             if (pM[s] > 0)
1172             {
1173                 tstart = ( thread   *pM[s]*pN[s]/plan->nthreads);
1174                 tend   = ((thread+1)*pM[s]*pN[s]/plan->nthreads);
1175                 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]);
1176             }
1177         }
1178         else
1179         {
1180             if (pN[s] > 0)
1181             {
1182                 tstart = ( thread   *pK[s]*pN[s]/plan->nthreads);
1183                 tend   = ((thread+1)*pK[s]*pN[s]/plan->nthreads);
1184                 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]);
1185             }
1186         }
1187
1188 #ifdef NOGMX
1189         if (times != NULL && thread == 0)
1190         {
1191             time_local += MPI_Wtime()-time;
1192         }
1193 #endif
1194         if ((plan->flags&FFT5D_DEBUG) && thread == 0)
1195         {
1196             print_localdata(lin, "%d %d: tranposed %d\n", s+1, plan);
1197         }
1198         /* ---------- END JOIN ------------ */
1199
1200         /*if (debug) print_localdata(lin, "%d %d: transposed x-z\n", N1, M0, K, ZYX, coor);*/
1201     }  /* for(s=0;s<2;s++) */
1202 #ifdef NOGMX
1203     if (times != NULL && thread == 0)
1204     {
1205         time = MPI_Wtime();
1206     }
1207 #endif
1208
1209     if (plan->flags&FFT5D_INPLACE)
1210     {
1211         lout = lin;                          /*in place currently not supported*/
1212
1213     }
1214     /*  ----------- FFT ----------- */
1215     tstart = (thread*pM[s]*pK[s]/plan->nthreads)*C[s];
1216     if ((plan->flags&FFT5D_REALCOMPLEX) && (plan->flags&FFT5D_BACKWARD))
1217     {
1218         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);
1219     }
1220     else
1221     {
1222         gmx_fft_many_1d(     p1d[s][thread], (plan->flags&FFT5D_BACKWARD) ? GMX_FFT_BACKWARD : GMX_FFT_FORWARD,               lin+tstart, lout+tstart);
1223     }
1224     /* ------------ END FFT ---------*/
1225
1226 #ifdef NOGMX
1227     if (times != NULL && thread == 0)
1228     {
1229         time_fft += MPI_Wtime()-time;
1230
1231         times->fft   += time_fft;
1232         times->local += time_local;
1233         times->mpi2  += time_mpi[1];
1234         times->mpi1  += time_mpi[0];
1235     }
1236 #endif
1237
1238     if ((plan->flags&FFT5D_DEBUG) && thread == 0)
1239     {
1240         print_localdata(lout, "%d %d: FFT %d\n", s, plan);
1241     }
1242 }
1243
1244 void fft5d_destroy(fft5d_plan plan)
1245 {
1246     int s, t;
1247
1248     for (s = 0; s < 3; s++)
1249     {
1250         if (plan->p1d[s])
1251         {
1252             for (t = 0; t < plan->nthreads; t++)
1253             {
1254                 gmx_many_fft_destroy(plan->p1d[s][t]);
1255             }
1256             free(plan->p1d[s]);
1257         }
1258         if (plan->iNin[s])
1259         {
1260             free(plan->iNin[s]);
1261             plan->iNin[s] = 0;
1262         }
1263         if (plan->oNin[s])
1264         {
1265             free(plan->oNin[s]);
1266             plan->oNin[s] = 0;
1267         }
1268         if (plan->iNout[s])
1269         {
1270             free(plan->iNout[s]);
1271             plan->iNout[s] = 0;
1272         }
1273         if (plan->oNout[s])
1274         {
1275             free(plan->oNout[s]);
1276             plan->oNout[s] = 0;
1277         }
1278     }
1279 #ifdef GMX_FFT_FFTW3
1280     FFTW_LOCK;
1281 #ifdef FFT5D_MPI_TRANSPOS
1282     for (s = 0; s < 2; s++)
1283     {
1284         FFTW(destroy_plan)(plan->mpip[s]);
1285     }
1286 #endif /* FFT5D_MPI_TRANSPOS */
1287     if (plan->p3d)
1288     {
1289         FFTW(destroy_plan)(plan->p3d);
1290     }
1291     FFTW_UNLOCK;
1292 #endif /* GMX_FFT_FFTW3 */
1293
1294     if (!(plan->flags&FFT5D_NOMALLOC))
1295     {
1296         sfree_aligned(plan->lin);
1297         sfree_aligned(plan->lout);
1298         if (plan->nthreads > 1)
1299         {
1300             sfree_aligned(plan->lout2);
1301             sfree_aligned(plan->lout3);
1302         }
1303     }
1304
1305 #ifdef FFT5D_THREADS
1306 #ifdef FFT5D_FFTW_THREADS
1307     /*FFTW(cleanup_threads)();*/
1308 #endif
1309 #endif
1310
1311     free(plan);
1312 }
1313
1314 /*Is this better than direct access of plan? enough data?
1315    here 0,1 reference divided by which processor grid dimension (not FFT step!)*/
1316 void fft5d_local_size(fft5d_plan plan, int* N1, int* M0, int* K0, int* K1, int** coor)
1317 {
1318     *N1 = plan->N[0];
1319     *M0 = plan->M[0];
1320     *K1 = plan->K[0];
1321     *K0 = plan->N[1];
1322
1323     *coor = plan->coor;
1324 }
1325
1326
1327 /*same as fft5d_plan_3d but with cartesian coordinator and automatic splitting
1328    of processor dimensions*/
1329 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)
1330 {
1331     MPI_Comm cart[2] = {0};
1332 #ifdef GMX_MPI
1333     int      size = 1, prank = 0;
1334     int      P[2];
1335     int      coor[2];
1336     int      wrap[] = {0, 0};
1337     MPI_Comm gcart;
1338     int      rdim1[] = {0, 1}, rdim2[] = {1, 0};
1339
1340     MPI_Comm_size(comm, &size);
1341     MPI_Comm_rank(comm, &prank);
1342
1343     if (P0 == 0)
1344     {
1345         P0 = lfactor(size);
1346     }
1347     if (size%P0 != 0)
1348     {
1349         if (prank == 0)
1350         {
1351             printf("FFT5D: WARNING: Number of ranks %d not evenly divisible by %d\n", size, P0);
1352         }
1353         P0 = lfactor(size);
1354     }
1355
1356     P[0] = P0; P[1] = size/P0; /*number of processors in the two dimensions*/
1357
1358     /*Difference between x-y-z regarding 2d decomposition is whether they are
1359        distributed along axis 1, 2 or both*/
1360
1361     MPI_Cart_create(comm, 2, P, wrap, 1, &gcart); /*parameter 4: value 1: reorder*/
1362     MPI_Cart_get(gcart, 2, P, wrap, coor);
1363     MPI_Cart_sub(gcart, rdim1, &cart[0]);
1364     MPI_Cart_sub(gcart, rdim2, &cart[1]);
1365 #else
1366     (void)P0;
1367     (void)comm;
1368 #endif
1369     return fft5d_plan_3d(NG, MG, KG, cart, flags, rlin, rlout, rlout2, rlout3, nthreads);
1370 }
1371
1372
1373
1374 /*prints in original coordinate system of data (as the input to FFT)*/
1375 void fft5d_compare_data(const t_complex* lin, const t_complex* in, fft5d_plan plan, int bothLocal, int normalize)
1376 {
1377     int  xs[3], xl[3], xc[3], NG[3];
1378     int  x, y, z, l;
1379     int *coor = plan->coor;
1380     int  ll   = 2; /*compare ll values per element (has to be 2 for complex)*/
1381     if ((plan->flags&FFT5D_REALCOMPLEX) && (plan->flags&FFT5D_BACKWARD))
1382     {
1383         ll = 1;
1384     }
1385
1386     compute_offsets(plan, xs, xl, xc, NG, 2);
1387     if (plan->flags&FFT5D_DEBUG)
1388     {
1389         printf("Compare2\n");
1390     }
1391     for (z = 0; z < xl[2]; z++)
1392     {
1393         for (y = 0; y < xl[1]; y++)
1394         {
1395             if (plan->flags&FFT5D_DEBUG)
1396             {
1397                 printf("%d %d: ", coor[0], coor[1]);
1398             }
1399             for (x = 0; x < xl[0]; x++)
1400             {
1401                 for (l = 0; l < ll; l++)   /*loop over real/complex parts*/
1402                 {
1403                     real a, b;
1404                     a = ((real*)lin)[(z*xs[2]+y*xs[1])*2+x*xs[0]*ll+l];
1405                     if (normalize)
1406                     {
1407                         a /= plan->rC[0]*plan->rC[1]*plan->rC[2];
1408                     }
1409                     if (!bothLocal)
1410                     {
1411                         b = ((real*)in)[((z+xc[2])*NG[0]*NG[1]+(y+xc[1])*NG[0])*2+(x+xc[0])*ll+l];
1412                     }
1413                     else
1414                     {
1415                         b = ((real*)in)[(z*xs[2]+y*xs[1])*2+x*xs[0]*ll+l];
1416                     }
1417                     if (plan->flags&FFT5D_DEBUG)
1418                     {
1419                         printf("%f %f, ", a, b);
1420                     }
1421                     else
1422                     {
1423                         if (fabs(a-b) > 2*NG[0]*NG[1]*NG[2]*GMX_REAL_EPS)
1424                         {
1425                             printf("result incorrect on %d,%d at %d,%d,%d: FFT5D:%f reference:%f\n", coor[0], coor[1], x, y, z, a, b);
1426                         }
1427 /*                        assert(fabs(a-b)<2*NG[0]*NG[1]*NG[2]*GMX_REAL_EPS);*/
1428                     }
1429                 }
1430                 if (plan->flags&FFT5D_DEBUG)
1431                 {
1432                     printf(",");
1433                 }
1434             }
1435             if (plan->flags&FFT5D_DEBUG)
1436             {
1437                 printf("\n");
1438             }
1439         }
1440     }
1441
1442 }