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