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