src/gromacs/fft/parallel_3dfft.cpp

   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 "parallel_3dfft.h"
  39
  40 #include <errno.h>
  41 #include <stdlib.h>
  42 #include <string.h>
  43
  44 #include "gromacs/fft/fft.h"
  45 #include "gromacs/fft/fft5d.h"
  46 #include "gromacs/math/gmxcomplex.h"
  47 #include "gromacs/utility/fatalerror.h"
  48 #include "gromacs/utility/gmxmpi.h"
  49 #include "gromacs/utility/smalloc.h"
  50
  51 struct gmx_parallel_3dfft  {
  52     fft5d_plan p1, p2;
  53 };
  54
  55 int
  56 gmx_parallel_3dfft_init   (gmx_parallel_3dfft_t     *    pfft_setup,
  57                            ivec                          ndata,
  58                            real     **                   real_data,
  59                            t_complex     **              complex_data,
  60                            MPI_Comm                      comm[2],
  61                            gmx_bool                      bReproducible,
  62                            int                           nthreads)
  63 {
  64     int        rN      = ndata[2], M = ndata[1], K = ndata[0];
  65     int        flags   = FFT5D_REALCOMPLEX | FFT5D_ORDER_YZ; /* FFT5D_DEBUG */
  66     MPI_Comm   rcomm[] = {comm[1], comm[0]};
  67     int        Nb, Mb, Kb;                                   /* dimension for backtransform (in starting order) */
  68     t_complex *buf1, *buf2;                                  /*intermediate buffers - used internally.*/
  69
  70     snew(*pfft_setup, 1);
  71     if (bReproducible)
  72     {
  73         flags |= FFT5D_NOMEASURE;
  74     }
  75
  76     if (!(flags&FFT5D_ORDER_YZ))
  77     {
  78         Nb = M; Mb = K; Kb = rN;
  79     }
  80     else
  81     {
  82         Nb = K; Mb = rN; Kb = M;  /* currently always true because ORDER_YZ always set */
  83     }
  84
  85     (*pfft_setup)->p1 = fft5d_plan_3d(rN, M, K, rcomm, flags, (t_complex**)real_data, complex_data, &buf1, &buf2, nthreads);
  86
  87     (*pfft_setup)->p2 = fft5d_plan_3d(Nb, Mb, Kb, rcomm,
  88                                       (flags|FFT5D_BACKWARD|FFT5D_NOMALLOC)^FFT5D_ORDER_YZ, complex_data, (t_complex**)real_data, &buf1, &buf2, nthreads);
  89
  90     return (*pfft_setup)->p1 != 0 && (*pfft_setup)->p2 != 0;
  91 }
  92
  93
  94 static int
  95 fft5d_limits(fft5d_plan                p,
  96              ivec                      local_ndata,
  97              ivec                      local_offset,
  98              ivec                      local_size)
  99 {
 100     local_offset[2] = 0;
 101     local_offset[1] = p->oM[0];  /*=p->coor[0]*p->MG/p->P[0]; */
 102     local_offset[0] = p->oK[0];  /*=p->coor[1]*p->KG/p->P[1]; */
 103
 104     local_ndata[2] = p->rC[0];
 105     local_ndata[1] = p->pM[0];
 106     local_ndata[0] = p->pK[0];
 107
 108     if ((!(p->flags&FFT5D_BACKWARD)) && (p->flags&FFT5D_REALCOMPLEX))
 109     {
 110         //C is length in multiples of complex local_size in multiples of real
 111         local_size[2] = p->C[0]*2;
 112     }
 113     else
 114     {
 115         local_size[2] = p->C[0];
 116     }
 117     local_size[1] = p->pM[0];
 118     local_size[0] = p->pK[0];
 119     return 0;
 120 }
 121
 122 int
 123 gmx_parallel_3dfft_real_limits(gmx_parallel_3dfft_t      pfft_setup,
 124                                ivec                      local_ndata,
 125                                ivec                      local_offset,
 126                                ivec                      local_size)
 127 {
 128     return fft5d_limits(pfft_setup->p1, local_ndata, local_offset, local_size);
 129 }
 130
 131 static void reorder_ivec_yzx(ivec v)
 132 {
 133     int tmp;
 134
 135     tmp   = v[0];
 136     v[XX] = v[2];
 137     v[ZZ] = v[1];
 138     v[YY] = tmp;
 139 }
 140
 141 int
 142 gmx_parallel_3dfft_complex_limits(gmx_parallel_3dfft_t      pfft_setup,
 143                                   ivec                      complex_order,
 144                                   ivec                      local_ndata,
 145                                   ivec                      local_offset,
 146                                   ivec                      local_size)
 147 {
 148     int ret;
 149
 150     /* For now everything is in-order, but prepare to save communication by avoiding transposes */
 151     complex_order[0] = 0;
 152     complex_order[1] = 1;
 153     complex_order[2] = 2;
 154
 155     ret = fft5d_limits(pfft_setup->p2, local_ndata, local_offset, local_size);
 156
 157     reorder_ivec_yzx(local_ndata);
 158     reorder_ivec_yzx(local_offset);
 159     reorder_ivec_yzx(local_size);
 160
 161     return ret;
 162 }
 163
 164
 165 int
 166 gmx_parallel_3dfft_execute(gmx_parallel_3dfft_t    pfft_setup,
 167                            enum gmx_fft_direction  dir,
 168                            int                     thread,
 169                            gmx_wallcycle_t         wcycle)
 170 {
 171     if ((!(pfft_setup->p1->flags&FFT5D_REALCOMPLEX)) ^ (dir == GMX_FFT_FORWARD || dir == GMX_FFT_BACKWARD))
 172     {
 173         gmx_fatal(FARGS, "Invalid transform. Plan and execution don't match regarding reel/complex");
 174     }
 175     if (dir == GMX_FFT_FORWARD || dir == GMX_FFT_REAL_TO_COMPLEX)
 176     {
 177         fft5d_execute(pfft_setup->p1, thread, wcycle);
 178     }
 179     else
 180     {
 181         fft5d_execute(pfft_setup->p2, thread, wcycle);
 182     }
 183     return 0;
 184 }
 185
 186 int
 187 gmx_parallel_3dfft_destroy(gmx_parallel_3dfft_t    pfft_setup)
 188 {
 189     if (pfft_setup)
 190     {
 191         fft5d_destroy(pfft_setup->p2);
 192         fft5d_destroy(pfft_setup->p1);
 193         sfree(pfft_setup);
 194     }
 195     return 0;
 196 }