1 /*
2  * Copyright (C) Quantum ESPRESSO group
3  *
4  * This file is distributed under the terms of the
5  * GNU General Public License. See the file `License'
6  * in the root directory of the present distribution,
7  * or http://www.gnu.org/copyleft/gpl.txt .
8  */
9 
10 /* Always compile FFTW beside external driver to be used as fallback
11  * hook for higher level FFT drivers in the FFTX library */
12 
13 #include "fftw_dp.h"
14 #include "fftw_sp.h"
15 
create_plan_1d(fftw_plan * p,int * n,int * idir)16 int create_plan_1d (fftw_plan *p, int *n, int *idir)
17 {
18    fftw_direction dir = ( (*idir < 0) ? FFTW_FORWARD : FFTW_BACKWARD );
19    *p = qe_fftw_create_plan(*n, dir, FFTW_ESTIMATE | FFTW_IN_PLACE);
20    if( *p == NULL ) fprintf(stderr," *** CREATE_PLAN: warning empty plan ***\n");
21 /*   printf(" pointer size = %d, value = %d\n", sizeof ( *p ), *p ); */
22    return 0;
23 }
24 
float_create_plan_1d(float_fftw_plan * p,int * n,int * idir)25 int float_create_plan_1d(float_fftw_plan* p, int* n, int* idir)
26 {
27 	fftw_direction dir = ((*idir < 0) ? FFTW_FORWARD : FFTW_BACKWARD);
28 	*p = qe_float_fftw_create_plan(*n, dir, FFTW_ESTIMATE | FFTW_IN_PLACE);
29 	if (*p == NULL) fprintf(stderr, " *** CREATE_PLAN: warning empty plan ***\n");
30 	/*   printf(" pointer size = %d, value = %d\n", sizeof ( *p ), *p ); */
31 	return 0;
32 }
33 
destroy_plan_1d(fftw_plan * p)34 int destroy_plan_1d (fftw_plan *p)
35 {
36    if ( *p != NULL ) qe_fftw_destroy_plan(*p);
37    else fprintf(stderr," *** DESTROY_PLAN: warning empty plan ***\n");
38    return 0;
39 }
40 
float_destroy_plan_1d(float_fftw_plan * p)41 int float_destroy_plan_1d(float_fftw_plan* p)
42 {
43 	if (*p != NULL) qe_float_fftw_destroy_plan(*p);
44 	else fprintf(stderr, " *** DESTROY_PLAN: warning empty plan ***\n");
45 	return 0;
46 }
47 
create_plan_2d(fftwnd_plan * p,int * n,int * m,int * idir)48 int create_plan_2d (fftwnd_plan *p, int *n, int *m, int *idir)
49 {
50    fftw_direction dir = ( (*idir < 0) ? FFTW_FORWARD : FFTW_BACKWARD );
51    *p = qe_fftw2d_create_plan(*m, *n, dir, FFTW_ESTIMATE | FFTW_IN_PLACE);
52    if( *p == NULL ) fprintf(stderr," *** CREATE_PLAN_2D: warning empty plan ***\n");
53 /*   printf(" pointer size = %d, value = %d\n", sizeof ( *p ), *p ); */
54    return 0;
55 }
56 
destroy_plan_2d(fftwnd_plan * p)57 int destroy_plan_2d (fftwnd_plan *p)
58 {
59    if ( *p != NULL ) qe_fftwnd_destroy_plan(*p);
60    else fprintf(stderr," *** DESTROY_PLAN_2D: warning empty plan ***\n");
61    return 0;
62 }
63 
create_plan_3d(fftwnd_plan * p,int * n,int * m,int * l,int * idir)64 int create_plan_3d (fftwnd_plan *p, int *n, int *m, int *l, int *idir)
65 {
66    fftw_direction dir = ( (*idir < 0) ? FFTW_FORWARD : FFTW_BACKWARD );
67    *p = qe_fftw3d_create_plan(*l, *m, *n, dir, FFTW_ESTIMATE | FFTW_IN_PLACE);
68    if( *p == NULL ) {
69 	fprintf(stderr," *** CREATE_PLAN_3D: warning empty plan ***\n");
70 	fprintf(stderr," *** input was (n,m,l,dir): %d %d %d %d ***\n", *l, *m, *n, *idir);
71    }
72 /*   printf(" pointer size = %d, value = %d\n", sizeof ( *p ), *p ); */
73    return 0;
74 }
75 
destroy_plan_3d(fftwnd_plan * p)76 int destroy_plan_3d (fftwnd_plan *p)
77 
78 {
79    if ( *p != NULL ) qe_fftwnd_destroy_plan(*p);
80    else fprintf(stderr," *** DESTROY_PLAN_3D: warning empty plan ***\n");
81    return 0;
82 }
83 
84 
fft_x_stick(fftw_plan * p,FFTW_COMPLEX * a,int * nx,int * ny,int * nz,int * ldx,int * ldy)85 int fft_x_stick
86 (fftw_plan *p, FFTW_COMPLEX *a, int *nx, int *ny, int *nz, int *ldx, int *ldy )
87 {
88 
89    int i;
90    int xstride, bigstride;
91    int xhowmany, xidist;
92 
93 /* trasform  along x and y */
94    bigstride = (*ldx) * (*ldy);
95 
96    xhowmany = (*ny);
97    xstride  = 1;
98    xidist   = (*ldx);
99 
100    /* ptr = (double *)a; */
101 
102    for(i = 0; i < *nz ; i++) {
103      /* trasform  along x */
104      fftw(*p,xhowmany,&a[i*bigstride],xstride,xidist,0,0,0);
105    }
106    return 0;
107 }
108 
fft_y_stick(fftw_plan * p,FFTW_COMPLEX * a,int * ny,int * ldx)109 int fft_y_stick
110    (fftw_plan *p, FFTW_COMPLEX *a, int *ny, int *ldx )
111 {
112    fftw(*p, 1, a, (*ldx), 1, 0, 0, 0);
113    return 0;
114 }
115 
116 
fft_z_stick(fftw_plan * p,FFTW_COMPLEX * zstick,int * ldz,int * nstick_l)117 int fft_z_stick
118    (fftw_plan *p, FFTW_COMPLEX *zstick, int *ldz, int *nstick_l)
119 {
120    int howmany, idist;
121    howmany = (*nstick_l) ;
122    idist   = (*ldz);
123    fftw(*p, howmany, zstick, 1, idist, 0, 0, 0);
124    return 0;
125 }
126 
fftw_inplace_drv_1d(fftw_plan * p,int * nfft,FFTW_COMPLEX * a,int * inca,int * idist)127 int fftw_inplace_drv_1d
128    (fftw_plan *p, int *nfft, FFTW_COMPLEX *a, int *inca, int *idist)
129 {
130    fftw(*p, (*nfft), a, (*inca), (*idist), 0, 0, 0);
131    return 0;
132 }
133 
fftw_inplace_drv_2d(fftwnd_plan * p,int * nfft,FFTW_COMPLEX * a,int * inca,int * idist)134 int fftw_inplace_drv_2d
135   ( fftwnd_plan *p, int *nfft, FFTW_COMPLEX *a, int *inca, int *idist)
136 {
137    fftwnd( *p, (*nfft), a, (*inca), (*idist), 0, 0, 0 );
138    return 0;
139 }
140 
fftw_inplace_drv_3d(fftwnd_plan * p,int * nfft,FFTW_COMPLEX * a,int * inca,int * idist)141 int fftw_inplace_drv_3d
142    ( fftwnd_plan *p, int *nfft, FFTW_COMPLEX *a, int *inca, int *idist)
143 {
144    fftwnd( *p, (*nfft), a, (*inca), (*idist), 0, 0, 0 );
145    return 0;
146 }
147 
fft_x_stick_single(fftw_plan * p,FFTW_COMPLEX * a,int * nx,int * ny,int * nz,int * ldx,int * ldy)148 int fft_x_stick_single
149 (fftw_plan *p, FFTW_COMPLEX *a, int *nx, int *ny, int *nz, int *ldx, int *ldy )
150 {
151 
152    int xstride, bigstride;
153    int xhowmany, xidist;
154 
155 /* trasform  along x and y */
156    bigstride = (*ldx) * (*ldy);
157 
158    xhowmany = (*ny);
159    xstride  = 1;
160    xidist   = (*ldx);
161 
162    fftw(*p,xhowmany,a,xstride,xidist,0,0,0);
163 
164    return 0;
165 }
166 
167 
fft_z_stick_single(fftw_plan * p,FFTW_COMPLEX * a,int * ldz)168 int fft_z_stick_single (fftw_plan *p, FFTW_COMPLEX *a, int *ldz)
169 {
170   fftw(*p, 1,a, 1, 0, 0, 0, 0);
171 
172   return 0;
173 }
174 
float_fft_z_stick_single(float_fftw_plan * p,FFTW_FLOAT_COMPLEX * a,int * ldz)175 int float_fft_z_stick_single(float_fftw_plan* p, FFTW_FLOAT_COMPLEX* a, int* ldz)
176 {
177 	float_fftw(*p, 1, a, 1, 0, 0, 0, 0);
178 
179 	return 0;
180 }
181 
182 /* Computing the N-Dimensional FFT
183 void fftwnd(fftwnd_plan plan, int howmany,
184             FFTW_COMPLEX *in, int istride, int idist,
185             FFTW_COMPLEX *out, int ostride, int odist);
186 */
187 
188 
189 /*
190 void fftw(fftw_plan plan, int howmany,
191           fftw_complex *in, int istride, int idist,
192           fftw_complex *out, int ostride, int odist);
193 
194 The function fftw computes the one-dimensional Fourier transform,
195 using a plan created by fftw_create_plan (See Section Plan Creation for
196 One-dimensional Transforms.) The function fftw_one provides a simplified
197 interface for the common case of single input array of stride 1.
198 
199 Arguments
200 plan    is the plan created by fftw_create_plan
201 howmany is the number of transforms fftw will compute. It is faster to
202    tell FFTW to compute many transforms, instead of simply calling fftw
203    many times.
204 in, istride and idist describe the input array(s).
205   There are howmany input arrays; the first one is pointed to by in,
206   the second one is pointed to by in + idist, and so on, up to
207   in + (howmany - 1) * idist. Each input array consists of complex numbers,
208   which are not necessarily contiguous in memory. Specifically, in[0] is
209   the first element of the first array, in[istride] is the second element
210   of the first array, and so on. In general, the i-th element of the j-th
211   input array will be in position in[i * istride + j * idist].
212 out, ostride and odist describe the output array(s).
213   The format is the same as for the input array.
214 
215 In-place transforms: If the plan specifies an in-place transform, ostride
216 and odist are always ignored. If out is NULL, out is ignored, too. Otherwise,
217 out is interpreted as a pointer to an array of n complex numbers, that FFTW
218 will use as temporary space to perform the in-place computation. out is used
219 as scratch space and its contents destroyed. In this case, out must be an
220 ordinary array whose elements are contiguous in memory (no striding).
221 
222 */
223