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