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