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, 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 "gromacs/fft/fft.h"
  39
  40 #include "config.h"
  41
  42 #include <errno.h>
  43 #include <stdio.h>
  44 #include <stdlib.h>
  45 #include <string.h>
  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 #ifndef 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 == NULL)
  74     {
  75         gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
  76         return EINVAL;
  77     }
  78     *pfft = NULL;
  79
  80     if ( (fft = (gmx_many_fft_t)malloc(sizeof(struct gmx_many_fft))) == NULL)
  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 == NULL)
 101     {
 102         gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
 103         return EINVAL;
 104     }
 105     *pfft = NULL;
 106
 107     if ( (fft = (gmx_many_fft_t)malloc(sizeof(struct gmx_many_fft))) == NULL)
 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 != NULL)
 168     {
 169         if (mfft->fft != NULL)
 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, done1, done2;
 185     int        i1, i1c, i2, i2c, kmi, max;
 186
 187     t_complex  tmp1, tmp2, tmp3;
 188     t_complex *data;
 189
 190     /* Use 500 bytes on stack to indicate moves.
 191      * This is just for optimization, it does not limit any dimensions.
 192      */
 193     char          move[500];
 194     int           nmove = 500;
 195
 196     if (nx < 2 || ny < 2)
 197     {
 198         if (in_data != out_data)
 199         {
 200             memcpy(out_data, in_data, sizeof(t_complex)*nx*ny);
 201         }
 202         return 0;
 203     }
 204
 205     /* Out-of-place transposes are easy */
 206     if (in_data != out_data)
 207     {
 208         for (i = 0; i < nx; i++)
 209         {
 210             for (j = 0; j < ny; j++)
 211             {
 212                 out_data[j*nx+i].re = in_data[i*ny+j].re;
 213                 out_data[j*nx+i].im = in_data[i*ny+j].im;
 214             }
 215         }
 216         return 0;
 217     }
 218
 219     /* In-place transform. in_data=out_data=data */
 220     data = in_data;
 221
 222     if (nx == ny)
 223     {
 224         /* trivial case, just swap elements */
 225         for (i = 0; i < nx; i++)
 226         {
 227             for (j = i+1; j < nx; j++)
 228             {
 229                 tmp1.re         = data[i*nx+j].re;
 230                 tmp1.im         = data[i*nx+j].im;
 231                 data[i*nx+j].re = data[j*nx+i].re;
 232                 data[i*nx+j].im = data[j*nx+i].im;
 233                 data[j*nx+i].re = tmp1.re;
 234                 data[j*nx+i].im = tmp1.im;
 235             }
 236         }
 237         return 0;
 238     }
 239
 240     for (i = 0; i < nmove; i++)
 241     {
 242         move[i] = 0;
 243     }
 244
 245     ncount = 2;
 246
 247     if (nx > 2 && ny > 2)
 248     {
 249         i = nx-1;
 250         j = ny-1;
 251         do
 252         {
 253             k = i % j;
 254             i = j;
 255             j = k;
 256         }
 257         while (k);
 258         ncount += i-1;
 259     }
 260
 261     n  = nx*ny;
 262     k  = n - 1;
 263     i  = 1;
 264     im = ny;
 265
 266     done1 = 0;
 267     do
 268     {
 269         i1      = i;
 270         kmi     = k-i;
 271         tmp1.re = data[i1].re;
 272         tmp1.im = data[i1].im;
 273         i1c     = kmi;
 274         tmp2.re = data[i1c].re;
 275         tmp2.im = data[i1c].im;
 276
 277         done2 = 0;
 278         do
 279         {
 280             i2  = ny*i1-k*(i1/nx);
 281             i2c = k-i2;
 282             if (i1 < nmove)
 283             {
 284                 move[i1] = 1;
 285             }
 286             if (i1c < nmove)
 287             {
 288                 move[i1c] = 1;
 289             }
 290             ncount += 2;
 291             if (i2 == i)
 292             {
 293                 done2 = 1;
 294             }
 295             else if (i2 == kmi)
 296             {
 297                 tmp3.re = tmp1.re;
 298                 tmp3.im = tmp1.im;
 299                 tmp1.re = tmp2.re;
 300                 tmp1.im = tmp2.im;
 301                 tmp2.re = tmp3.re;
 302                 tmp2.im = tmp3.im;
 303                 done2   = 1;
 304             }
 305             else
 306             {
 307                 data[i1].re  = data[i2].re;
 308                 data[i1].im  = data[i2].im;
 309                 data[i1c].re = data[i2c].re;
 310                 data[i1c].im = data[i2c].im;
 311                 i1           = i2;
 312                 i1c          = i2c;
 313             }
 314         }
 315         while (!done2);
 316
 317         data[i1].re  = tmp1.re;
 318         data[i1].im  = tmp1.im;
 319         data[i1c].re = tmp2.re;
 320         data[i1c].im = tmp2.im;
 321
 322         if (ncount >= n)
 323         {
 324             done1 = 1;
 325         }
 326         else
 327         {
 328             done2 = 0;
 329             do
 330             {
 331                 max = k-i;
 332                 i++;
 333                 im += ny;
 334                 if (im > k)
 335                 {
 336                     im -= k;
 337                 }
 338                 i2 = im;
 339                 if (i != i2)
 340                 {
 341                     if (i >= nmove)
 342                     {
 343                         while (i2 > i && i2 < max)
 344                         {
 345                             i1 = i2;
 346                             i2 = ny*i1-k*(i1/nx);
 347                         }
 348                         if (i2 == i)
 349                         {
 350                             done2 = 1;
 351                         }
 352                     }
 353                     else if (!move[i])
 354                     {
 355                         done2 = 1;
 356                     }
 357                 }
 358             }
 359             while (!done2);
 360         }
 361     }
 362     while (!done1);
 363
 364     return 0;
 365 }