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