src/gromacs/fft/fft.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2003 Erik Lindahl, David van der Spoel, University of Groningen.
   5  * Copyright (c) 2013,2014,2015,2017,2018, 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 "fft.h"
  39
  40 #include "config.h"
  41
  42 #include <cerrno>
  43 #include <cstdio>
  44 #include <cstdlib>
  45 #include <cstring>
  46
  47 #include "gromacs/math/gmxcomplex.h"
  48 #include "gromacs/utility/fatalerror.h"
  49 #include "gromacs/utility/real.h"
  50
  51 /* This file contains common fft utility functions, but not
  52  * the actual transform implementations. Check the
  53  * files like fft_fftw3.c or fft_mkl.c for that.
  54  */
  55
  56 #if !GMX_FFT_FFTW3
  57
  58 struct gmx_many_fft {
  59     int       howmany;
  60     int       dist;
  61     gmx_fft_t fft;
  62 };
  63
  64 typedef struct gmx_many_fft* gmx_many_fft_t;
  65
  66 int
  67 gmx_fft_init_many_1d(gmx_fft_t *        pfft,
  68                      int                nx,
  69                      int                howmany,
  70                      gmx_fft_flag       flags)
  71 {
  72     gmx_many_fft_t fft;
  73     if (pfft == nullptr)
  74     {
  75         gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
  76         return EINVAL;
  77     }
  78     *pfft = nullptr;
  79
  80     if ( (fft = (gmx_many_fft_t)malloc(sizeof(struct gmx_many_fft))) == nullptr)
  81     {
  82         return ENOMEM;
  83     }
  84
  85     gmx_fft_init_1d(&fft->fft, nx, flags);
  86     fft->howmany = howmany;
  87     fft->dist    = 2*nx;
  88
  89     *pfft = (gmx_fft_t)fft;
  90     return 0;
  91 }
  92
  93 int
  94 gmx_fft_init_many_1d_real(gmx_fft_t *        pfft,
  95                           int                nx,
  96                           int                howmany,
  97                           gmx_fft_flag       flags)
  98 {
  99     gmx_many_fft_t fft;
 100     if (pfft == nullptr)
 101     {
 102         gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
 103         return EINVAL;
 104     }
 105     *pfft = nullptr;
 106
 107     if ( (fft = (gmx_many_fft_t)malloc(sizeof(struct gmx_many_fft))) == nullptr)
 108     {
 109         return ENOMEM;
 110     }
 111
 112     gmx_fft_init_1d_real(&fft->fft, nx, flags);
 113     fft->howmany = howmany;
 114     fft->dist    = 2*(nx/2+1);
 115
 116     *pfft = (gmx_fft_t)fft;
 117     return 0;
 118 }
 119
 120 int
 121 gmx_fft_many_1d     (gmx_fft_t                  fft,
 122                      enum gmx_fft_direction     dir,
 123                      void *                     in_data,
 124                      void *                     out_data)
 125 {
 126     gmx_many_fft_t mfft = (gmx_many_fft_t)fft;
 127     int            i, ret;
 128     for (i = 0; i < mfft->howmany; i++)
 129     {
 130         ret = gmx_fft_1d(mfft->fft, dir, in_data, out_data);
 131         if (ret != 0)
 132         {
 133             return ret;
 134         }
 135         in_data  = (real*)in_data+mfft->dist;
 136         out_data = (real*)out_data+mfft->dist;
 137     }
 138     return 0;
 139 }
 140
 141 int
 142 gmx_fft_many_1d_real     (gmx_fft_t                  fft,
 143                           enum gmx_fft_direction     dir,
 144                           void *                     in_data,
 145                           void *                     out_data)
 146 {
 147     gmx_many_fft_t mfft = (gmx_many_fft_t)fft;
 148     int            i, ret;
 149     for (i = 0; i < mfft->howmany; i++)
 150     {
 151         ret = gmx_fft_1d_real(mfft->fft, dir, in_data, out_data);
 152         if (ret != 0)
 153         {
 154             return ret;
 155         }
 156         in_data  = (real*)in_data+mfft->dist;
 157         out_data = (real*)out_data+mfft->dist;
 158     }
 159     return 0;
 160 }
 161
 162
 163 void
 164 gmx_many_fft_destroy(gmx_fft_t    fft)
 165 {
 166     gmx_many_fft_t mfft = (gmx_many_fft_t)fft;
 167     if (mfft != nullptr)
 168     {
 169         if (mfft->fft != nullptr)
 170         {
 171             gmx_fft_destroy(mfft->fft);
 172         }
 173         free(mfft);
 174     }
 175 }
 176
 177 #endif //not GMX_FFT_FFTW3
 178
 179 int gmx_fft_transpose_2d(t_complex *          in_data,
 180                          t_complex *          out_data,
 181                          int                  nx,
 182                          int                  ny)
 183 {
 184     int        i, j, k, im, n, ncount;
 185     bool       done1, done2;
 186     int        i1, i1c, i2, i2c, kmi, max;
 187
 188     t_complex  tmp1, tmp2, tmp3;
 189     t_complex *data;
 190
 191     /* Use 500 bytes on stack to indicate moves.
 192      * This is just for optimization, it does not limit any dimensions.
 193      */
 194     char          move[500];
 195     int           nmove = 500;
 196
 197     if (nx < 2 || ny < 2)
 198     {
 199         if (in_data != out_data)
 200         {
 201             memcpy(out_data, in_data, sizeof(t_complex)*nx*ny);
 202         }
 203         return 0;
 204     }
 205
 206     /* Out-of-place transposes are easy */
 207     if (in_data != out_data)
 208     {
 209         for (i = 0; i < nx; i++)
 210         {
 211             for (j = 0; j < ny; j++)
 212             {
 213                 out_data[j*nx+i].re = in_data[i*ny+j].re;
 214                 out_data[j*nx+i].im = in_data[i*ny+j].im;
 215             }
 216         }
 217         return 0;
 218     }
 219
 220     /* In-place transform. in_data=out_data=data */
 221     data = in_data;
 222
 223     if (nx == ny)
 224     {
 225         /* trivial case, just swap elements */
 226         for (i = 0; i < nx; i++)
 227         {
 228             for (j = i+1; j < nx; j++)
 229             {
 230                 tmp1.re         = data[i*nx+j].re;
 231                 tmp1.im         = data[i*nx+j].im;
 232                 data[i*nx+j].re = data[j*nx+i].re;
 233                 data[i*nx+j].im = data[j*nx+i].im;
 234                 data[j*nx+i].re = tmp1.re;
 235                 data[j*nx+i].im = tmp1.im;
 236             }
 237         }
 238         return 0;
 239     }
 240
 241     for (i = 0; i < nmove; i++)
 242     {
 243         move[i] = 0;
 244     }
 245
 246     ncount = 2;
 247
 248     if (nx > 2 && ny > 2)
 249     {
 250         i = nx-1;
 251         j = ny-1;
 252         do
 253         {
 254             k = i % j;
 255             i = j;
 256             j = k;
 257         }
 258         while (k);
 259         ncount += i-1;
 260     }
 261
 262     n  = nx*ny;
 263     k  = n - 1;
 264     i  = 1;
 265     im = ny;
 266
 267     done1 = false;
 268     do
 269     {
 270         i1      = i;
 271         kmi     = k-i;
 272         tmp1.re = data[i1].re;
 273         tmp1.im = data[i1].im;
 274         i1c     = kmi;
 275         tmp2.re = data[i1c].re;
 276         tmp2.im = data[i1c].im;
 277
 278         done2 = false;
 279         do
 280         {
 281             i2  = ny*i1-k*(i1/nx);
 282             i2c = k-i2;
 283             if (i1 < nmove)
 284             {
 285                 move[i1] = 1;
 286             }
 287             if (i1c < nmove)
 288             {
 289                 move[i1c] = 1;
 290             }
 291             ncount += 2;
 292             if (i2 == i)
 293             {
 294                 done2 = true;
 295             }
 296             else if (i2 == kmi)
 297             {
 298                 tmp3.re = tmp1.re;
 299                 tmp3.im = tmp1.im;
 300                 tmp1.re = tmp2.re;
 301                 tmp1.im = tmp2.im;
 302                 tmp2.re = tmp3.re;
 303                 tmp2.im = tmp3.im;
 304                 done2   = true;
 305             }
 306             else
 307             {
 308                 data[i1].re  = data[i2].re;
 309                 data[i1].im  = data[i2].im;
 310                 data[i1c].re = data[i2c].re;
 311                 data[i1c].im = data[i2c].im;
 312                 i1           = i2;
 313                 i1c          = i2c;
 314             }
 315         }
 316         while (!done2);
 317
 318         data[i1].re  = tmp1.re;
 319         data[i1].im  = tmp1.im;
 320         data[i1c].re = tmp2.re;
 321         data[i1c].im = tmp2.im;
 322
 323         if (ncount >= n)
 324         {
 325             done1 = true;
 326         }
 327         else
 328         {
 329             done2 = false;
 330             do
 331             {
 332                 max = k-i;
 333                 i++;
 334                 im += ny;
 335                 if (im > k)
 336                 {
 337                     im -= k;
 338                 }
 339                 i2 = im;
 340                 if (i != i2)
 341                 {
 342                     if (i >= nmove)
 343                     {
 344                         while (i2 > i && i2 < max)
 345                         {
 346                             i1 = i2;
 347                             i2 = ny*i1-k*(i1/nx);
 348                         }
 349                         if (i2 == i)
 350                         {
 351                             done2 = true;
 352                         }
 353                     }
 354                     else if (!move[i])
 355                     {
 356                         done2 = true;
 357                     }
 358                 }
 359             }
 360             while (!done2);
 361         }
 362     }
 363     while (!done1);
 364
 365     return 0;
 366 }