src/mdlib/gmx_fft_fftpack.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Gromacs 4.0                         Copyright (c) 1991-2003
   5  * David van der Spoel, Erik Lindahl, University of Groningen.
   6  * Copyright (c) 2012, by the GROMACS development team, led by
   7  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   8  * others, as listed in the AUTHORS file in the top-level source
   9  * directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #ifdef GMX_FFT_FFTPACK
  42
  43 #include <math.h>
  44 #include <stdlib.h>
  45 #include <string.h>
  46 #include <errno.h>
  47
  48
  49 #include "gmx_fft.h"
  50 #include "gmx_fatal.h"
  51 #include "fftpack.h"
  52
  53 /** Contents of the FFTPACK fft datatype.
  54  *
  55  *  FFTPACK only does 1d transforms, so we use a pointers to another fft for
  56  *  the transform in the next dimension.
  57  * Thus, a 3d-structure contains a pointer to a 2d one, which in turns contains
  58  * a pointer to a 1d. The 1d structure has next==NULL.
  59  */
  60 struct gmx_fft
  61 {
  62     int            ndim;     /**< Dimensions, including our subdimensions.  */
  63     int            n;        /**< Number of points in this dimension.       */
  64     int            ifac[15]; /**< 15 bytes needed for cfft and rfft         */
  65     struct gmx_fft *next;    /**< Pointer to next dimension, or NULL.       */
  66     real *         work;     /**< 1st 4n reserved for cfft, 1st 2n for rfft */
  67 };
  68
  69 #include <math.h>
  70 #include <stdio.h>
  71
  72
  73 int
  74 gmx_fft_init_1d(gmx_fft_t *        pfft,
  75                 int                nx,
  76                 int                flags)
  77 {
  78     gmx_fft_t    fft;
  79
  80     if(pfft==NULL)
  81     {
  82         gmx_fatal(FARGS,"Invalid FFT opaque type pointer.");
  83         return EINVAL;
  84     }
  85     *pfft = NULL;
  86
  87     if( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
  88     {
  89         return ENOMEM;
  90     }
  91
  92     fft->next = NULL;
  93     fft->n    = nx;
  94
  95     /* Need 4*n storage for 1D complex FFT */
  96     if( (fft->work = (real *)malloc(sizeof(real)*(4*nx))) == NULL)
  97     {
  98         free(fft);
  99         return ENOMEM;
 100     }
 101
 102     if(fft->n>1)
 103         fftpack_cffti1(nx,fft->work,fft->ifac);
 104
 105     *pfft = fft;
 106     return 0;
 107 };
 108
 109
 110
 111 int
 112 gmx_fft_init_1d_real(gmx_fft_t *        pfft,
 113                      int                nx,
 114                      int                flags)
 115 {
 116     gmx_fft_t    fft;
 117
 118     if(pfft==NULL)
 119     {
 120         gmx_fatal(FARGS,"Invalid FFT opaque type pointer.");
 121         return EINVAL;
 122     }
 123     *pfft = NULL;
 124
 125     if( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
 126     {
 127         return ENOMEM;
 128     }
 129
 130     fft->next = NULL;
 131     fft->n    = nx;
 132
 133     /* Need 2*n storage for 1D real FFT */
 134     if((fft->work = (real *)malloc(sizeof(real)*(2*nx)))==NULL)
 135     {
 136         free(fft);
 137         return ENOMEM;
 138     }
 139
 140     if(fft->n>1)
 141         fftpack_rffti1(nx,fft->work,fft->ifac);
 142
 143     *pfft = fft;
 144     return 0;
 145 }
 146
 147
 148
 149 int
 150 gmx_fft_init_2d(gmx_fft_t *        pfft,
 151                 int                nx,
 152                 int                ny,
 153                 int                flags)
 154 {
 155     gmx_fft_t     fft;
 156     int           rc;
 157
 158     if(pfft==NULL)
 159     {
 160         gmx_fatal(FARGS,"Invalid FFT opaque type pointer.");
 161         return EINVAL;
 162     }
 163     *pfft = NULL;
 164
 165     /* Create the X transform */
 166     if( (rc = gmx_fft_init_1d(&fft,nx,flags)) != 0)
 167     {
 168         return rc;
 169     }
 170
 171     /* Create Y transform as a link from X */
 172     if( (rc=gmx_fft_init_1d(&(fft->next),ny,flags)) != 0)
 173     {
 174         free(fft);
 175         return rc;
 176     }
 177
 178     *pfft = fft;
 179     return 0;
 180 };
 181
 182
 183 int
 184 gmx_fft_init_2d_real(gmx_fft_t *        pfft,
 185                      int                nx,
 186                      int                ny,
 187                      int                flags)
 188 {
 189     gmx_fft_t     fft;
 190     int           nyc = (ny/2 + 1);
 191     int           rc;
 192
 193     if(pfft==NULL)
 194     {
 195         gmx_fatal(FARGS,"Invalid FFT opaque type pointer.");
 196         return EINVAL;
 197     }
 198     *pfft = NULL;
 199
 200     /* Create the X transform */
 201     if( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
 202     {
 203         return ENOMEM;
 204     }
 205
 206     fft->n    = nx;
 207
 208     /* Need 4*nx storage for 1D complex FFT, and another
 209      * 2*nx*nyc elements for complex-to-real storage in our high-level routine.
 210      */
 211     if( (fft->work = (real *)malloc(sizeof(real)*(4*nx+2*nx*nyc))) == NULL)
 212     {
 213         free(fft);
 214         return ENOMEM;
 215     }
 216     fftpack_cffti1(nx,fft->work,fft->ifac);
 217
 218     /* Create real Y transform as a link from X */
 219     if( (rc=gmx_fft_init_1d_real(&(fft->next),ny,flags)) != 0)
 220     {
 221         free(fft);
 222         return rc;
 223     }
 224
 225     *pfft = fft;
 226     return 0;
 227 }
 228
 229
 230 int
 231 gmx_fft_init_3d(gmx_fft_t *        pfft,
 232                 int                nx,
 233                 int                ny,
 234                 int                nz,
 235                 int                flags)
 236 {
 237     gmx_fft_t     fft;
 238     int           rc;
 239
 240     if(pfft==NULL)
 241     {
 242         gmx_fatal(FARGS,"Invalid FFT opaque type pointer.");
 243         return EINVAL;
 244     }
 245     *pfft = NULL;
 246
 247     /* Create the X transform */
 248
 249     if( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
 250     {
 251         return ENOMEM;
 252     }
 253
 254     fft->n    = nx;
 255
 256     /* Need 4*nx storage for 1D complex FFT, and another
 257      * 2*nz elements for gmx_fft_transpose_2d_nelem() storage.
 258      */
 259     if( (fft->work = (real *)malloc(sizeof(real)*(4*nx+2*nz))) == NULL)
 260     {
 261         free(fft);
 262         return ENOMEM;
 263     }
 264
 265     fftpack_cffti1(nx,fft->work,fft->ifac);
 266
 267
 268     /* Create 2D Y/Z transforms as a link from X */
 269     if( (rc=gmx_fft_init_2d(&(fft->next),ny,nz,flags)) != 0)
 270     {
 271         free(fft);
 272         return rc;
 273     }
 274
 275     *pfft = fft;
 276     return 0;
 277 };
 278
 279
 280 int
 281 gmx_fft_init_3d_real(gmx_fft_t *        pfft,
 282                      int                nx,
 283                      int                ny,
 284                      int                nz,
 285                      int                flags)
 286 {
 287     gmx_fft_t     fft;
 288     int           nzc = (nz/2 + 1);
 289     int           rc;
 290
 291     if(pfft==NULL)
 292     {
 293         gmx_fatal(FARGS,"Invalid FFT opaque type pointer.");
 294         return EINVAL;
 295     }
 296     *pfft = NULL;
 297
 298     /* Create the X transform */
 299     if( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
 300     {
 301         return ENOMEM;
 302     }
 303
 304     fft->n    = nx;
 305
 306     /* Need 4*nx storage for 1D complex FFT, another
 307      * 2*nx*ny*nzc elements to copy the entire 3D matrix when
 308      * doing out-of-place complex-to-real FFTs, and finally
 309      * 2*nzc elements for transpose work space.
 310      */
 311     if( (fft->work = (real *)malloc(sizeof(real)*(4*nx+2*nx*ny*nzc+2*nzc))) == NULL)
 312     {
 313         free(fft);
 314         return ENOMEM;
 315     }
 316     fftpack_cffti1(nx,fft->work,fft->ifac);
 317
 318     /* Create 2D real Y/Z transform as a link from X */
 319     if( (rc=gmx_fft_init_2d_real(&(fft->next),ny,nz,flags)) != 0)
 320     {
 321         free(fft);
 322         return rc;
 323     }
 324
 325     *pfft = fft;
 326     return 0;
 327 }
 328
 329
 330 int
 331 gmx_fft_1d               (gmx_fft_t                  fft,
 332                           enum gmx_fft_direction     dir,
 333                           void *                     in_data,
 334                           void *                     out_data)
 335 {
 336     int             i,n;
 337     real *    p1;
 338     real *    p2;
 339
 340     n=fft->n;
 341
 342     if(n==1)
 343     {
 344         p1 = (real *)in_data;
 345         p2 = (real *)out_data;
 346         p2[0] = p1[0];
 347         p2[1] = p1[1];
 348     }
 349
 350     /* FFTPACK only does in-place transforms, so emulate out-of-place
 351      * by copying data to the output array first.
 352      */
 353     if( in_data != out_data )
 354     {
 355         p1 = (real *)in_data;
 356         p2 = (real *)out_data;
 357
 358         /* n complex = 2*n real elements */
 359         for(i=0;i<2*n;i++)
 360         {
 361             p2[i] = p1[i];
 362         }
 363     }
 364
 365     /* Elements 0   .. 2*n-1 in work are used for ffac values,
 366      * Elements 2*n .. 4*n-1 are internal FFTPACK work space.
 367      */
 368
 369     if(dir == GMX_FFT_FORWARD)
 370     {
 371         fftpack_cfftf1(n,(real *)out_data,fft->work+2*n,fft->work,fft->ifac, -1);
 372     }
 373     else if(dir == GMX_FFT_BACKWARD)
 374     {
 375         fftpack_cfftf1(n,(real *)out_data,fft->work+2*n,fft->work,fft->ifac, 1);
 376     }
 377     else
 378     {
 379         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
 380         return EINVAL;
 381     }
 382
 383     return 0;
 384 }
 385
 386
 387
 388 int
 389 gmx_fft_1d_real          (gmx_fft_t                  fft,
 390                           enum gmx_fft_direction     dir,
 391                           void *                     in_data,
 392                           void *                     out_data)
 393 {
 394     int           i,n;
 395     real *  p1;
 396     real *  p2;
 397
 398     n = fft->n;
 399
 400     if(n==1)
 401     {
 402         p1 = (real *)in_data;
 403         p2 = (real *)out_data;
 404         p2[0] = p1[0];
 405         if(dir == GMX_FFT_REAL_TO_COMPLEX)
 406             p2[1] = 0.0;
 407     }
 408
 409     if(dir == GMX_FFT_REAL_TO_COMPLEX)
 410     {
 411         /* FFTPACK only does in-place transforms, so emulate out-of-place
 412          * by copying data to the output array first. This works fine, since
 413          * the complex array must be larger than the real.
 414          */
 415         if( in_data != out_data )
 416         {
 417             p1 = (real *)in_data;
 418             p2 = (real *)out_data;
 419
 420             for(i=0;i<2*(n/2+1);i++)
 421             {
 422                 p2[i] = p1[i];
 423             }
 424         }
 425
 426         /* Elements 0 ..   n-1 in work are used for ffac values,
 427          * Elements n .. 2*n-1 are internal FFTPACK work space.
 428          */
 429         fftpack_rfftf1(n,(real *)out_data,fft->work+n,fft->work,fft->ifac);
 430
 431         /*
 432          * FFTPACK has a slightly more compact storage than we, time to
 433          * convert it: ove most of the array one step up to make room for
 434          * zero imaginary parts.
 435          */
 436         p2 = (real *)out_data;
 437         for(i=n-1;i>0;i--)
 438         {
 439             p2[i+1] = p2[i];
 440         }
 441         /* imaginary zero freq. */
 442         p2[1] = 0;
 443
 444         /* Is n even? */
 445         if( (n & 0x1) == 0 )
 446         {
 447             p2[n+1] = 0;
 448         }
 449
 450     }
 451     else if(dir == GMX_FFT_COMPLEX_TO_REAL)
 452     {
 453         /* FFTPACK only does in-place transforms, and we cannot just copy
 454          * input to output first here since our real array is smaller than
 455          * the complex one. However, since the FFTPACK complex storage format
 456          * is more compact than ours (2 reals) it will fit, so compact it
 457          * and copy on-the-fly to the output array.
 458          */
 459         p1 = (real *) in_data;
 460         p2 = (real *)out_data;
 461
 462         p2[0] = p1[0];
 463         for(i=1;i<n;i++)
 464         {
 465             p2[i] = p1[i+1];
 466         }
 467         fftpack_rfftb1(n,(real *)out_data,fft->work+n,fft->work,fft->ifac);
 468     }
 469     else
 470     {
 471         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
 472         return EINVAL;
 473     }
 474
 475     return 0;
 476 }
 477
 478
 479 int
 480 gmx_fft_2d               (gmx_fft_t                  fft,
 481                           enum gmx_fft_direction     dir,
 482                           void *                     in_data,
 483                           void *                     out_data)
 484 {
 485     int                i,nx,ny;
 486     t_complex *    data;
 487
 488     nx = fft->n;
 489     ny = fft->next->n;
 490
 491     /* FFTPACK only does in-place transforms, so emulate out-of-place
 492      * by copying data to the output array first.
 493      * For 2D there is likely enough data to benefit from memcpy().
 494      */
 495     if( in_data != out_data )
 496     {
 497         memcpy(out_data,in_data,sizeof(t_complex)*nx*ny);
 498     }
 499
 500     /* Much easier to do pointer arithmetic when base has the correct type */
 501     data = (t_complex *)out_data;
 502
 503     /* y transforms */
 504     for(i=0;i<nx;i++)
 505     {
 506         gmx_fft_1d(fft->next,dir,data+i*ny,data+i*ny);
 507     }
 508
 509     /* Transpose in-place to get data in place for x transform now */
 510     gmx_fft_transpose_2d(data,data,nx,ny);
 511
 512     /* x transforms */
 513     for(i=0;i<ny;i++)
 514     {
 515         gmx_fft_1d(fft,dir,data+i*nx,data+i*nx);
 516     }
 517
 518     /* Transpose in-place to get data back in original order */
 519     gmx_fft_transpose_2d(data,data,ny,nx);
 520
 521     return 0;
 522 }
 523
 524
 525
 526 int
 527 gmx_fft_2d_real          (gmx_fft_t                  fft,
 528                           enum gmx_fft_direction     dir,
 529                           void *                     in_data,
 530                           void *                     out_data)
 531 {
 532     int                i,j,nx,ny,nyc;
 533     t_complex *    data;
 534     real *       work;
 535     real *       p1;
 536     real *       p2;
 537
 538     nx=fft->n;
 539     ny=fft->next->n;
 540     /* Number of complex elements in y direction */
 541     nyc=(ny/2+1);
 542
 543     work = fft->work+4*nx;
 544
 545     if(dir==GMX_FFT_REAL_TO_COMPLEX)
 546     {
 547         /* If we are doing an in-place transform the 2D array is already
 548          * properly padded by the user, and we are all set.
 549          *
 550          * For out-of-place there is no array padding, but FFTPACK only
 551          * does in-place FFTs internally, so we need to start by copying
 552          * data from the input to the padded (larger) output array.
 553          */
 554         if( in_data != out_data )
 555         {
 556             p1 = (real *)in_data;
 557             p2 = (real *)out_data;
 558
 559             for(i=0;i<nx;i++)
 560             {
 561                 for(j=0;j<ny;j++)
 562                 {
 563                     p2[i*nyc*2+j] = p1[i*ny+j];
 564                 }
 565             }
 566         }
 567         data = (t_complex *)out_data;
 568
 569         /* y real-to-complex FFTs */
 570         for(i=0;i<nx;i++)
 571         {
 572             gmx_fft_1d_real(fft->next,GMX_FFT_REAL_TO_COMPLEX,data+i*nyc,data+i*nyc);
 573         }
 574
 575         /* Transform to get X data in place */
 576         gmx_fft_transpose_2d(data,data,nx,nyc);
 577
 578         /* Complex-to-complex X FFTs */
 579         for(i=0;i<nyc;i++)
 580         {
 581             gmx_fft_1d(fft,GMX_FFT_FORWARD,data+i*nx,data+i*nx);
 582         }
 583
 584         /* Transpose back */
 585         gmx_fft_transpose_2d(data,data,nyc,nx);
 586
 587     }
 588     else if(dir==GMX_FFT_COMPLEX_TO_REAL)
 589     {
 590         /* An in-place complex-to-real transform is straightforward,
 591          * since the output array must be large enough for the padding to fit.
 592          *
 593          * For out-of-place complex-to-real transforms we cannot just copy
 594          * data to the output array, since it is smaller than the input.
 595          * In this case there's nothing to do but employing temporary work data,
 596          * starting at work+4*nx and using nx*nyc*2 elements.
 597          */
 598         if(in_data != out_data)
 599         {
 600             memcpy(work,in_data,sizeof(t_complex)*nx*nyc);
 601             data = (t_complex *)work;
 602         }
 603         else
 604         {
 605             /* in-place */
 606             data = (t_complex *)out_data;
 607         }
 608
 609         /* Transpose to get X arrays */
 610         gmx_fft_transpose_2d(data,data,nx,nyc);
 611
 612         /* Do X iFFTs */
 613         for(i=0;i<nyc;i++)
 614         {
 615             gmx_fft_1d(fft,GMX_FFT_BACKWARD,data+i*nx,data+i*nx);
 616         }
 617
 618         /* Transpose to get Y arrays */
 619         gmx_fft_transpose_2d(data,data,nyc,nx);
 620
 621         /* Do Y iFFTs */
 622         for(i=0;i<nx;i++)
 623         {
 624             gmx_fft_1d_real(fft->next,GMX_FFT_COMPLEX_TO_REAL,data+i*nyc,data+i*nyc);
 625         }
 626
 627         if( in_data != out_data )
 628         {
 629             /* Output (pointed to by data) is now in padded format.
 630              * Pack it into out_data if we were doing an out-of-place transform.
 631              */
 632             p1 = (real *)data;
 633             p2 = (real *)out_data;
 634
 635             for(i=0;i<nx;i++)
 636             {
 637                 for(j=0;j<ny;j++)
 638                 {
 639                     p2[i*ny+j] = p1[i*nyc*2+j];
 640                 }
 641             }
 642         }
 643     }
 644     else
 645     {
 646         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
 647         return EINVAL;
 648     }
 649
 650     return 0;
 651 }
 652
 653
 654
 655 int
 656 gmx_fft_3d          (gmx_fft_t                  fft,
 657                      enum gmx_fft_direction     dir,
 658                      void *                     in_data,
 659                      void *                     out_data)
 660 {
 661     int              i,nx,ny,nz,rc;
 662     t_complex *  data;
 663     t_complex *  work;
 664     nx=fft->n;
 665     ny=fft->next->n;
 666     nz=fft->next->next->n;
 667
 668     /* First 4*nx positions are FFTPACK workspace, then ours starts */
 669     work = (t_complex *)(fft->work+4*nx);
 670
 671     /* FFTPACK only does in-place transforms, so emulate out-of-place
 672      * by copying data to the output array first.
 673      * For 3D there is likely enough data to benefit from memcpy().
 674      */
 675     if( in_data != out_data )
 676     {
 677         memcpy(out_data,in_data,sizeof(t_complex)*nx*ny*nz);
 678     }
 679
 680     /* Much easier to do pointer arithmetic when base has the correct type */
 681     data = (t_complex *)out_data;
 682
 683     /* Perform z transforms */
 684     for(i=0;i<nx*ny;i++)
 685         gmx_fft_1d(fft->next->next,dir,data+i*nz,data+i*nz);
 686
 687     /* For each X slice, transpose the y & z dimensions inside the slice */
 688     for(i=0;i<nx;i++)
 689     {
 690         gmx_fft_transpose_2d(data+i*ny*nz,data+i*ny*nz,ny,nz);
 691     }
 692
 693     /* Array is now (nx,nz,ny) - perform y transforms */
 694     for(i=0;i<nx*nz;i++)
 695     {
 696         gmx_fft_1d(fft->next,dir,data+i*ny,data+i*ny);
 697     }
 698
 699     /* Transpose back to (nx,ny,nz) */
 700     for(i=0;i<nx;i++)
 701     {
 702         gmx_fft_transpose_2d(data+i*ny*nz,data+i*ny*nz,nz,ny);
 703     }
 704
 705     /* Transpose entire x & y slices to go from
 706      * (nx,ny,nz) to (ny,nx,nz).
 707      * Use work data elements 4*n .. 4*n+2*nz-1.
 708      */
 709     rc=gmx_fft_transpose_2d_nelem(data,data,nx,ny,nz,work);
 710     if( rc != 0)
 711     {
 712         gmx_fatal(FARGS,"Cannot transpose X & Y/Z in gmx_fft_3d().");
 713         return rc;
 714     }
 715
 716     /* Then go from (ny,nx,nz) to (ny,nz,nx) */
 717     for(i=0;i<ny;i++)
 718     {
 719         gmx_fft_transpose_2d(data+i*nx*nz,data+i*nx*nz,nx,nz);
 720     }
 721
 722     /* Perform x transforms */
 723     for(i=0;i<ny*nz;i++)
 724     {
 725         gmx_fft_1d(fft,dir,data+i*nx,data+i*nx);
 726     }
 727
 728     /* Transpose back from (ny,nz,nx) to (ny,nx,nz) */
 729     for(i=0;i<ny;i++)
 730     {
 731         gmx_fft_transpose_2d(data+i*nz*nx,data+i*nz*nx,nz,nx);
 732     }
 733
 734     /* Transpose from (ny,nx,nz) to (nx,ny,nz)
 735      * Use work data elements 4*n .. 4*n+2*nz-1.
 736      */
 737     rc = gmx_fft_transpose_2d_nelem(data,data,ny,nx,nz,work);
 738     if( rc != 0)
 739     {
 740         gmx_fatal(FARGS,"Cannot transpose Y/Z & X in gmx_fft_3d().");
 741         return rc;
 742     }
 743
 744     return 0;
 745 }
 746
 747
 748 int
 749 gmx_fft_3d_real          (gmx_fft_t                  fft,
 750                           enum gmx_fft_direction     dir,
 751                           void *                     in_data,
 752                           void *                     out_data)
 753 {
 754     int              i,j,k;
 755     int              nx,ny,nz,nzc,rc;
 756     t_complex *  data;
 757     t_complex *  work_transp;
 758     t_complex *  work_c2r;
 759     real *     p1;
 760     real *     p2;
 761
 762     nx=fft->n;
 763     ny=fft->next->n;
 764     nz=fft->next->next->n;
 765     nzc=(nz/2+1);
 766
 767
 768     /* First 4*nx positions are FFTPACK workspace, then ours starts.
 769      * We have 2*nx*ny*nzc elements for temp complex-to-real storage when
 770      * doing out-of-place transforms, and another 2*nzc for transpose data.
 771      */
 772     work_c2r    = (t_complex *)(fft->work+4*nx);
 773     work_transp = (t_complex *)(fft->work+4*nx+2*nx*ny*nzc);
 774
 775     /* Much easier to do pointer arithmetic when base has the correct type */
 776     data = (t_complex *)out_data;
 777
 778     if(dir==GMX_FFT_REAL_TO_COMPLEX)
 779     {
 780         /* FFTPACK only does in-place transforms, so emulate out-of-place
 781          * by copying data to the output array first. This is guaranteed to
 782          * work for real-to-complex since complex data is larger than the real.
 783          * For 3D there is likely enough data to benefit from memcpy().
 784          */
 785         if( in_data != out_data )
 786         {
 787             p1 = (real *)in_data;
 788             p2 = (real *)out_data;
 789
 790             for(i=0;i<nx;i++)
 791             {
 792                 for(j=0;j<ny;j++)
 793                 {
 794                     for(k=0;k<nz;k++)
 795                     {
 796                         p2[(i*ny+j)*2*nzc+k] = p1[(i*ny+j)*nz+k];
 797                     }
 798                 }
 799             }
 800         }
 801         data = (t_complex *)out_data;
 802
 803         /* Transform the Y/Z slices real-to-complex */
 804         for(i=0;i<nx;i++)
 805         {
 806             gmx_fft_2d_real(fft->next,dir,data+i*ny*nzc,data+i*ny*nzc);
 807         }
 808
 809         /* Transpose x & y slices to go from
 810          * (nx,ny,nzc) to (ny,nx,nzc).
 811          */
 812         rc=gmx_fft_transpose_2d_nelem(data,data,nx,ny,nzc,work_transp);
 813         if( rc != 0)
 814         {
 815             gmx_fatal(FARGS,"Cannot transpose X & Y/Z gmx_fft_3d_real().");
 816             return rc;
 817         }
 818
 819         /* Then transpose from (ny,nx,nzc) to (ny,nzc,nx) */
 820         for(i=0;i<ny;i++)
 821         {
 822             gmx_fft_transpose_2d(data+i*nx*nzc,data+i*nx*nzc,nx,nzc);
 823         }
 824
 825         /* Perform x transforms */
 826         for(i=0;i<ny*nzc;i++)
 827         {
 828             gmx_fft_1d(fft,GMX_FFT_FORWARD,data+i*nx,data+i*nx);
 829         }
 830
 831         /* Transpose from (ny,nzc,nx) back to (ny,nx,nzc) */
 832         for(i=0;i<ny;i++)
 833         {
 834             gmx_fft_transpose_2d(data+i*nzc*nx,data+i*nzc*nx,nzc,nx);
 835         }
 836
 837         /* Transpose back from (ny,nx,nzc) to (nx,ny,nz) */
 838         rc=gmx_fft_transpose_2d_nelem(data,data,ny,nx,nzc,work_transp);
 839         if( rc != 0)
 840         {
 841             gmx_fatal(FARGS,"Cannot transpose Y/Z & X in gmx_fft_3d_real().");
 842             return rc;
 843         }
 844
 845     }
 846     else if(dir==GMX_FFT_COMPLEX_TO_REAL)
 847     {
 848         /* An in-place complex-to-real transform is straightforward,
 849          * since the output array must be large enough for the padding to fit.
 850          *
 851          * For out-of-place complex-to-real transforms we cannot just copy
 852          * data to the output array, since it is smaller than the input.
 853          * In this case there's nothing to do but employing temporary work data.
 854          */
 855         if(in_data != out_data)
 856         {
 857             memcpy(work_c2r,in_data,sizeof(t_complex)*nx*ny*nzc);
 858             data = (t_complex *)work_c2r;
 859         }
 860         else
 861         {
 862             /* in-place */
 863             data = (t_complex *)out_data;
 864         }
 865
 866         /* Transpose x & y slices to go from
 867         * (nx,ny,nz) to (ny,nx,nz).
 868         */
 869         gmx_fft_transpose_2d_nelem(data,data,nx,ny,nzc,work_transp);
 870
 871         /* Then go from (ny,nx,nzc) to (ny,nzc,nx) */
 872         for(i=0;i<ny;i++)
 873         {
 874             gmx_fft_transpose_2d(data+i*nx*nzc,data+i*nx*nzc,nx,nzc);
 875         }
 876
 877
 878         /* Perform x transforms */
 879         for(i=0;i<ny*nzc;i++)
 880         {
 881             gmx_fft_1d(fft,GMX_FFT_BACKWARD,data+i*nx,data+i*nx);
 882         }
 883
 884         /* Transpose back from (ny,nzc,nx) to (ny,nx,nzc) */
 885         for(i=0;i<ny;i++)
 886         {
 887             gmx_fft_transpose_2d(data+i*nzc*nx,data+i*nzc*nx,nzc,nx);
 888         }
 889
 890         /* Transpose back from (ny,nx,nzc) to (nx,ny,nz) */
 891         gmx_fft_transpose_2d_nelem(data,data,ny,nx,nzc,work_transp);
 892
 893
 894         /* Do 2D complex-to-real */
 895         for(i=0;i<nx;i++)
 896         {
 897             gmx_fft_2d_real(fft->next,dir,data+i*ny*nzc,data+i*ny*nzc);
 898         }
 899
 900         if( in_data != out_data )
 901         {
 902             /* Output (pointed to by data) is now in padded format.
 903              * Pack it into out_data if we were doing an out-of-place transform.
 904              */
 905             p1 = (real *)data;
 906             p2 = (real *)out_data;
 907
 908             for(i=0;i<nx;i++)
 909             {
 910                 for(j=0;j<ny;j++)
 911                 {
 912                     for(k=0;k<nz;k++)
 913                     {
 914                         p2[(i*ny+j)*nz+k] = p1[(i*ny+j)*nzc*2+k];
 915                     }
 916                 }
 917             }
 918         }
 919
 920     }
 921     else
 922     {
 923         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
 924         return EINVAL;
 925     }
 926
 927     return 0;
 928 }
 929
 930
 931
 932
 933 void
 934 gmx_fft_destroy(gmx_fft_t      fft)
 935 {
 936     if(fft != NULL)
 937     {
 938         free(fft->work);
 939         if(fft->next != NULL)
 940             gmx_fft_destroy(fft->next);
 941         free(fft);
 942     }
 943 }
 944 #endif /* GMX_FFT_FFTPACK */