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