src/gromacs/fft/parallel_3dfft.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2005 David van der Spoel, Erik Lindahl, University of Groningen.
   5  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36 #include "gmxpre.h"
  37
  38 #include "config.h"
  39
  40 #include <stdlib.h>
  41 #include <string.h>
  42 #include <errno.h>
  43
  44 #include "gromacs/fft/fft.h"
  45 #include "gromacs/fft/parallel_3dfft.h"
  46 #include "gromacs/math/gmxcomplex.h"
  47 #include "gromacs/utility/gmxmpi.h"
  48
  49 #include "gromacs/utility/smalloc.h"
  50 #include "gromacs/utility/fatalerror.h"
  51
  52 #include "fft5d.h"
  53
  54 struct gmx_parallel_3dfft  {
  55     fft5d_plan p1, p2;
  56 };
  57
  58 int
  59 gmx_parallel_3dfft_init   (gmx_parallel_3dfft_t     *    pfft_setup,
  60                            ivec                          ndata,
  61                            real     **                   real_data,
  62                            t_complex     **              complex_data,
  63                            MPI_Comm                      comm[2],
  64                            gmx_bool                      bReproducible,
  65                            int                           nthreads)
  66 {
  67     int        rN      = ndata[2], M = ndata[1], K = ndata[0];
  68     int        flags   = FFT5D_REALCOMPLEX | FFT5D_ORDER_YZ; /* FFT5D_DEBUG */
  69     MPI_Comm   rcomm[] = {comm[1], comm[0]};
  70     int        Nb, Mb, Kb;                                   /* dimension for backtransform (in starting order) */
  71     t_complex *buf1, *buf2;                                  /*intermediate buffers - used internally.*/
  72
  73     snew(*pfft_setup, 1);
  74     if (bReproducible)
  75     {
  76         flags |= FFT5D_NOMEASURE;
  77     }
  78
  79     if (!(flags&FFT5D_ORDER_YZ))
  80     {
  81         Nb = M; Mb = K; Kb = rN;
  82     }
  83     else
  84     {
  85         Nb = K; Mb = rN; Kb = M;  /* currently always true because ORDER_YZ always set */
  86     }
  87
  88     (*pfft_setup)->p1 = fft5d_plan_3d(rN, M, K, rcomm, flags, (t_complex**)real_data, complex_data, &buf1, &buf2, nthreads);
  89
  90     (*pfft_setup)->p2 = fft5d_plan_3d(Nb, Mb, Kb, rcomm,
  91                                       (flags|FFT5D_BACKWARD|FFT5D_NOMALLOC)^FFT5D_ORDER_YZ, complex_data, (t_complex**)real_data, &buf1, &buf2, nthreads);
  92
  93     return (*pfft_setup)->p1 != 0 && (*pfft_setup)->p2 != 0;
  94 }
  95
  96
  97 static int
  98 fft5d_limits(fft5d_plan                p,
  99              ivec                      local_ndata,
 100              ivec                      local_offset,
 101              ivec                      local_size)
 102 {
 103     local_offset[2] = 0;
 104     local_offset[1] = p->oM[0];  /*=p->coor[0]*p->MG/p->P[0]; */
 105     local_offset[0] = p->oK[0];  /*=p->coor[1]*p->KG/p->P[1]; */
 106
 107     local_ndata[2] = p->rC[0];
 108     local_ndata[1] = p->pM[0];
 109     local_ndata[0] = p->pK[0];
 110
 111     if ((!(p->flags&FFT5D_BACKWARD)) && (p->flags&FFT5D_REALCOMPLEX))
 112     {
 113         //C is length in multiples of complex local_size in multiples of real
 114         local_size[2] = p->C[0]*2;
 115     }
 116     else
 117     {
 118         local_size[2] = p->C[0];
 119     }
 120     local_size[1] = p->pM[0];
 121     local_size[0] = p->pK[0];
 122     return 0;
 123 }
 124
 125 int
 126 gmx_parallel_3dfft_real_limits(gmx_parallel_3dfft_t      pfft_setup,
 127                                ivec                      local_ndata,
 128                                ivec                      local_offset,
 129                                ivec                      local_size)
 130 {
 131     return fft5d_limits(pfft_setup->p1, local_ndata, local_offset, local_size);
 132 }
 133
 134 static void reorder_ivec_yzx(ivec v)
 135 {
 136     real tmp;
 137
 138     tmp   = v[0];
 139     v[XX] = v[2];
 140     v[ZZ] = v[1];
 141     v[YY] = tmp;
 142 }
 143
 144 int
 145 gmx_parallel_3dfft_complex_limits(gmx_parallel_3dfft_t      pfft_setup,
 146                                   ivec                      complex_order,
 147                                   ivec                      local_ndata,
 148                                   ivec                      local_offset,
 149                                   ivec                      local_size)
 150 {
 151     int ret;
 152
 153     /* For now everything is in-order, but prepare to save communication by avoiding transposes */
 154     complex_order[0] = 0;
 155     complex_order[1] = 1;
 156     complex_order[2] = 2;
 157
 158     ret = fft5d_limits(pfft_setup->p2, local_ndata, local_offset, local_size);
 159
 160     reorder_ivec_yzx(local_ndata);
 161     reorder_ivec_yzx(local_offset);
 162     reorder_ivec_yzx(local_size);
 163
 164     return ret;
 165 }
 166
 167
 168 int
 169 gmx_parallel_3dfft_execute(gmx_parallel_3dfft_t    pfft_setup,
 170                            enum gmx_fft_direction  dir,
 171                            int                     thread,
 172                            gmx_wallcycle_t         wcycle)
 173 {
 174     if ((!(pfft_setup->p1->flags&FFT5D_REALCOMPLEX)) ^ (dir == GMX_FFT_FORWARD || dir == GMX_FFT_BACKWARD))
 175     {
 176         gmx_fatal(FARGS, "Invalid transform. Plan and execution don't match regarding reel/complex");
 177     }
 178     if (dir == GMX_FFT_FORWARD || dir == GMX_FFT_REAL_TO_COMPLEX)
 179     {
 180         fft5d_execute(pfft_setup->p1, thread, wcycle);
 181     }
 182     else
 183     {
 184         fft5d_execute(pfft_setup->p2, thread, wcycle);
 185     }
 186     return 0;
 187 }
 188
 189 int
 190 gmx_parallel_3dfft_destroy(gmx_parallel_3dfft_t    pfft_setup)
 191 {
 192     if (pfft_setup)
 193     {
 194         fft5d_destroy(pfft_setup->p2);
 195         fft5d_destroy(pfft_setup->p1);
 196         sfree(pfft_setup);
 197     }
 198     return 0;
 199 }