Octopus
nfft_f.c
Go to the documentation of this file.
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#include <complex.h>
24#include <math.h>
25#include <stdio.h>
26#include <stdlib.h>
27#include <string.h>
28
29#ifdef HAVE_NFFT
30#include "nfft3.h"
31// #include "nfft3util.h" // No longer used in NFFT3.3.2 and later https://github.com/NFFT/nfft/commit/9f91543bbc7271d1a8cc9c21cfc87287f921bdfe
32#else
33typedef int nfft_plan;
34#endif
35
36// NFFT FUNCTIONS
37void FC_FUNC(oct_nfft_init_1d, OCT_NFFT_INIT_1D)(nfft_plan *plan, int *N1,
38 int *M) {
39#ifdef HAVE_NFFT
40 nfft_init_1d(plan, *N1, *M);
41#endif
42}
43
44void FC_FUNC(oct_nfft_init_2d, OCT_NFFT_INIT_2D)(nfft_plan *plan, int *N1,
45 int *N2, int *M) {
46#ifdef HAVE_NFFT
47 nfft_init_2d(plan, *N1, *N2, *M);
48#endif
49}
50
51void FC_FUNC(oct_nfft_init_3d, OCT_NFFT_INIT_3D)(nfft_plan *plan, int *N1,
52 int *N2, int *N3, int *M) {
53#ifdef HAVE_NFFT
54 nfft_init_3d(plan, *N1, *N2, *N3, *M);
55#endif
56}
57
58void FC_FUNC(oct_nfft_init_guru,
59 OCT_NFFT_INIT_GURU)(nfft_plan *ths, int *d, int *N, int *M, int *n,
60 int *m, unsigned *nfft_flags,
61 unsigned *fftw_flags) {
62#ifdef HAVE_NFFT
63 nfft_init_guru(ths, *d, N, *M, n, *m, *nfft_flags, *fftw_flags);
64#endif
65}
66
67void FC_FUNC(oct_nfft_check, OCT_NFFT_CHECK)(nfft_plan *ths) {
68#ifdef HAVE_NFFT
69 nfft_check(ths);
70#endif
71}
72
73void FC_FUNC(oct_nfft_finalize, OCT_NFFT_FINALIZE)(nfft_plan *ths) {
74#ifdef HAVE_NFFT
75 nfft_finalize(ths);
76#endif
77}
78
79void FC_FUNC(oct_nfft_trafo, OCT_NFFT_TRAFO)(nfft_plan *ths) {
80#ifdef HAVE_NFFT
81 nfft_trafo(ths);
82#endif
83}
84
85void FC_FUNC(oct_nfft_adjoint, OCT_NFFT_ADJOINT)(nfft_plan *ths) {
86#ifdef HAVE_NFFT
87 nfft_adjoint(ths);
88#endif
89}
90
91void FC_FUNC(oct_nfft_precompute_one_psi_1d,
92 OCT_NFFT_PRECOMPUTE_ONE_PSI_1D)(nfft_plan *plan, int *m,
93 double *X1) {
94
95#ifdef HAVE_NFFT
96 int ii;
97 int M = *m;
98
99 for (ii = 0; ii < M; ii++)
100 plan->x[ii] = X1[ii];
101
102#ifdef HAVE_NFFT_3_3
103 if (plan->flags & PRE_ONE_PSI)
104 nfft_precompute_one_psi(plan);
105#else
106 if (plan->nfft_flags & PRE_ONE_PSI)
107 nfft_precompute_one_psi(plan);
108#endif
109#endif
110}
111
112void FC_FUNC(oct_nfft_precompute_one_psi_2d,
113 OCT_NFFT_PRECOMPUTE_ONE_PSI_2D)(nfft_plan *plan, int *M,
114 double *X1, double *X2) {
115
116#ifdef HAVE_NFFT
117 int ii;
118 int jj;
119
120 for (ii = 0; ii < M[0]; ii++) {
121 for (jj = 0; jj < M[1]; jj++) {
122 plan->x[2 * (M[1] * ii + jj) + 0] = X1[ii];
123 plan->x[2 * (M[1] * ii + jj) + 1] = X2[jj];
124 }
125 }
126
127#ifdef HAVE_NFFT_3_3
128 if (plan->flags & PRE_ONE_PSI)
129 nfft_precompute_one_psi(plan);
130#else
131 if (plan->nfft_flags & PRE_ONE_PSI)
132 nfft_precompute_one_psi(plan);
133#endif
134#endif
135}
136
137void FC_FUNC(oct_nfft_precompute_one_psi_3d,
138 OCT_NFFT_PRECOMPUTE_ONE_PSI_3D)(nfft_plan *plan, int *M,
139 double *X1, double *X2,
140 double *X3) {
141#ifdef HAVE_NFFT
142 int ii, jj, kk;
143
144 for (ii = 0; ii < M[0]; ii++) {
145 for (jj = 0; jj < M[1]; jj++) {
146 for (kk = 0; kk < M[2]; kk++) {
147 plan->x[3 * (M[1] * M[2] * ii + M[2] * jj + kk) + 0] = X1[ii];
148 plan->x[3 * (M[1] * M[2] * ii + M[2] * jj + kk) + 1] = X2[jj];
149 plan->x[3 * (M[1] * M[2] * ii + M[2] * jj + kk) + 2] = X3[kk];
150 }
151 }
152 }
153#ifdef HAVE_NFFT_3_3
154 if (plan->flags & PRE_ONE_PSI)
155 nfft_precompute_one_psi(plan);
156#else
157 if (plan->nfft_flags & PRE_ONE_PSI)
158 nfft_precompute_one_psi(plan);
159#endif
160#endif
161}
162
163// Type dependent functions
164
165// ********** COMPLEX ************
166void FC_FUNC(zoct_set_f, ZOCT_SET_F)(nfft_plan *plan, int *M, int *DIM,
167 double complex *VAL, int *IX, int *IY,
168 int *IZ) {
169
170#ifdef HAVE_NFFT
171 int dim = *DIM;
172 int ix = *IX;
173 int iy = *IY;
174 int iz = *IZ;
175 double complex val = *VAL;
176
177 switch (dim) {
178 case 1:
179 plan->f[ix - 1] = val;
180 break;
181 case 2:
182 plan->f[(ix - 1) * M[1] + (iy - 1)] = val;
183 break;
184 case 3:
185 plan->f[(ix - 1) * M[1] * M[2] + (iy - 1) * M[2] + (iz - 1)] = val;
186 break;
187 }
188#endif
189}
190
191void FC_FUNC(zoct_get_f, ZOCT_GET_F)(nfft_plan *plan, int *M, int *DIM,
192 double complex *val, int *IX, int *IY,
193 int *IZ) {
194
195#ifdef HAVE_NFFT
196 int dim = *DIM;
197 int ix = *IX;
198 int iy = *IY;
199 int iz = *IZ;
200
201 switch (dim) {
202 case 1:
203 *val = plan->f[ix - 1];
204 break;
205 case 2:
206 *val = plan->f[(ix - 1) * M[1] + (iy - 1)];
207 break;
208 case 3:
209 *val = plan->f[(ix - 1) * M[1] * M[2] + (iy - 1) * M[2] + (iz - 1)];
210 break;
211 }
212#endif
213}
214
215void FC_FUNC(zoct_set_f_hat, ZOCT_SET_F_HAT)(nfft_plan *plan, int *DIM,
216 double complex *VAL, int *IX,
217 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 double complex val = *VAL;
225
226 switch (dim) {
227 case 1:
228 plan->f_hat[ix - 1] = val;
229 break;
230 case 2:
231 plan->f_hat[(ix - 1) * plan->N[1] + (iy - 1)] = val;
232 break;
233 case 3:
234 plan->f_hat[(ix - 1) * plan->N[1] * plan->N[2] + (iy - 1) * plan->N[2] +
235 (iz - 1)] = val;
236 break;
237 }
238#endif
239}
240
241void FC_FUNC(zoct_get_f_hat, ZOCT_GET_F_HAT)(nfft_plan *plan, int *DIM,
242 double complex *val, int *IX,
243 int *IY, int *IZ) {
244
245#ifdef HAVE_NFFT
246 int dim = *DIM;
247 int ix = *IX;
248 int iy = *IY;
249 int iz = *IZ;
250
251 switch (dim) {
252 case 1:
253 *val = plan->f_hat[ix - 1];
254 break;
255 case 2:
256 *val = plan->f_hat[(ix - 1) * plan->N[1] + (iy - 1)];
257 break;
258 case 3:
259 *val = plan->f_hat[(ix - 1) * plan->N[1] * plan->N[2] +
260 (iy - 1) * plan->N[2] + (iz - 1)];
261 break;
262 }
263#endif
264}
265
266// ********** DOUBLE ************
267void FC_FUNC(doct_set_f, DOCT_SET_F)(nfft_plan *plan, int *M, int *DIM,
268 double *VAL, int *IX, int *IY, int *IZ) {
269
270#ifdef HAVE_NFFT
271 int dim = *DIM;
272 int ix = *IX;
273 int iy = *IY;
274 int iz = *IZ;
275 double val = *VAL;
276
277 switch (dim) {
278 case 1:
279 plan->f[ix - 1] = val;
280 break;
281 case 2:
282 plan->f[(ix - 1) * M[1] + (iy - 1)] = val;
283 break;
284 case 3:
285 plan->f[(ix - 1) * M[1] * M[2] + (iy - 1) * M[2] + (iz - 1)] = val;
286 break;
287 }
288#endif
289}
290
291void FC_FUNC(doct_get_f, DOCT_GET_F)(nfft_plan *plan, int *M, int *DIM,
292 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
300 switch (dim) {
301 case 1:
302 *val = plan->f[ix - 1];
303 break;
304 case 2:
305 *val = plan->f[(ix - 1) * M[1] + (iy - 1)];
306 break;
307 case 3:
308 *val = plan->f[(ix - 1) * M[1] * M[2] + (iy - 1) * M[2] + (iz - 1)];
309 break;
310 }
311#endif
312}
313
314void FC_FUNC(doct_set_f_hat, DOCT_SET_F_HAT)(nfft_plan *plan, int *DIM,
315 double *VAL, int *IX, int *IY,
316 int *IZ) {
317
318#ifdef HAVE_NFFT
319 int dim = *DIM;
320 int ix = *IX;
321 int iy = *IY;
322 int iz = *IZ;
323 double val = *VAL;
324
325 switch (dim) {
326 case 1:
327 plan->f_hat[ix - 1] = val;
328 break;
329 case 2:
330 plan->f_hat[(ix - 1) * plan->N[1] + (iy - 1)] = val;
331 break;
332 case 3:
333 plan->f_hat[(ix - 1) * plan->N[1] * plan->N[2] + (iy - 1) * plan->N[2] +
334 (iz - 1)] = val;
335 break;
336 }
337#endif
338}
339
340void FC_FUNC(doct_get_f_hat, DOCT_GET_F_HAT)(nfft_plan *plan, int *DIM,
341 double *val, int *IX, int *IY,
342 int *IZ) {
343
344#ifdef HAVE_NFFT
345 int dim = *DIM;
346 int ix = *IX;
347 int iy = *IY;
348 int iz = *IZ;
349
350 switch (dim) {
351 case 1:
352 *val = plan->f_hat[ix - 1];
353 break;
354 case 2:
355 *val = plan->f_hat[(ix - 1) * plan->N[1] + (iy - 1)];
356 break;
357 case 3:
358 *val = plan->f_hat[(ix - 1) * plan->N[1] * plan->N[2] +
359 (iy - 1) * plan->N[2] + (iz - 1)];
360 break;
361 }
362#endif
363}