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