src/gromacs/fft/fft.c

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