src/gromacs/fft/fft_fftpack.c

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