src/gromacs/fft/fft_fftpack.cpp

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