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.
   6  * Copyright (c) 2019,2020, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "fft.h"
  40
  41 #include "config.h"
  42
  43 #include <cerrno>
  44 #include <cstdio>
  45 #include <cstdlib>
  46 #include <cstring>
  47
  48 #include "gromacs/math/gmxcomplex.h"
  49 #include "gromacs/utility/fatalerror.h"
  50 #include "gromacs/utility/real.h"
  51
  52 /* This file contains common fft utility functions, but not
  53  * the actual transform implementations. Check the
  54  * files like fft_fftw3.c or fft_mkl.c for that.
  55  */
  56
  57 #if !GMX_FFT_FFTW3 && !GMX_FFT_ARMPL_FFTW3
  58
  59 struct gmx_many_fft
  60 {
  61     int       howmany;
  62     int       dist;
  63     gmx_fft_t fft;
  64 };
  65
  66 typedef struct gmx_many_fft* gmx_many_fft_t;
  67
  68 int gmx_fft_init_many_1d(gmx_fft_t* pfft, int nx, int howmany, gmx_fft_flag flags)
  69 {
  70     gmx_many_fft_t fft;
  71     if (pfft == nullptr)
  72     {
  73         gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
  74         return EINVAL;
  75     }
  76     *pfft = nullptr;
  77
  78     if ((fft = static_cast<gmx_many_fft_t>(malloc(sizeof(struct gmx_many_fft)))) == nullptr)
  79     {
  80         return ENOMEM;
  81     }
  82
  83     gmx_fft_init_1d(&fft->fft, nx, flags);
  84     fft->howmany = howmany;
  85     fft->dist    = 2 * nx;
  86
  87     *pfft = reinterpret_cast<gmx_fft_t>(fft);
  88     return 0;
  89 }
  90
  91 int gmx_fft_init_many_1d_real(gmx_fft_t* pfft, int nx, int howmany, gmx_fft_flag flags)
  92 {
  93     gmx_many_fft_t fft;
  94     if (pfft == nullptr)
  95     {
  96         gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
  97         return EINVAL;
  98     }
  99     *pfft = nullptr;
 100
 101     if ((fft = static_cast<gmx_many_fft_t>(malloc(sizeof(struct gmx_many_fft)))) == nullptr)
 102     {
 103         return ENOMEM;
 104     }
 105
 106     gmx_fft_init_1d_real(&fft->fft, nx, flags);
 107     fft->howmany = howmany;
 108     fft->dist    = 2 * (nx / 2 + 1);
 109
 110     *pfft = reinterpret_cast<gmx_fft_t>(fft);
 111     return 0;
 112 }
 113
 114 int gmx_fft_many_1d(gmx_fft_t fft, enum gmx_fft_direction dir, void* in_data, void* out_data)
 115 {
 116     gmx_many_fft_t mfft = reinterpret_cast<gmx_many_fft_t>(fft);
 117     int            i, ret;
 118     for (i = 0; i < mfft->howmany; i++)
 119     {
 120         ret = gmx_fft_1d(mfft->fft, dir, in_data, out_data);
 121         if (ret != 0)
 122         {
 123             return ret;
 124         }
 125         in_data  = static_cast<real*>(in_data) + mfft->dist;
 126         out_data = static_cast<real*>(out_data) + mfft->dist;
 127     }
 128     return 0;
 129 }
 130
 131 int gmx_fft_many_1d_real(gmx_fft_t fft, enum gmx_fft_direction dir, void* in_data, void* out_data)
 132 {
 133     gmx_many_fft_t mfft = reinterpret_cast<gmx_many_fft_t>(fft);
 134     int            i, ret;
 135     for (i = 0; i < mfft->howmany; i++)
 136     {
 137         ret = gmx_fft_1d_real(mfft->fft, dir, in_data, out_data);
 138         if (ret != 0)
 139         {
 140             return ret;
 141         }
 142         in_data  = static_cast<real*>(in_data) + mfft->dist;
 143         out_data = static_cast<real*>(out_data) + mfft->dist;
 144     }
 145     return 0;
 146 }
 147
 148
 149 void gmx_many_fft_destroy(gmx_fft_t fft)
 150 {
 151     gmx_many_fft_t mfft = reinterpret_cast<gmx_many_fft_t>(fft);
 152     if (mfft != nullptr)
 153     {
 154         if (mfft->fft != nullptr)
 155         {
 156             gmx_fft_destroy(mfft->fft);
 157         }
 158         free(mfft);
 159     }
 160 }
 161
 162 #endif // not GMX_FFT_FFTW3
 163
 164 int gmx_fft_transpose_2d(t_complex* in_data, t_complex* out_data, int nx, int ny)
 165 {
 166     int  i, j, k, im, n, ncount;
 167     bool done1, done2;
 168     int  i1, i1c, i2, i2c, kmi, max;
 169
 170     t_complex  tmp1, tmp2, tmp3;
 171     t_complex* data;
 172
 173     /* Use 500 bytes on stack to indicate moves.
 174      * This is just for optimization, it does not limit any dimensions.
 175      */
 176     char move[500];
 177     int  nmove = 500;
 178
 179     if (nx < 2 || ny < 2)
 180     {
 181         if (in_data != out_data)
 182         {
 183             memcpy(out_data, in_data, sizeof(t_complex) * nx * ny);
 184         }
 185         return 0;
 186     }
 187
 188     /* Out-of-place transposes are easy */
 189     if (in_data != out_data)
 190     {
 191         for (i = 0; i < nx; i++)
 192         {
 193             for (j = 0; j < ny; j++)
 194             {
 195                 out_data[j * nx + i].re = in_data[i * ny + j].re;
 196                 out_data[j * nx + i].im = in_data[i * ny + j].im;
 197             }
 198         }
 199         return 0;
 200     }
 201
 202     /* In-place transform. in_data=out_data=data */
 203     data = in_data;
 204
 205     if (nx == ny)
 206     {
 207         /* trivial case, just swap elements */
 208         for (i = 0; i < nx; i++)
 209         {
 210             for (j = i + 1; j < nx; j++)
 211             {
 212                 tmp1.re             = data[i * nx + j].re;
 213                 tmp1.im             = data[i * nx + j].im;
 214                 data[i * nx + j].re = data[j * nx + i].re;
 215                 data[i * nx + j].im = data[j * nx + i].im;
 216                 data[j * nx + i].re = tmp1.re;
 217                 data[j * nx + i].im = tmp1.im;
 218             }
 219         }
 220         return 0;
 221     }
 222
 223     for (i = 0; i < nmove; i++)
 224     {
 225         move[i] = 0;
 226     }
 227
 228     ncount = 2;
 229
 230     if (nx > 2 && ny > 2)
 231     {
 232         i = nx - 1;
 233         j = ny - 1;
 234         do
 235         {
 236             k = i % j;
 237             i = j;
 238             j = k;
 239         } while (k != 0);
 240         ncount += i - 1;
 241     }
 242
 243     n  = nx * ny;
 244     k  = n - 1;
 245     i  = 1;
 246     im = ny;
 247
 248     done1 = false;
 249     do
 250     {
 251         i1      = i;
 252         kmi     = k - i;
 253         tmp1.re = data[i1].re;
 254         tmp1.im = data[i1].im;
 255         i1c     = kmi;
 256         tmp2.re = data[i1c].re;
 257         tmp2.im = data[i1c].im;
 258
 259         done2 = false;
 260         do
 261         {
 262             i2  = ny * i1 - k * (i1 / nx);
 263             i2c = k - i2;
 264             if (i1 < nmove)
 265             {
 266                 move[i1] = 1;
 267             }
 268             if (i1c < nmove)
 269             {
 270                 move[i1c] = 1;
 271             }
 272             ncount += 2;
 273             if (i2 == i)
 274             {
 275                 done2 = true;
 276             }
 277             else if (i2 == kmi)
 278             {
 279                 tmp3.re = tmp1.re;
 280                 tmp3.im = tmp1.im;
 281                 tmp1.re = tmp2.re;
 282                 tmp1.im = tmp2.im;
 283                 tmp2.re = tmp3.re;
 284                 tmp2.im = tmp3.im;
 285                 done2   = true;
 286             }
 287             else
 288             {
 289                 data[i1].re  = data[i2].re;
 290                 data[i1].im  = data[i2].im;
 291                 data[i1c].re = data[i2c].re;
 292                 data[i1c].im = data[i2c].im;
 293                 i1           = i2;
 294                 i1c          = i2c;
 295             }
 296         } while (!done2);
 297
 298         data[i1].re  = tmp1.re;
 299         data[i1].im  = tmp1.im;
 300         data[i1c].re = tmp2.re;
 301         data[i1c].im = tmp2.im;
 302
 303         if (ncount >= n)
 304         {
 305             done1 = true;
 306         }
 307         else
 308         {
 309             done2 = false;
 310             do
 311             {
 312                 max = k - i;
 313                 i++;
 314                 im += ny;
 315                 if (im > k)
 316                 {
 317                     im -= k;
 318                 }
 319                 i2 = im;
 320                 if (i != i2)
 321                 {
 322                     if (i >= nmove)
 323                     {
 324                         while (i2 > i && i2 < max)
 325                         {
 326                             i1 = i2;
 327                             i2 = ny * i1 - k * (i1 / nx);
 328                         }
 329                         if (i2 == i)
 330                         {
 331                             done2 = true;
 332                         }
 333                     }
 334                     else if (!move[i])
 335                     {
 336                         done2 = true;
 337                     }
 338                 }
 339             } while (!done2);
 340         }
 341     } while (!done1);
 342
 343     return 0;
 344 }