2 * This file is part of the GROMACS molecular simulation package.
4 * Gromacs 4.0 Copyright (c) 1991-2003
5 * David van der Spoel, Erik Lindahl, University of Groningen.
6 * Copyright (c) 2012,2013, 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.
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.
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.
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.
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.
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.
41 #ifdef GMX_FFT_FFTPACK
50 #include "gmx_fatal.h"
53 /** Contents of the FFTPACK fft datatype.
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.
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 */
74 gmx_fft_init_1d(gmx_fft_t * pfft,
82 gmx_fatal(FARGS,"Invalid FFT opaque type pointer.");
87 if( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
95 /* Need 4*n storage for 1D complex FFT */
96 if( (fft->work = (real *)malloc(sizeof(real)*(4*nx))) == NULL)
103 fftpack_cffti1(nx,fft->work,fft->ifac);
112 gmx_fft_init_1d_real(gmx_fft_t * pfft,
120 gmx_fatal(FARGS,"Invalid FFT opaque type pointer.");
125 if( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
133 /* Need 2*n storage for 1D real FFT */
134 if((fft->work = (real *)malloc(sizeof(real)*(2*nx)))==NULL)
141 fftpack_rffti1(nx,fft->work,fft->ifac);
150 gmx_fft_init_2d(gmx_fft_t * pfft,
160 gmx_fatal(FARGS,"Invalid FFT opaque type pointer.");
165 /* Create the X transform */
166 if( (rc = gmx_fft_init_1d(&fft,nx,flags)) != 0)
171 /* Create Y transform as a link from X */
172 if( (rc=gmx_fft_init_1d(&(fft->next),ny,flags)) != 0)
184 gmx_fft_init_2d_real(gmx_fft_t * pfft,
190 int nyc = (ny/2 + 1);
195 gmx_fatal(FARGS,"Invalid FFT opaque type pointer.");
200 /* Create the X transform */
201 if( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
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.
211 if( (fft->work = (real *)malloc(sizeof(real)*(4*nx+2*nx*nyc))) == NULL)
216 fftpack_cffti1(nx,fft->work,fft->ifac);
218 /* Create real Y transform as a link from X */
219 if( (rc=gmx_fft_init_1d_real(&(fft->next),ny,flags)) != 0)
231 gmx_fft_init_3d(gmx_fft_t * pfft,
242 gmx_fatal(FARGS,"Invalid FFT opaque type pointer.");
247 /* Create the X transform */
249 if( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
256 /* Need 4*nx storage for 1D complex FFT, and another
257 * 2*nz elements for gmx_fft_transpose_2d_nelem() storage.
259 if( (fft->work = (real *)malloc(sizeof(real)*(4*nx+2*nz))) == NULL)
265 fftpack_cffti1(nx,fft->work,fft->ifac);
268 /* Create 2D Y/Z transforms as a link from X */
269 if( (rc=gmx_fft_init_2d(&(fft->next),ny,nz,flags)) != 0)
281 gmx_fft_init_3d_real(gmx_fft_t * pfft,
288 int nzc = (nz/2 + 1);
293 gmx_fatal(FARGS,"Invalid FFT opaque type pointer.");
298 /* Create the X transform */
299 if( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
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.
311 if( (fft->work = (real *)malloc(sizeof(real)*(4*nx+2*nx*ny*nzc+2*nzc))) == NULL)
316 fftpack_cffti1(nx,fft->work,fft->ifac);
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)
331 gmx_fft_1d (gmx_fft_t fft,
332 enum gmx_fft_direction dir,
344 p1 = (real *)in_data;
345 p2 = (real *)out_data;
350 /* FFTPACK only does in-place transforms, so emulate out-of-place
351 * by copying data to the output array first.
353 if( in_data != out_data )
355 p1 = (real *)in_data;
356 p2 = (real *)out_data;
358 /* n complex = 2*n real elements */
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.
369 if(dir == GMX_FFT_FORWARD)
371 fftpack_cfftf1(n,(real *)out_data,fft->work+2*n,fft->work,fft->ifac, -1);
373 else if(dir == GMX_FFT_BACKWARD)
375 fftpack_cfftf1(n,(real *)out_data,fft->work+2*n,fft->work,fft->ifac, 1);
379 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
389 gmx_fft_1d_real (gmx_fft_t fft,
390 enum gmx_fft_direction dir,
402 p1 = (real *)in_data;
403 p2 = (real *)out_data;
405 if(dir == GMX_FFT_REAL_TO_COMPLEX)
409 if(dir == GMX_FFT_REAL_TO_COMPLEX)
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.
415 if( in_data != out_data )
417 p1 = (real *)in_data;
418 p2 = (real *)out_data;
420 for(i=0;i<2*(n/2+1);i++)
426 /* Elements 0 .. n-1 in work are used for ffac values,
427 * Elements n .. 2*n-1 are internal FFTPACK work space.
429 fftpack_rfftf1(n,(real *)out_data,fft->work+n,fft->work,fft->ifac);
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.
436 p2 = (real *)out_data;
441 /* imaginary zero freq. */
451 else if(dir == GMX_FFT_COMPLEX_TO_REAL)
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.
459 p1 = (real *) in_data;
460 p2 = (real *)out_data;
467 fftpack_rfftb1(n,(real *)out_data,fft->work+n,fft->work,fft->ifac);
471 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
480 gmx_fft_2d (gmx_fft_t fft,
481 enum gmx_fft_direction dir,
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().
495 if( in_data != out_data )
497 memcpy(out_data,in_data,sizeof(t_complex)*nx*ny);
500 /* Much easier to do pointer arithmetic when base has the correct type */
501 data = (t_complex *)out_data;
506 gmx_fft_1d(fft->next,dir,data+i*ny,data+i*ny);
509 /* Transpose in-place to get data in place for x transform now */
510 gmx_fft_transpose_2d(data,data,nx,ny);
515 gmx_fft_1d(fft,dir,data+i*nx,data+i*nx);
518 /* Transpose in-place to get data back in original order */
519 gmx_fft_transpose_2d(data,data,ny,nx);
527 gmx_fft_2d_real (gmx_fft_t fft,
528 enum gmx_fft_direction dir,
540 /* Number of complex elements in y direction */
543 work = fft->work+4*nx;
545 if(dir==GMX_FFT_REAL_TO_COMPLEX)
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.
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.
554 if( in_data != out_data )
556 p1 = (real *)in_data;
557 p2 = (real *)out_data;
563 p2[i*nyc*2+j] = p1[i*ny+j];
567 data = (t_complex *)out_data;
569 /* y real-to-complex FFTs */
572 gmx_fft_1d_real(fft->next,GMX_FFT_REAL_TO_COMPLEX,data+i*nyc,data+i*nyc);
575 /* Transform to get X data in place */
576 gmx_fft_transpose_2d(data,data,nx,nyc);
578 /* Complex-to-complex X FFTs */
581 gmx_fft_1d(fft,GMX_FFT_FORWARD,data+i*nx,data+i*nx);
585 gmx_fft_transpose_2d(data,data,nyc,nx);
588 else if(dir==GMX_FFT_COMPLEX_TO_REAL)
590 /* An in-place complex-to-real transform is straightforward,
591 * since the output array must be large enough for the padding to fit.
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.
598 if(in_data != out_data)
600 memcpy(work,in_data,sizeof(t_complex)*nx*nyc);
601 data = (t_complex *)work;
606 data = (t_complex *)out_data;
609 /* Transpose to get X arrays */
610 gmx_fft_transpose_2d(data,data,nx,nyc);
615 gmx_fft_1d(fft,GMX_FFT_BACKWARD,data+i*nx,data+i*nx);
618 /* Transpose to get Y arrays */
619 gmx_fft_transpose_2d(data,data,nyc,nx);
624 gmx_fft_1d_real(fft->next,GMX_FFT_COMPLEX_TO_REAL,data+i*nyc,data+i*nyc);
627 if( in_data != out_data )
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.
633 p2 = (real *)out_data;
639 p2[i*ny+j] = p1[i*nyc*2+j];
646 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
656 gmx_fft_3d (gmx_fft_t fft,
657 enum gmx_fft_direction dir,
666 nz=fft->next->next->n;
668 /* First 4*nx positions are FFTPACK workspace, then ours starts */
669 work = (t_complex *)(fft->work+4*nx);
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().
675 if( in_data != out_data )
677 memcpy(out_data,in_data,sizeof(t_complex)*nx*ny*nz);
680 /* Much easier to do pointer arithmetic when base has the correct type */
681 data = (t_complex *)out_data;
683 /* Perform z transforms */
685 gmx_fft_1d(fft->next->next,dir,data+i*nz,data+i*nz);
687 /* For each X slice, transpose the y & z dimensions inside the slice */
690 gmx_fft_transpose_2d(data+i*ny*nz,data+i*ny*nz,ny,nz);
693 /* Array is now (nx,nz,ny) - perform y transforms */
696 gmx_fft_1d(fft->next,dir,data+i*ny,data+i*ny);
699 /* Transpose back to (nx,ny,nz) */
702 gmx_fft_transpose_2d(data+i*ny*nz,data+i*ny*nz,nz,ny);
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.
709 rc=gmx_fft_transpose_2d_nelem(data,data,nx,ny,nz,work);
712 gmx_fatal(FARGS,"Cannot transpose X & Y/Z in gmx_fft_3d().");
716 /* Then go from (ny,nx,nz) to (ny,nz,nx) */
719 gmx_fft_transpose_2d(data+i*nx*nz,data+i*nx*nz,nx,nz);
722 /* Perform x transforms */
725 gmx_fft_1d(fft,dir,data+i*nx,data+i*nx);
728 /* Transpose back from (ny,nz,nx) to (ny,nx,nz) */
731 gmx_fft_transpose_2d(data+i*nz*nx,data+i*nz*nx,nz,nx);
734 /* Transpose from (ny,nx,nz) to (nx,ny,nz)
735 * Use work data elements 4*n .. 4*n+2*nz-1.
737 rc = gmx_fft_transpose_2d_nelem(data,data,ny,nx,nz,work);
740 gmx_fatal(FARGS,"Cannot transpose Y/Z & X in gmx_fft_3d().");
749 gmx_fft_3d_real (gmx_fft_t fft,
750 enum gmx_fft_direction dir,
757 t_complex * work_transp;
758 t_complex * work_c2r;
764 nz=fft->next->next->n;
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.
772 work_c2r = (t_complex *)(fft->work+4*nx);
773 work_transp = (t_complex *)(fft->work+4*nx+2*nx*ny*nzc);
775 /* Much easier to do pointer arithmetic when base has the correct type */
776 data = (t_complex *)out_data;
778 if(dir==GMX_FFT_REAL_TO_COMPLEX)
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().
785 if( in_data != out_data )
787 p1 = (real *)in_data;
788 p2 = (real *)out_data;
796 p2[(i*ny+j)*2*nzc+k] = p1[(i*ny+j)*nz+k];
801 data = (t_complex *)out_data;
803 /* Transform the Y/Z slices real-to-complex */
806 gmx_fft_2d_real(fft->next,dir,data+i*ny*nzc,data+i*ny*nzc);
809 /* Transpose x & y slices to go from
810 * (nx,ny,nzc) to (ny,nx,nzc).
812 rc=gmx_fft_transpose_2d_nelem(data,data,nx,ny,nzc,work_transp);
815 gmx_fatal(FARGS,"Cannot transpose X & Y/Z gmx_fft_3d_real().");
819 /* Then transpose from (ny,nx,nzc) to (ny,nzc,nx) */
822 gmx_fft_transpose_2d(data+i*nx*nzc,data+i*nx*nzc,nx,nzc);
825 /* Perform x transforms */
826 for(i=0;i<ny*nzc;i++)
828 gmx_fft_1d(fft,GMX_FFT_FORWARD,data+i*nx,data+i*nx);
831 /* Transpose from (ny,nzc,nx) back to (ny,nx,nzc) */
834 gmx_fft_transpose_2d(data+i*nzc*nx,data+i*nzc*nx,nzc,nx);
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);
841 gmx_fatal(FARGS,"Cannot transpose Y/Z & X in gmx_fft_3d_real().");
846 else if(dir==GMX_FFT_COMPLEX_TO_REAL)
848 /* An in-place complex-to-real transform is straightforward,
849 * since the output array must be large enough for the padding to fit.
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.
855 if(in_data != out_data)
857 memcpy(work_c2r,in_data,sizeof(t_complex)*nx*ny*nzc);
858 data = (t_complex *)work_c2r;
863 data = (t_complex *)out_data;
866 /* Transpose x & y slices to go from
867 * (nx,ny,nz) to (ny,nx,nz).
869 gmx_fft_transpose_2d_nelem(data,data,nx,ny,nzc,work_transp);
871 /* Then go from (ny,nx,nzc) to (ny,nzc,nx) */
874 gmx_fft_transpose_2d(data+i*nx*nzc,data+i*nx*nzc,nx,nzc);
878 /* Perform x transforms */
879 for(i=0;i<ny*nzc;i++)
881 gmx_fft_1d(fft,GMX_FFT_BACKWARD,data+i*nx,data+i*nx);
884 /* Transpose back from (ny,nzc,nx) to (ny,nx,nzc) */
887 gmx_fft_transpose_2d(data+i*nzc*nx,data+i*nzc*nx,nzc,nx);
890 /* Transpose back from (ny,nx,nzc) to (nx,ny,nz) */
891 gmx_fft_transpose_2d_nelem(data,data,ny,nx,nzc,work_transp);
894 /* Do 2D complex-to-real */
897 gmx_fft_2d_real(fft->next,dir,data+i*ny*nzc,data+i*ny*nzc);
900 if( in_data != out_data )
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.
906 p2 = (real *)out_data;
914 p2[(i*ny+j)*nz+k] = p1[(i*ny+j)*nzc*2+k];
923 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
934 gmx_fft_destroy(gmx_fft_t fft)
939 if(fft->next != NULL)
940 gmx_fft_destroy(fft->next);
944 #endif /* GMX_FFT_FFTPACK */