src/gromacs/fft/fft_fftpack.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2003 David van der Spoel, Erik Lindahl, 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 <math.h>
  44 #include <stdlib.h>
  45 #include <string.h>
  46
  47 #include "gromacs/utility/fatalerror.h"
  48 #include "gromacs/utility/real.h"
  49
  50 #include "external/fftpack/fftpack.h"
  51
  52 /*! \internal
  53  * \brief
  54  * Contents of the FFTPACK fft datatype.
  55  *
  56  * Note that this is one of several possible implementations of gmx_fft_t.
  57  *
  58  *  FFTPACK only does 1d transforms, so we use a pointers to another fft for
  59  *  the transform in the next dimension.
  60  * Thus, a 3d-structure contains a pointer to a 2d one, which in turns contains
  61  * a pointer to a 1d. The 1d structure has next==NULL.
  62  */
  63 #ifdef DOXYGEN
  64 struct gmx_fft_fftpack
  65 #else
  66 struct gmx_fft
  67 #endif
  68 {
  69     int             ndim;     /**< Dimensions, including our subdimensions.  */
  70     int             n;        /**< Number of points in this dimension.       */
  71     int             ifac[15]; /**< 15 bytes needed for cfft and rfft         */
  72     struct gmx_fft *next;     /**< Pointer to next dimension, or NULL.       */
  73     real  *         work;     /**< 1st 4n reserved for cfft, 1st 2n for rfft */
  74 };
  75
  76 #include <math.h>
  77 #include <stdio.h>
  78
  79
  80 int
  81 gmx_fft_init_1d(gmx_fft_t *        pfft,
  82                 int                nx,
  83                 int gmx_unused     flags)
  84 {
  85     gmx_fft_t    fft;
  86
  87     if (pfft == NULL)
  88     {
  89         gmx_fatal(FARGS, "Invalid FFT opaque type pointer.");
  90         return EINVAL;
  91     }
  92     *pfft = NULL;
  93
  94     if ( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
  95     {
  96         return ENOMEM;
  97     }
  98
  99     fft->next = NULL;
 100     fft->n    = nx;
 101
 102     /* Need 4*n storage for 1D complex FFT */
 103     if ( (fft->work = (real *)malloc(sizeof(real)*(4*nx))) == NULL)
 104     {
 105         free(fft);
 106         return ENOMEM;
 107     }
 108
 109     if (fft->n > 1)
 110     {
 111         fftpack_cffti1(nx, fft->work, fft->ifac);
 112     }
 113
 114     *pfft = fft;
 115     return 0;
 116 };
 117
 118
 119
 120 int
 121 gmx_fft_init_1d_real(gmx_fft_t *        pfft,
 122                      int                nx,
 123                      int gmx_unused     flags)
 124 {
 125     gmx_fft_t    fft;
 126
 127     if (pfft == NULL)
 128     {
 129         gmx_fatal(FARGS, "Invalid FFT opaque type pointer.");
 130         return EINVAL;
 131     }
 132     *pfft = NULL;
 133
 134     if ( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
 135     {
 136         return ENOMEM;
 137     }
 138
 139     fft->next = NULL;
 140     fft->n    = nx;
 141
 142     /* Need 2*n storage for 1D real FFT */
 143     if ((fft->work = (real *)malloc(sizeof(real)*(2*nx))) == NULL)
 144     {
 145         free(fft);
 146         return ENOMEM;
 147     }
 148
 149     if (fft->n > 1)
 150     {
 151         fftpack_rffti1(nx, fft->work, fft->ifac);
 152     }
 153
 154     *pfft = fft;
 155     return 0;
 156 }
 157
 158 int
 159 gmx_fft_init_2d_real(gmx_fft_t *        pfft,
 160                      int                nx,
 161                      int                ny,
 162                      int                flags)
 163 {
 164     gmx_fft_t     fft;
 165     int           nyc = (ny/2 + 1);
 166     int           rc;
 167
 168     if (pfft == NULL)
 169     {
 170         gmx_fatal(FARGS, "Invalid FFT opaque type pointer.");
 171         return EINVAL;
 172     }
 173     *pfft = NULL;
 174
 175     /* Create the X transform */
 176     if ( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
 177     {
 178         return ENOMEM;
 179     }
 180
 181     fft->n    = nx;
 182
 183     /* Need 4*nx storage for 1D complex FFT, and another
 184      * 2*nx*nyc elements for complex-to-real storage in our high-level routine.
 185      */
 186     if ( (fft->work = (real *)malloc(sizeof(real)*(4*nx+2*nx*nyc))) == NULL)
 187     {
 188         free(fft);
 189         return ENOMEM;
 190     }
 191     fftpack_cffti1(nx, fft->work, fft->ifac);
 192
 193     /* Create real Y transform as a link from X */
 194     if ( (rc = gmx_fft_init_1d_real(&(fft->next), ny, flags)) != 0)
 195     {
 196         free(fft);
 197         return rc;
 198     }
 199
 200     *pfft = fft;
 201     return 0;
 202 }
 203
 204
 205 int
 206 gmx_fft_1d               (gmx_fft_t                  fft,
 207                           enum gmx_fft_direction     dir,
 208                           void *                     in_data,
 209                           void *                     out_data)
 210 {
 211     int             i, n;
 212     real       *    p1;
 213     real       *    p2;
 214
 215     n = fft->n;
 216
 217     if (n == 1)
 218     {
 219         p1    = (real *)in_data;
 220         p2    = (real *)out_data;
 221         p2[0] = p1[0];
 222         p2[1] = p1[1];
 223     }
 224
 225     /* FFTPACK only does in-place transforms, so emulate out-of-place
 226      * by copying data to the output array first.
 227      */
 228     if (in_data != out_data)
 229     {
 230         p1 = (real *)in_data;
 231         p2 = (real *)out_data;
 232
 233         /* n complex = 2*n real elements */
 234         for (i = 0; i < 2*n; i++)
 235         {
 236             p2[i] = p1[i];
 237         }
 238     }
 239
 240     /* Elements 0   .. 2*n-1 in work are used for ffac values,
 241      * Elements 2*n .. 4*n-1 are internal FFTPACK work space.
 242      */
 243
 244     if (dir == GMX_FFT_FORWARD)
 245     {
 246         fftpack_cfftf1(n, (real *)out_data, fft->work+2*n, fft->work, fft->ifac, -1);
 247     }
 248     else if (dir == GMX_FFT_BACKWARD)
 249     {
 250         fftpack_cfftf1(n, (real *)out_data, fft->work+2*n, fft->work, fft->ifac, 1);
 251     }
 252     else
 253     {
 254         gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
 255         return EINVAL;
 256     }
 257
 258     return 0;
 259 }
 260
 261
 262
 263 int
 264 gmx_fft_1d_real          (gmx_fft_t                  fft,
 265                           enum gmx_fft_direction     dir,
 266                           void *                     in_data,
 267                           void *                     out_data)
 268 {
 269     int           i, n;
 270     real       *  p1;
 271     real       *  p2;
 272
 273     n = fft->n;
 274
 275     if (n == 1)
 276     {
 277         p1    = (real *)in_data;
 278         p2    = (real *)out_data;
 279         p2[0] = p1[0];
 280         if (dir == GMX_FFT_REAL_TO_COMPLEX)
 281         {
 282             p2[1] = 0.0;
 283         }
 284     }
 285
 286     if (dir == GMX_FFT_REAL_TO_COMPLEX)
 287     {
 288         /* FFTPACK only does in-place transforms, so emulate out-of-place
 289          * by copying data to the output array first. This works fine, since
 290          * the complex array must be larger than the real.
 291          */
 292         if (in_data != out_data)
 293         {
 294             p1 = (real *)in_data;
 295             p2 = (real *)out_data;
 296
 297             for (i = 0; i < 2*(n/2+1); i++)
 298             {
 299                 p2[i] = p1[i];
 300             }
 301         }
 302
 303         /* Elements 0 ..   n-1 in work are used for ffac values,
 304          * Elements n .. 2*n-1 are internal FFTPACK work space.
 305          */
 306         fftpack_rfftf1(n, (real *)out_data, fft->work+n, fft->work, fft->ifac);
 307
 308         /*
 309          * FFTPACK has a slightly more compact storage than we, time to
 310          * convert it: ove most of the array one step up to make room for
 311          * zero imaginary parts.
 312          */
 313         p2 = (real *)out_data;
 314         for (i = n-1; i > 0; i--)
 315         {
 316             p2[i+1] = p2[i];
 317         }
 318         /* imaginary zero freq. */
 319         p2[1] = 0;
 320
 321         /* Is n even? */
 322         if ( (n & 0x1) == 0)
 323         {
 324             p2[n+1] = 0;
 325         }
 326
 327     }
 328     else if (dir == GMX_FFT_COMPLEX_TO_REAL)
 329     {
 330         /* FFTPACK only does in-place transforms, and we cannot just copy
 331          * input to output first here since our real array is smaller than
 332          * the complex one. However, since the FFTPACK complex storage format
 333          * is more compact than ours (2 reals) it will fit, so compact it
 334          * and copy on-the-fly to the output array.
 335          */
 336         p1 = (real *) in_data;
 337         p2 = (real *)out_data;
 338
 339         p2[0] = p1[0];
 340         for (i = 1; i < n; i++)
 341         {
 342             p2[i] = p1[i+1];
 343         }
 344         fftpack_rfftb1(n, (real *)out_data, fft->work+n, fft->work, fft->ifac);
 345     }
 346     else
 347     {
 348         gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
 349         return EINVAL;
 350     }
 351
 352     return 0;
 353 }
 354
 355
 356 int
 357 gmx_fft_2d_real          (gmx_fft_t                  fft,
 358                           enum gmx_fft_direction     dir,
 359                           void *                     in_data,
 360                           void *                     out_data)
 361 {
 362     int                i, j, nx, ny, nyc;
 363     t_complex     *    data;
 364     real       *       work;
 365     real       *       p1;
 366     real       *       p2;
 367
 368     nx = fft->n;
 369     ny = fft->next->n;
 370     /* Number of complex elements in y direction */
 371     nyc = (ny/2+1);
 372
 373     work = fft->work+4*nx;
 374
 375     if (dir == GMX_FFT_REAL_TO_COMPLEX)
 376     {
 377         /* If we are doing an in-place transform the 2D array is already
 378          * properly padded by the user, and we are all set.
 379          *
 380          * For out-of-place there is no array padding, but FFTPACK only
 381          * does in-place FFTs internally, so we need to start by copying
 382          * data from the input to the padded (larger) output array.
 383          */
 384         if (in_data != out_data)
 385         {
 386             p1 = (real *)in_data;
 387             p2 = (real *)out_data;
 388
 389             for (i = 0; i < nx; i++)
 390             {
 391                 for (j = 0; j < ny; j++)
 392                 {
 393                     p2[i*nyc*2+j] = p1[i*ny+j];
 394                 }
 395             }
 396         }
 397         data = (t_complex *)out_data;
 398
 399         /* y real-to-complex FFTs */
 400         for (i = 0; i < nx; i++)
 401         {
 402             gmx_fft_1d_real(fft->next, GMX_FFT_REAL_TO_COMPLEX, data+i*nyc, data+i*nyc);
 403         }
 404
 405         /* Transform to get X data in place */
 406         gmx_fft_transpose_2d(data, data, nx, nyc);
 407
 408         /* Complex-to-complex X FFTs */
 409         for (i = 0; i < nyc; i++)
 410         {
 411             gmx_fft_1d(fft, GMX_FFT_FORWARD, data+i*nx, data+i*nx);
 412         }
 413
 414         /* Transpose back */
 415         gmx_fft_transpose_2d(data, data, nyc, nx);
 416
 417     }
 418     else if (dir == GMX_FFT_COMPLEX_TO_REAL)
 419     {
 420         /* An in-place complex-to-real transform is straightforward,
 421          * since the output array must be large enough for the padding to fit.
 422          *
 423          * For out-of-place complex-to-real transforms we cannot just copy
 424          * data to the output array, since it is smaller than the input.
 425          * In this case there's nothing to do but employing temporary work data,
 426          * starting at work+4*nx and using nx*nyc*2 elements.
 427          */
 428         if (in_data != out_data)
 429         {
 430             memcpy(work, in_data, sizeof(t_complex)*nx*nyc);
 431             data = (t_complex *)work;
 432         }
 433         else
 434         {
 435             /* in-place */
 436             data = (t_complex *)out_data;
 437         }
 438
 439         /* Transpose to get X arrays */
 440         gmx_fft_transpose_2d(data, data, nx, nyc);
 441
 442         /* Do X iFFTs */
 443         for (i = 0; i < nyc; i++)
 444         {
 445             gmx_fft_1d(fft, GMX_FFT_BACKWARD, data+i*nx, data+i*nx);
 446         }
 447
 448         /* Transpose to get Y arrays */
 449         gmx_fft_transpose_2d(data, data, nyc, nx);
 450
 451         /* Do Y iFFTs */
 452         for (i = 0; i < nx; i++)
 453         {
 454             gmx_fft_1d_real(fft->next, GMX_FFT_COMPLEX_TO_REAL, data+i*nyc, data+i*nyc);
 455         }
 456
 457         if (in_data != out_data)
 458         {
 459             /* Output (pointed to by data) is now in padded format.
 460              * Pack it into out_data if we were doing an out-of-place transform.
 461              */
 462             p1 = (real *)data;
 463             p2 = (real *)out_data;
 464
 465             for (i = 0; i < nx; i++)
 466             {
 467                 for (j = 0; j < ny; j++)
 468                 {
 469                     p2[i*ny+j] = p1[i*nyc*2+j];
 470                 }
 471             }
 472         }
 473     }
 474     else
 475     {
 476         gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
 477         return EINVAL;
 478     }
 479
 480     return 0;
 481 }
 482
 483 void
 484 gmx_fft_destroy(gmx_fft_t      fft)
 485 {
 486     if (fft != NULL)
 487     {
 488         free(fft->work);
 489         if (fft->next != NULL)
 490         {
 491             gmx_fft_destroy(fft->next);
 492         }
 493         free(fft);
 494     }
 495 }
 496
 497 void gmx_fft_cleanup()
 498 {
 499 }
 500
 501 const char *gmx_fft_get_version_info()
 502 {
 503     return "fftpack (built-in)";
 504 }