1 /*
2 Copyright (C) 2002 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
17 02110-1301, USA.
18
19 */
20
21 #include <config.h>
22
23
24 #include <stdio.h>
25 #include <math.h>
26 #include <stdlib.h>
27 #include <string.h>
28 #include <complex.h>
29
30 #ifdef HAVE_NFFT
31 #include "nfft3util.h"
32 #include "nfft3.h"
33 #else
34 typedef int nfft_plan;
35 #endif
36
37 // NFFT FUNCTIONS
FC_FUNC(oct_nfft_init_1d,OCT_NFFT_INIT_1D)38 void FC_FUNC(oct_nfft_init_1d,OCT_NFFT_INIT_1D)
39 (nfft_plan *plan, int *N1, int *M){
40 #ifdef HAVE_NFFT
41 nfft_init_1d(plan, *N1, *M);
42 #endif
43 }
44
FC_FUNC(oct_nfft_init_2d,OCT_NFFT_INIT_2D)45 void FC_FUNC(oct_nfft_init_2d,OCT_NFFT_INIT_2D)
46 (nfft_plan *plan, int *N1, int *N2, int *M){
47 #ifdef HAVE_NFFT
48 nfft_init_2d(plan, *N1, *N2, *M);
49 #endif
50 }
51
FC_FUNC(oct_nfft_init_3d,OCT_NFFT_INIT_3D)52 void FC_FUNC(oct_nfft_init_3d,OCT_NFFT_INIT_3D)
53 (nfft_plan *plan, int *N1, int *N2, int *N3, int *M){
54 #ifdef HAVE_NFFT
55 nfft_init_3d(plan, *N1, *N2, *N3, *M);
56 #endif
57 }
58
FC_FUNC(oct_nfft_init_guru,OCT_NFFT_INIT_GURU)59 void FC_FUNC(oct_nfft_init_guru,OCT_NFFT_INIT_GURU)
60 (nfft_plan *ths, int *d, int *N, int *M, int *n, int *m, unsigned *nfft_flags, unsigned *fftw_flags){
61 #ifdef HAVE_NFFT
62 nfft_init_guru (ths, *d, N, *M, n, *m, *nfft_flags, *fftw_flags);
63 #endif
64 }
65
FC_FUNC(oct_nfft_check,OCT_NFFT_CHECK)66 void FC_FUNC(oct_nfft_check, OCT_NFFT_CHECK)(nfft_plan *ths){
67 #ifdef HAVE_NFFT
68 nfft_check (ths);
69 #endif
70 }
71
FC_FUNC(oct_nfft_finalize,OCT_NFFT_FINALIZE)72 void FC_FUNC(oct_nfft_finalize,OCT_NFFT_FINALIZE)(nfft_plan *ths){
73 #ifdef HAVE_NFFT
74 nfft_finalize (ths);
75 #endif
76 }
77
FC_FUNC(oct_nfft_trafo,OCT_NFFT_TRAFO)78 void FC_FUNC(oct_nfft_trafo,OCT_NFFT_TRAFO)(nfft_plan *ths){
79 #ifdef HAVE_NFFT
80 nfft_trafo (ths);
81 #endif
82 }
83
FC_FUNC(oct_nfft_adjoint,OCT_NFFT_ADJOINT)84 void FC_FUNC(oct_nfft_adjoint,OCT_NFFT_ADJOINT)(nfft_plan *ths){
85 #ifdef HAVE_NFFT
86 nfft_adjoint (ths);
87 #endif
88 }
89
FC_FUNC(oct_nfft_precompute_one_psi_1d,OCT_NFFT_PRECOMPUTE_ONE_PSI_1D)90 void FC_FUNC(oct_nfft_precompute_one_psi_1d, OCT_NFFT_PRECOMPUTE_ONE_PSI_1D)(nfft_plan *plan, int *m, double *X1){
91
92 #ifdef HAVE_NFFT
93 int ii;
94 int M = *m;
95
96 for (ii=0;ii< M;ii++) plan->x[ii] = X1[ii];
97
98 if(plan->nfft_flags & PRE_ONE_PSI) nfft_precompute_one_psi(plan);
99 #endif
100 }
101
FC_FUNC(oct_nfft_precompute_one_psi_2d,OCT_NFFT_PRECOMPUTE_ONE_PSI_2D)102 void FC_FUNC(oct_nfft_precompute_one_psi_2d, OCT_NFFT_PRECOMPUTE_ONE_PSI_2D)(nfft_plan *plan, int *M, double* X1, double* X2){
103
104 #ifdef HAVE_NFFT
105 int ii;
106 int jj;
107
108 for (ii=0; ii< M[0]; ii++){
109 for (jj=0; jj< M[1]; jj++){
110 plan->x[2*(M[1] * ii + jj) + 0] = X1[ii];
111 plan->x[2*(M[1] * ii + jj) + 1] = X2[jj];
112 }
113 }
114
115
116 if(plan->nfft_flags & PRE_ONE_PSI)nfft_precompute_one_psi(plan);
117 #endif
118 }
119
FC_FUNC(oct_nfft_precompute_one_psi_3d,OCT_NFFT_PRECOMPUTE_ONE_PSI_3D)120 void FC_FUNC(oct_nfft_precompute_one_psi_3d, OCT_NFFT_PRECOMPUTE_ONE_PSI_3D)
121 (nfft_plan *plan, int *M, double* X1, double* X2, double* X3){
122 #ifdef HAVE_NFFT
123 int ii,jj,kk;
124
125 for (ii=0;ii< M[0];ii++){
126 for (jj=0;jj< M[1];jj++){
127 for (kk=0;kk< M[2];kk++){
128 plan->x[3*(M[1]*M[2]*ii + M[2]*jj + kk) + 0] = X1[ii];
129 plan->x[3*(M[1]*M[2]*ii + M[2]*jj + kk) + 1] = X2[jj];
130 plan->x[3*(M[1]*M[2]*ii + M[2]*jj + kk) + 2] = X3[kk];
131 }
132 }
133 }
134
135 if(plan->nfft_flags & PRE_ONE_PSI) nfft_precompute_one_psi(plan);
136 #endif
137 }
138
139 // Type dependent functions
140
141
142 // ********** COMPLEX ************
FC_FUNC(zoct_set_f,ZOCT_SET_F)143 void FC_FUNC(zoct_set_f, ZOCT_SET_F)
144 (nfft_plan *plan, int *M, int *DIM, double complex *VAL, int *IX, int *IY, int *IZ){
145
146 #ifdef HAVE_NFFT
147 int dim = *DIM;
148 int ix = *IX;
149 int iy = *IY;
150 int iz = *IZ;
151 double complex val = *VAL;
152
153 switch (dim){
154 case 1:
155 plan->f[ix-1] = val;
156 break;
157 case 2:
158 plan->f[(ix-1)*M[1] + (iy-1)] = val;
159 break;
160 case 3:
161 plan->f[(ix-1)*M[1]*M[2] + (iy-1)*M[2] + (iz-1)] = val;
162 break;
163 }
164 #endif
165 }
166
FC_FUNC(zoct_get_f,ZOCT_GET_F)167 void FC_FUNC(zoct_get_f, ZOCT_GET_F)
168 (nfft_plan *plan, int *M, int *DIM, double complex *val, int *IX, int *IY, int *IZ){
169
170 #ifdef HAVE_NFFT
171 int dim = *DIM;
172 int ix = *IX;
173 int iy = *IY;
174 int iz = *IZ;
175
176
177 switch (dim){
178 case 1:
179 *val = plan->f[ix-1];
180 break;
181 case 2:
182 *val = plan->f[(ix-1)*M[1] + (iy-1)];
183 break;
184 case 3:
185 *val = plan->f[(ix-1)*M[1]*M[2] + (iy-1)*M[2] + (iz-1)];
186 break;
187 }
188 #endif
189 }
190
191
FC_FUNC(zoct_set_f_hat,ZOCT_SET_F_HAT)192 void FC_FUNC(zoct_set_f_hat, ZOCT_SET_F_HAT)
193 (nfft_plan *plan, int *DIM, double complex *VAL, int *IX, int *IY, int *IZ){
194
195 #ifdef HAVE_NFFT
196 int dim = *DIM;
197 int ix = *IX;
198 int iy = *IY;
199 int iz = *IZ;
200 double complex val = *VAL;
201
202 switch (dim){
203 case 1:
204 plan->f_hat[ix-1] = val;
205 break;
206 case 2:
207 plan->f_hat[(ix-1)*plan->N[1] + (iy-1)] = val;
208 break;
209 case 3:
210 plan->f_hat[(ix-1)*plan->N[1]*plan->N[2] + (iy-1)*plan->N[2] + (iz-1)] = val;
211 break;
212 }
213 #endif
214 }
215
FC_FUNC(zoct_get_f_hat,ZOCT_GET_F_HAT)216 void FC_FUNC(zoct_get_f_hat, ZOCT_GET_F_HAT)
217 (nfft_plan *plan, int *DIM, double complex *val, int *IX, int *IY, int *IZ){
218
219 #ifdef HAVE_NFFT
220 int dim = *DIM;
221 int ix = *IX;
222 int iy = *IY;
223 int iz = *IZ;
224
225
226 switch (dim){
227 case 1:
228 *val = plan->f_hat[ix-1];
229 break;
230 case 2:
231 *val = plan->f_hat[(ix-1)*plan->N[1] + (iy-1)];
232 break;
233 case 3:
234 *val = plan->f_hat[(ix-1)*plan->N[1]*plan->N[2] + (iy-1)*plan->N[2] + (iz-1)];
235 break;
236 }
237 #endif
238 }
239
240 // ********** DOUBLE ************
FC_FUNC(doct_set_f,DOCT_SET_F)241 void FC_FUNC(doct_set_f, DOCT_SET_F)
242 (nfft_plan *plan, int *M, int *DIM, double *VAL, int *IX, int *IY, int *IZ){
243
244 #ifdef HAVE_NFFT
245 int dim = *DIM;
246 int ix = *IX;
247 int iy = *IY;
248 int iz = *IZ;
249 double val = *VAL;
250
251
252 switch (dim){
253 case 1:
254 plan->f[ix-1] = val;
255 break;
256 case 2:
257 plan->f[(ix-1)*M[1] + (iy-1)] = val;
258 break;
259 case 3:
260 plan->f[(ix-1)*M[1]*M[2] + (iy-1)*M[2] + (iz-1)] = val;
261 break;
262 }
263 #endif
264 }
265
FC_FUNC(doct_get_f,DOCT_GET_F)266 void FC_FUNC(doct_get_f, DOCT_GET_F)
267 (nfft_plan *plan, int *M, int *DIM, double *val, int *IX, int *IY, int *IZ){
268
269 #ifdef HAVE_NFFT
270 int dim = *DIM;
271 int ix = *IX;
272 int iy = *IY;
273 int iz = *IZ;
274
275
276 switch (dim){
277 case 1:
278 *val = plan->f[ix-1];
279 break;
280 case 2:
281 *val = plan->f[(ix-1)*M[1] + (iy-1)];
282 break;
283 case 3:
284 *val = plan->f[(ix-1)*M[1]*M[2] + (iy-1)*M[2] + (iz-1)];
285 break;
286 }
287 #endif
288 }
289
290
FC_FUNC(doct_set_f_hat,DOCT_SET_F_HAT)291 void FC_FUNC(doct_set_f_hat, DOCT_SET_F_HAT)
292 (nfft_plan *plan, int *DIM, double *VAL, int *IX, int *IY, int *IZ){
293
294 #ifdef HAVE_NFFT
295 int dim = *DIM;
296 int ix = *IX;
297 int iy = *IY;
298 int iz = *IZ;
299 double val = *VAL;
300
301 switch (dim){
302 case 1:
303 plan->f_hat[ix-1] = val;
304 break;
305 case 2:
306 plan->f_hat[(ix-1)*plan->N[1] + (iy-1)] = val;
307 break;
308 case 3:
309 plan->f_hat[(ix-1)*plan->N[1]*plan->N[2] + (iy-1)*plan->N[2] + (iz-1)] = val;
310 break;
311 }
312 #endif
313 }
314
FC_FUNC(doct_get_f_hat,DOCT_GET_F_HAT)315 void FC_FUNC(doct_get_f_hat, DOCT_GET_F_HAT)
316 (nfft_plan *plan, int *DIM, double *val, int *IX, int *IY, int *IZ){
317
318 #ifdef HAVE_NFFT
319 int dim = *DIM;
320 int ix = *IX;
321 int iy = *IY;
322 int iz = *IZ;
323
324
325 switch (dim){
326 case 1:
327 *val = plan->f_hat[ix-1];
328 break;
329 case 2:
330 *val = plan->f_hat[(ix-1)*plan->N[1] + (iy-1)];
331 break;
332 case 3:
333 *val = plan->f_hat[(ix-1)*plan->N[1]*plan->N[2] + (iy-1)*plan->N[2] + (iz-1)];
334 break;
335 }
336 #endif
337 }
338
339
340