NFFT  3.4.1
nsfft.c
1 /*
2  * Copyright (c) 2002, 2017 Jens Keiner, Stefan Kunis, Daniel Potts
3  *
4  * This program is free software; you can redistribute it and/or modify it under
5  * the terms of the GNU General Public License as published by the Free Software
6  * Foundation; either version 2 of the License, or (at your option) any later
7  * version.
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12  * details.
13  *
14  * You should have received a copy of the GNU General Public License along with
15  * this program; if not, write to the Free Software Foundation, Inc., 51
16  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  */
18 #include "config.h"
19 
20 #include <stdio.h>
21 #include <math.h>
22 #include <string.h>
23 #include <stdlib.h>
24 #ifdef HAVE_COMPLEX_H
25 #include <complex.h>
26 #endif
27 #include "nfft3.h"
28 #include "infft.h"
29 
30 #define NSFTT_DISABLE_TEST
31 
32 /* computes a 2d ndft by 1d nfft along the dimension 1 times
33  1d ndft along dimension 0
34 */
35 static void short_nfft_trafo_2d(nfft_plan* ths, nfft_plan* plan_1d)
36 {
37  int j,k0;
38  double omega;
39 
40  for(j=0;j<ths->M_total;j++)
41  {
42  ths->f[j]= 0;
43  plan_1d->x[j] = ths->x[ths->d * j + 1];
44  }
45 
46  for(k0=0;k0<ths->N[0];k0++) /* for shorties */
47  {
48  plan_1d->f_hat = ths->f_hat + k0*ths->N[1];
49 
50  nfft_trafo(plan_1d);
51 
52  for(j=0;j<ths->M_total;j++)
53  {
54  omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0];
55  ths->f[j] += plan_1d->f[j] * cexp( - I*2*KPI*omega);
56  }
57  }
58 }
59 
60 static void short_nfft_adjoint_2d(nfft_plan* ths, nfft_plan* plan_1d)
61 {
62  int j,k0;
63  double omega;
64 
65  for(j=0;j<ths->M_total;j++)
66  plan_1d->x[j] = ths->x[ths->d * j + 1];
67 
68  for(k0=0;k0<ths->N[0];k0++) /* for shorties */
69  {
70  for(j=0;j<ths->M_total;j++)
71  {
72  omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0];
73  plan_1d->f[j] = ths->f[j] * cexp( + _Complex_I*2*KPI*omega);
74  }
75 
76  plan_1d->f_hat = ths->f_hat + k0*ths->N[1];
77 
78  nfft_adjoint(plan_1d);
79  }
80 }
81 
82 /* computes a 3d ndft by 1d nfft along the dimension 2 times
83  2d ndft along dimension 0,1
84 */
85 static void short_nfft_trafo_3d_1(nfft_plan* ths, nfft_plan* plan_1d)
86 {
87  int j,k0,k1;
88  double omega;
89 
90  for(j=0;j<ths->M_total;j++)
91  {
92  ths->f[j] = 0;
93  plan_1d->x[j] = ths->x[ths->d * j + 2];
94  }
95 
96  for(k0=0;k0<ths->N[0];k0++) /* for shorties */
97  for(k1=0;k1<ths->N[1];k1++)
98  {
99  plan_1d->f_hat = ths->f_hat + (k0*ths->N[1]+k1)*ths->N[2];
100 
101  nfft_trafo(plan_1d);
102 
103  for(j=0;j<ths->M_total;j++)
104  {
105  omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0]
106  + ((double)(k1 - ths->N[1]/2)) * ths->x[ths->d * j + 1];
107  ths->f[j] += plan_1d->f[j] * cexp( - I*2*KPI*omega);
108  }
109  }
110 }
111 
112 static void short_nfft_adjoint_3d_1(nfft_plan* ths, nfft_plan* plan_1d)
113 {
114  int j,k0,k1;
115  double omega;
116 
117  for(j=0;j<ths->M_total;j++)
118  plan_1d->x[j] = ths->x[ths->d * j + 2];
119 
120  for(k0=0;k0<ths->N[0];k0++) /* for shorties */
121  for(k1=0;k1<ths->N[1];k1++)
122  {
123  for(j=0;j<ths->M_total;j++)
124  {
125  omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0]
126  + ((double)(k1 - ths->N[1]/2)) * ths->x[ths->d * j + 1];
127  plan_1d->f[j] = ths->f[j] * cexp( + _Complex_I*2*KPI*omega);
128  }
129 
130  plan_1d->f_hat = ths->f_hat + (k0*ths->N[1]+k1)*ths->N[2];
131 
132  nfft_adjoint(plan_1d);
133  }
134 }
135 
136 /* computes a 3d ndft by 2d nfft along the dimension 1,2 times
137  1d ndft along dimension 0
138 */
139 static void short_nfft_trafo_3d_2(nfft_plan* ths, nfft_plan* plan_2d)
140 {
141  int j,k0;
142  double omega;
143 
144  for(j=0;j<ths->M_total;j++)
145  {
146  ths->f[j] = 0;
147  plan_2d->x[2*j+0] = ths->x[ths->d * j + 1];
148  plan_2d->x[2*j+1] = ths->x[ths->d * j + 2];
149  }
150 
151  for(k0=0;k0<ths->N[0];k0++) /* for shorties */
152  {
153  plan_2d->f_hat = ths->f_hat + k0*ths->N[1]*ths->N[2];
154 
155  nfft_trafo(plan_2d);
156 
157  for(j=0;j<ths->M_total;j++)
158  {
159  omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0];
160  ths->f[j] += plan_2d->f[j] * cexp( - I*2*KPI*omega);
161  }
162  }
163 }
164 
165 static void short_nfft_adjoint_3d_2(nfft_plan* ths, nfft_plan* plan_2d)
166 {
167  int j,k0;
168  double omega;
169 
170  for(j=0;j<ths->M_total;j++)
171  {
172  plan_2d->x[2*j+0] = ths->x[ths->d * j + 1];
173  plan_2d->x[2*j+1] = ths->x[ths->d * j + 2];
174  }
175 
176  for(k0=0;k0<ths->N[0];k0++) /* for shorties */
177  {
178  for(j=0;j<ths->M_total;j++)
179  {
180  omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0];
181  plan_2d->f[j] = ths->f[j] * cexp( + _Complex_I*2*KPI*omega);
182  }
183 
184  plan_2d->f_hat = ths->f_hat + k0*ths->N[1]*ths->N[2];
185 
186  nfft_adjoint(plan_2d);
187  }
188 }
189 
190 /*---------------------------------------------------------------------------*/
191 
192 #ifdef GAUSSIAN
193 static int index_sparse_to_full_direct_2d(int J, int k)
194 {
195  int N=X(exp2i)(J+2); /* number of full coeffs */
196  int N_B=X(exp2i)(J); /* number in each sparse block */
197 
198  int j=k/N_B; /* consecutive number of Block */
199  int r=j/4; /* level of block */
200 
201  int i, o, a, b,s,l,m1,m2;
202  int k1,k2;
203 
204  if (k>=(J+4)*X(exp2i)(J+1))
205  {
206  printf("Fehler!\n");
207  return(-1);
208  }
209  else
210  {
211  if (r>(J+1)/2) /* center block */
212  {
213  i=k-4*((J+1)/2+1)*N_B;
214  a=X(exp2i)(J/2+1);
215  m1=i/a;
216  m2=i%a;
217  k1=N/2-a/2+m1;
218  k2=N/2-a/2+m2;
219  }
220  else /* no center block */
221  {
222  i=k-j*N_B; /* index in specific block */
223  o=j%4; /* kind of specific block */
224  a=X(exp2i)(r);
225  b=X(exp2i)(J-r);
226  l=MAX(a,b); /* long dimension of block */
227  s=MIN(a,b); /* short dimension of block */
228  m1=i/l;
229  m2=i%l;
230 
231  switch(o)
232  {
233  case 0:
234  {
235  k1=N/2-a/2 ;
236  k2=N/2+ b ;
237 
238  if (b>=a)
239  {
240  k1+=m1;
241  k2+=m2;
242  }
243  else
244  {
245  k1+=m2;
246  k2+=m1;
247  }
248  break;
249  }
250  case 1:
251  {
252  k1=N/2+ b ;
253  k2=N/2-a/2 ;
254 
255  if (b>a)
256  {
257  k1+=m2;
258  k2+=m1;
259  }
260  else
261  {
262  k1+=m1;
263  k2+=m2;
264  }
265  break;
266  }
267  case 2:
268  {
269  k1=N/2-a/2 ;
270  k2=N/2-2*b ;
271 
272  if (b>=a)
273  {
274  k1+=m1;
275  k2+=m2;
276  }
277  else
278  {
279  k1+=m2;
280  k2+=m1;
281  }
282  break;
283  }
284  case 3:
285  {
286  k1=N/2-2*b ;
287  k2=N/2-a/2 ;
288 
289  if (b>a)
290  {
291  k1+=m2;
292  k2+=m1;
293  }
294  else
295  {
296  k1+=m1;
297  k2+=m2;
298  }
299  break;
300  }
301  default:
302  {
303  k1=-1;
304  k2=-1;
305  }
306  }
307  }
308  //printf("m1=%d, m2=%d\n",m1,m2);
309  return(k1*N+k2);
310  }
311 }
312 #endif
313 
314 static inline int index_sparse_to_full_2d(nsfft_plan *ths, int k)
315 {
316  /* only by lookup table */
317  if( k < ths->N_total)
318  return ths->index_sparse_to_full[k];
319  else
320  return -1;
321 }
322 
323 #ifndef NSFTT_DISABLE_TEST
324 static int index_full_to_sparse_2d(int J, int k)
325 {
326  int N=X(exp2i)(J+2); /* number of full coeffs */
327  int N_B=X(exp2i)(J); /* number in each sparse block */
328 
329  int k1=k/N-N/2; /* coordinates in the full grid */
330  int k2=k%N-N/2; /* k1: row, k2: column */
331 
332  int r,a,b;
333 
334  a=X(exp2i)(J/2+1);
335 
336  if ( (k1>=-(a/2)) && (k1<a/2) && (k2>=(-a/2)) && (k2<a/2) )
337  {
338  return(4*((J+1)/2+1)*N_B+(k1+a/2)*a+(k2+a/2));
339  }
340 
341  for (r=0; r<=(J+1)/2; r++)
342  {
343  b=X(exp2i)(r);
344  a=X(exp2i)(J-r);
345  if ( (k1>=-(b/2)) && (k1<(b+1)/2) && (k2>=a) && (k2<2*a) )
346  {
347  if (a>=b)
348  return((4*r+0)*N_B+(k1+b/2)*a+(k2-a));
349  else
350  return((4*r+0)*N_B+(k2-a)*b+(k1+b/2));
351  }
352  else if ( (k1>=a) && (k1<2*a) && (k2>=-(b/2)) && (k2<(b+1)/2) )
353  {
354  if (a>b)
355  return((4*r+1)*N_B+(k2+b/2)*a+(k1-a));
356  else
357  return((4*r+1)*N_B+(k1-a)*b+(k2+b/2));
358  }
359  else if ( (k1>=-(b/2)) && (k1<(b+1)/2) && (k2>=-2*a) && (k2<-a) )
360  {
361  if (a>=b)
362  return((4*r+2)*N_B+(k1+b/2)*a+(k2+2*a));
363  else
364  return((4*r+2)*N_B+(k2+2*a)*b+(k1+b/2));
365  }
366  else if ( (k1>=-2*a) && (k1<-a) && (k2>=-(b/2)) && (k2<(b+1)/2) )
367  {
368  if (a>b)
369  return((4*r+3)*N_B+(k2+b/2)*a+(k1+2*a));
370  else
371  return((4*r+3)*N_B+(k1+2*a)*b+(k2+b/2));
372  }
373  }
374 
375  return(-1);
376 }
377 #endif
378 
379 #ifdef GAUSSIAN
380 static void init_index_sparse_to_full_2d(nsfft_plan *ths)
381 {
382  int k_S;
383 
384  for (k_S=0; k_S<ths->N_total; k_S++)
385  ths->index_sparse_to_full[k_S]=index_sparse_to_full_direct_2d(ths->J, k_S);
386 }
387 #endif
388 
389 #ifdef GAUSSIAN
390 static inline int index_sparse_to_full_3d(nsfft_plan *ths, int k)
391 {
392  /* only by lookup table */
393  if( k < ths->N_total)
394  return ths->index_sparse_to_full[k];
395  else
396  return -1;
397 }
398 #endif
399 
400 #ifndef NSFTT_DISABLE_TEST
401 static int index_full_to_sparse_3d(int J, int k)
402 {
403  int N=X(exp2i)(J+2); /* length of the full grid */
404  int N_B_r; /* size of a sparse block in level r */
405  int sum_N_B_less_r; /* sum N_B_r */
406 
407  int r,a,b;
408 
409  int k3=(k%N)-N/2; /* coordinates in the full grid */
410  int k2=((k/N)%N)-N/2;
411  int k1=k/(N*N)-N/2;
412 
413  a=X(exp2i)(J/2+1); /* length of center block */
414 
415  if((k1>=-(a/2)) && (k1<a/2) && (k2>=(-a/2)) && (k2<a/2) && (k3>=(-a/2)) &&
416  (k3<a/2))
417  {
418  return(6*X(exp2i)(J)*(X(exp2i)((J+1)/2+1)-1)+((k1+a/2)*a+(k2+a/2))*a+
419  (k3+a/2));
420  }
421 
422  sum_N_B_less_r=0;
423  for (r=0; r<=(J+1)/2; r++)
424  {
425  a=X(exp2i)(J-r);
426  b=X(exp2i)(r);
427 
428  N_B_r=a*b*b;
429 
430  /* right - rear - top - left - front - bottom */
431  if ((k1>=a) && (k1<2*a) && (k2>=-(b/2)) && (k2<(b+1)/2) &&
432  (k3>=-(b/2)) && (k3<(b+1)/2)) /* right */
433  {
434  if(a>b)
435  return sum_N_B_less_r+N_B_r*0 + ((k2+b/2)*b+k3+b/2)*a + (k1-a);
436  else
437  return sum_N_B_less_r+N_B_r*0 + ((k1-a)*b+(k2+b/2))*b + (k3+b/2);
438  }
439  else if ((k2>=a) && (k2<2*a) && (k1>=-(b/2)) && (k1<(b+1)/2) &&
440  (k3>=-(b/2)) && (k3<(b+1)/2)) /* rear */
441  {
442  if(a>b)
443  return sum_N_B_less_r+N_B_r*1 + ((k1+b/2)*b+k3+b/2)*a + (k2-a);
444  else if (a==b)
445  return sum_N_B_less_r+N_B_r*1 + ((k1+b/2)*b+(k2-a))*a + (k3+b/2);
446  else
447  return sum_N_B_less_r+N_B_r*1 + ((k2-a)*b+(k1+b/2))*b + (k3+b/2);
448  }
449  else if ((k3>=a) && (k3<2*a) && (k1>=-(b/2)) && (k1<(b+1)/2) &&
450  (k2>=-(b/2)) && (k2<(b+1)/2)) /* top */
451  {
452  if(a>=b)
453  return sum_N_B_less_r+N_B_r*2 + ((k1+b/2)*b+k2+b/2)*a + (k3-a);
454  else
455  return sum_N_B_less_r+N_B_r*2 + ((k3-a)*b+(k1+b/2))*b + (k2+b/2);
456  }
457 
458  else if ((k1>=-2*a) && (k1<-a) && (k2>=-(b/2)) && (k2<(b+1)/2) &&
459  (k3>=-(b/2)) && (k3<(b+1)/2)) /* left */
460  {
461  if(a>b)
462  return sum_N_B_less_r+N_B_r*3 + ((k2+b/2)*b+k3+b/2)*a + (k1+2*a);
463  else
464  return sum_N_B_less_r+N_B_r*3 + ((k1+2*a)*b+(k2+b/2))*b + (k3+b/2);
465  }
466  else if ((k2>=-2*a) && (k2<-a) && (k1>=-(b/2)) && (k1<(b+1)/2) &&
467  (k3>=-(b/2)) && (k3<(b+1)/2)) /* front */
468  {
469  if(a>b)
470  return sum_N_B_less_r+N_B_r*4 + ((k1+b/2)*b+k3+b/2)*a + (k2+2*a);
471  else if (a==b)
472  return sum_N_B_less_r+N_B_r*4 + ((k1+b/2)*b+(k2+2*a))*a + (k3+b/2);
473  else
474  return sum_N_B_less_r+N_B_r*4 + ((k2+2*a)*b+(k1+b/2))*b + (k3+b/2);
475  }
476  else if ((k3>=-2*a) && (k3<-a) && (k1>=-(b/2)) && (k1<(b+1)/2) &&
477  (k2>=-(b/2)) && (k2<(b+1)/2)) /* bottom */
478  {
479  if(a>=b)
480  return sum_N_B_less_r+N_B_r*5 + ((k1+b/2)*b+k2+b/2)*a + (k3+2*a);
481  else
482  return sum_N_B_less_r+N_B_r*5 + ((k3+2*a)*b+(k1+b/2))*b + (k2+b/2);
483  }
484 
485  sum_N_B_less_r+=6*N_B_r;
486  } /* for(r) */
487 
488  return(-1);
489 }
490 #endif
491 
492 #ifdef GAUSSIAN
493 static void init_index_sparse_to_full_3d(nsfft_plan *ths)
494 {
495  int k1,k2,k3,k_s,r;
496  int a,b;
497  int N=X(exp2i)(ths->J+2); /* length of the full grid */
498  int Nc=ths->center_nfft_plan->N[0]; /* length of the center block */
499 
500  for (k_s=0, r=0; r<=(ths->J+1)/2; r++)
501  {
502  a=X(exp2i)(ths->J-r);
503  b=X(exp2i)(r);
504 
505  /* right - rear - top - left - front - bottom */
506 
507  /* right */
508  if(a>b)
509  for(k2=-b/2;k2<(b+1)/2;k2++)
510  for(k3=-b/2;k3<(b+1)/2;k3++)
511  for(k1=a; k1<2*a; k1++,k_s++)
512  ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
513  else
514  for(k1=a; k1<2*a; k1++)
515  for(k2=-b/2;k2<(b+1)/2;k2++)
516  for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
517  ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
518 
519  /* rear */
520  if(a>b)
521  for(k1=-b/2;k1<(b+1)/2;k1++)
522  for(k3=-b/2;k3<(b+1)/2;k3++)
523  for(k2=a; k2<2*a; k2++,k_s++)
524  ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
525  else if(a==b)
526  for(k1=-b/2;k1<(b+1)/2;k1++)
527  for(k2=a; k2<2*a; k2++)
528  for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
529  ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
530  else
531  for(k2=a; k2<2*a; k2++)
532  for(k1=-b/2;k1<(b+1)/2;k1++)
533  for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
534  ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
535 
536  /* top */
537  if(a>=b)
538  for(k1=-b/2;k1<(b+1)/2;k1++)
539  for(k2=-b/2;k2<(b+1)/2;k2++)
540  for(k3=a; k3<2*a; k3++,k_s++)
541  ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
542  else
543  for(k3=a; k3<2*a; k3++)
544  for(k1=-b/2;k1<(b+1)/2;k1++)
545  for(k2=-b/2;k2<(b+1)/2;k2++,k_s++)
546  ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
547 
548  /* left */
549  if(a>b)
550  for(k2=-b/2;k2<(b+1)/2;k2++)
551  for(k3=-b/2;k3<(b+1)/2;k3++)
552  for(k1=-2*a; k1<-a; k1++,k_s++)
553  ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
554  else
555  for(k1=-2*a; k1<-a; k1++)
556  for(k2=-b/2;k2<(b+1)/2;k2++)
557  for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
558  ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
559 
560  /* front */
561  if(a>b)
562  for(k1=-b/2;k1<(b+1)/2;k1++)
563  for(k3=-b/2;k3<(b+1)/2;k3++)
564  for(k2=-2*a; k2<-a; k2++,k_s++)
565  ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
566  else if(a==b)
567  for(k1=-b/2;k1<(b+1)/2;k1++)
568  for(k2=-2*a; k2<-a; k2++)
569  for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
570  ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
571  else
572  for(k2=-2*a; k2<-a; k2++)
573  for(k1=-b/2;k1<(b+1)/2;k1++)
574  for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
575  ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
576 
577  /* top */
578  if(a>=b)
579  for(k1=-b/2;k1<(b+1)/2;k1++)
580  for(k2=-b/2;k2<(b+1)/2;k2++)
581  for(k3=-2*a; k3<-a; k3++,k_s++)
582  ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
583  else
584  for(k3=-2*a; k3<-a; k3++)
585  for(k1=-b/2;k1<(b+1)/2;k1++)
586  for(k2=-b/2;k2<(b+1)/2;k2++,k_s++)
587  ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
588  }
589 
590  /* center */
591  for(k1=-Nc/2;k1<Nc/2;k1++)
592  for(k2=-Nc/2;k2<Nc/2;k2++)
593  for(k3=-Nc/2; k3<Nc/2; k3++,k_s++)
594  ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
595 }
596 #endif
597 
598 /* copies ths->f_hat to ths_plan->f_hat */
599 void nsfft_cp(nsfft_plan *ths, nfft_plan *ths_full_plan)
600 {
601  int k;
602 
603  /* initialize f_hat with zero values */
604  memset(ths_full_plan->f_hat, 0, ths_full_plan->N_total*sizeof(double _Complex));
605 
606  /* copy values at hyperbolic grid points */
607  for(k=0;k<ths->N_total;k++)
608  ths_full_plan->f_hat[ths->index_sparse_to_full[k]]=ths->f_hat[k];
609 
610  /* copy nodes */
611  memcpy(ths_full_plan->x,ths->act_nfft_plan->x,ths->M_total*ths->d*sizeof(double));
612 }
613 
614 #ifndef NSFTT_DISABLE_TEST
615 /* test copy_sparse_to_full */
616 static void test_copy_sparse_to_full_2d(nsfft_plan *ths, nfft_plan *ths_full_plan)
617 {
618  int r;
619  int k1, k2;
620  int a,b;
621  const int J=ths->J; /* N=2^J */
622  const int N=ths_full_plan->N[0]; /* size of full NFFT */
623  const int N_B=X(exp2i)(J); /* size of small blocks */
624 
625  /* copy sparse plan to full plan */
626  nsfft_cp(ths, ths_full_plan);
627 
628  /* show blockwise f_hat */
629  printf("f_hat blockwise\n");
630  for (r=0; r<=(J+1)/2; r++){
631  a=X(exp2i)(J-r); b=X(exp2i)(r);
632 
633  printf("top\n");
634  for (k1=0; k1<a; k1++){
635  for (k2=0; k2<b; k2++){
636  printf("(%1.1f,%1.1f) ", creal(ths->f_hat[(4*r+1)*N_B+ k1*b+k2]),
637  cimag(ths->f_hat[(4*r+1)*N_B+ k1*b+k2]));
638  }
639  printf("\n");
640  }
641 
642  printf("bottom\n");
643  for (k1=0; k1<a; k1++){
644  for (k2=0; k2<b; k2++){
645  printf("(%1.1f,%1.1f) ", creal(ths->f_hat[(4*r+3)*N_B+ k1*b+k2]),
646  cimag(ths->f_hat[(4*r+3)*N_B+ k1*b+k2]));
647  }
648  printf("\n");
649  }
650 
651  printf("right\n");
652  for (k2=0; k2<b; k2++){
653  for (k1=0; k1<a; k1++){
654  printf("(%1.1f,%1.1f) ", creal(ths->f_hat[(4*r+0)*N_B+ k2*a+k1]),
655  cimag(ths->f_hat[(4*r+0)*N_B+ k2*a+k1]));
656  }
657  printf("\n");
658  }
659 
660  printf("left\n");
661  for (k2=0; k2<b; k2++){
662  for (k1=0; k1<a; k1++){
663  printf("(%1.1f,%1.1f) ", creal(ths->f_hat[(4*r+2)*N_B+ k2*a+k1]),
664  cimag(ths->f_hat[(4*r+2)*N_B+ k2*a+k1]));
665  }
666  printf("\n");
667  }
668  }
669 
670  return;
671  /* show full f_hat */
672  printf("full f_hat\n");
673  for (k1=0;k1<N;k1++){
674  for (k2=0;k2<N;k2++){
675  printf("(%1.1f,%1.1f) ", creal(ths_full_plan->f_hat[k1*N+k2]),
676  cimag(ths_full_plan->f_hat[k1*N+k2]));
677  }
678  printf("\n");
679  }
680 }
681 #endif
682 
683 #ifndef NSFTT_DISABLE_TEST
684 static void test_sparse_to_full_2d(nsfft_plan* ths)
685 {
686  int k_S,k1,k2;
687  int N=X(exp2i)(ths->J+2);
688 
689  printf("N=%d\n\n",N);
690 
691  for(k1=0;k1<N;k1++)
692  for(k2=0;k2<N;k2++)
693  {
694  k_S=index_full_to_sparse_2d(ths->J, k1*N+k2);
695  if(k_S!=-1)
696  printf("(%+d, %+d)\t= %+d \t= %+d = %+d \n",k1-N/2,k2-N/2,
697  k1*N+k2, k_S, ths->index_sparse_to_full[k_S]);
698  }
699 }
700 #endif
701 
702 #ifndef NSFTT_DISABLE_TEST
703 static void test_sparse_to_full_3d(nsfft_plan* ths)
704 {
705  int k_S,k1,k2,k3;
706  int N=X(exp2i)(ths->J+2);
707 
708  printf("N=%d\n\n",N);
709 
710  for(k1=0;k1<N;k1++)
711  for(k2=0;k2<N;k2++)
712  for(k3=0;k3<N;k3++)
713  {
714  k_S=index_full_to_sparse_3d(ths->J, (k1*N+k2)*N+k3);
715  if(k_S!=-1)
716  printf("(%d, %d, %d)\t= %d \t= %d = %d \n",k1-N/2,k2-N/2,k3-N/2,
717  (k1*N+k2)*N+k3,k_S, ths->index_sparse_to_full[k_S]);
718  }
719 }
720 #endif
721 
722 
723 void nsfft_init_random_nodes_coeffs(nsfft_plan *ths)
724 {
725  int j;
726 
727  /* init frequencies */
728  nfft_vrand_unit_complex(ths->f_hat, ths->N_total);
729 
730  /* init nodes */
731  nfft_vrand_shifted_unit_double(ths->act_nfft_plan->x, ths->d * ths->M_total);
732 
733  if(ths->d==2)
734  for(j=0;j<ths->M_total;j++)
735  {
736  ths->x_transposed[2*j+0]=ths->act_nfft_plan->x[2*j+1];
737  ths->x_transposed[2*j+1]=ths->act_nfft_plan->x[2*j+0];
738  }
739  else /* this->d==3 */
740  for(j=0;j<ths->M_total;j++)
741  {
742  ths->x_102[3*j+0]=ths->act_nfft_plan->x[3*j+1];
743  ths->x_102[3*j+1]=ths->act_nfft_plan->x[3*j+0];
744  ths->x_102[3*j+2]=ths->act_nfft_plan->x[3*j+2];
745 
746  ths->x_201[3*j+0]=ths->act_nfft_plan->x[3*j+2];
747  ths->x_201[3*j+1]=ths->act_nfft_plan->x[3*j+0];
748  ths->x_201[3*j+2]=ths->act_nfft_plan->x[3*j+1];
749 
750  ths->x_120[3*j+0]=ths->act_nfft_plan->x[3*j+1];
751  ths->x_120[3*j+1]=ths->act_nfft_plan->x[3*j+2];
752  ths->x_120[3*j+2]=ths->act_nfft_plan->x[3*j+0];
753 
754  ths->x_021[3*j+0]=ths->act_nfft_plan->x[3*j+0];
755  ths->x_021[3*j+1]=ths->act_nfft_plan->x[3*j+2];
756  ths->x_021[3*j+2]=ths->act_nfft_plan->x[3*j+1];
757  }
758 }
759 
760 static void nsdft_trafo_2d(nsfft_plan *ths)
761 {
762  int j,k_S,k_L,k0,k1;
763  double omega;
764  int N=X(exp2i)(ths->J+2);
765 
766  memset(ths->f,0,ths->M_total*sizeof(double _Complex));
767 
768  for(k_S=0;k_S<ths->N_total;k_S++)
769  {
770  k_L=ths->index_sparse_to_full[k_S];
771  k0=k_L / N;
772  k1=k_L % N;
773 
774  for(j=0;j<ths->M_total;j++)
775  {
776  omega =
777  ((double)(k0 - N/2)) * ths->act_nfft_plan->x[2 * j + 0] +
778  ((double)(k1 - N/2)) * ths->act_nfft_plan->x[2 * j + 1];
779  ths->f[j] += ths->f_hat[k_S] * cexp( - I*2*KPI*omega);
780  }
781  }
782 } /* void nsdft_trafo_2d */
783 
784 static void nsdft_trafo_3d(nsfft_plan *ths)
785 {
786  int j,k_S,k0,k1,k2;
787  double omega;
788  int N=X(exp2i)(ths->J+2);
789  int k_L;
790 
791  memset(ths->f,0,ths->M_total*sizeof(double _Complex));
792 
793  for(k_S=0;k_S<ths->N_total;k_S++)
794  {
795  k_L=ths->index_sparse_to_full[k_S];
796 
797  k0=k_L/(N*N);
798  k1=(k_L/N)%N;
799  k2=k_L%N;
800 
801  for(j=0;j<ths->M_total;j++)
802  {
803  omega =
804  ((double)(k0 - N/2)) * ths->act_nfft_plan->x[3 * j + 0] +
805  ((double)(k1 - N/2)) * ths->act_nfft_plan->x[3 * j + 1] +
806  ((double)(k2 - N/2)) * ths->act_nfft_plan->x[3 * j + 2];
807  ths->f[j] += ths->f_hat[k_S] * cexp( - I*2*KPI*omega);
808  }
809  }
810 } /* void nsdft_trafo_3d */
811 
812 void nsfft_trafo_direct(nsfft_plan *ths)
813 {
814  if(ths->d==2)
815  nsdft_trafo_2d(ths);
816  else
817  nsdft_trafo_3d(ths);
818 }
819 
820 static void nsdft_adjoint_2d(nsfft_plan *ths)
821 {
822  int j,k_S,k_L,k0,k1;
823  double omega;
824  int N=X(exp2i)(ths->J+2);
825 
826  memset(ths->f_hat,0,ths->N_total*sizeof(double _Complex));
827 
828  for(k_S=0;k_S<ths->N_total;k_S++)
829  {
830  k_L=ths->index_sparse_to_full[k_S];
831  k0=k_L / N;
832  k1=k_L % N;
833 
834  for(j=0;j<ths->M_total;j++)
835  {
836  omega =
837  ((double)(k0 - N/2)) * ths->act_nfft_plan->x[2 * j + 0] +
838  ((double)(k1 - N/2)) * ths->act_nfft_plan->x[2 * j + 1];
839  ths->f_hat[k_S] += ths->f[j] * cexp( + _Complex_I*2*KPI*omega);
840  }
841  }
842 } /* void nsdft_adjoint_2d */
843 
844 static void nsdft_adjoint_3d(nsfft_plan *ths)
845 {
846  int j,k_S,k0,k1,k2;
847  double omega;
848  int N=X(exp2i)(ths->J+2);
849  int k_L;
850 
851  memset(ths->f_hat,0,ths->N_total*sizeof(double _Complex));
852 
853  for(k_S=0;k_S<ths->N_total;k_S++)
854  {
855  k_L=ths->index_sparse_to_full[k_S];
856 
857  k0=k_L/(N*N);
858  k1=(k_L/N)%N;
859  k2=k_L%N;
860 
861  for(j=0;j<ths->M_total;j++)
862  {
863  omega =
864  ((double)(k0 - N/2)) * ths->act_nfft_plan->x[3 * j + 0] +
865  ((double)(k1 - N/2)) * ths->act_nfft_plan->x[3 * j + 1] +
866  ((double)(k2 - N/2)) * ths->act_nfft_plan->x[3 * j + 2];
867  ths->f_hat[k_S] += ths->f[j] * cexp( + _Complex_I*2*KPI*omega);
868  }
869  }
870 } /* void nsdft_adjoint_3d */
871 
872 void nsfft_adjoint_direct(nsfft_plan *ths)
873 {
874  if(ths->d==2)
875  nsdft_adjoint_2d(ths);
876  else
877  nsdft_adjoint_3d(ths);
878 }
879 
880 static void nsfft_trafo_2d(nsfft_plan *ths)
881 {
882  int r,rr,j;
883  double temp;
884 
885  int M=ths->M_total;
886  int J=ths->J;
887 
888  /* center */
889  ths->center_nfft_plan->f_hat=ths->f_hat+4*((J+1)/2+1)*X(exp2i)(J);
890 
891  if (ths->center_nfft_plan->N[0]<=ths->center_nfft_plan->m)
892  nfft_trafo_direct(ths->center_nfft_plan);
893  else
894  nfft_trafo(ths->center_nfft_plan);
895 
896  for (j=0; j<M; j++)
897  ths->f[j] = ths->center_nfft_plan->f[j];
898 
899  for(rr=0;rr<=(J+1)/2;rr++)
900  {
901  r=MIN(rr,J-rr);
902  ths->act_nfft_plan->my_fftw_plan1 = ths->set_fftw_plan1[r];
903  ths->act_nfft_plan->N[0]=X(exp2i)(r); ths->act_nfft_plan->n[0]=ths->sigma*ths->act_nfft_plan->N[0];
904  ths->act_nfft_plan->N[1]=X(exp2i)(J-r); ths->act_nfft_plan->n[1]=ths->sigma*ths->act_nfft_plan->N[1];
905 
906  /*printf("%d x %d\n",ths->act_nfft_plan->N[0],ths->act_nfft_plan->N[1]);*/
907 
908  temp=-3.0*KPI*X(exp2i)(J-rr);
909 
910  /* right */
911  ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+0)*X(exp2i)(J);
912 
913  if(r<rr)
914  RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
915 
916  if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
917  if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
918  nfft_trafo_direct(ths->act_nfft_plan);
919  else
920  short_nfft_trafo_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
921  else
922  nfft_trafo(ths->act_nfft_plan);
923 
924  if(r<rr)
925  RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
926 
927  for (j=0; j<M; j++)
928  ths->f[j] += ths->act_nfft_plan->f[j] *
929  cexp( + _Complex_I*temp*ths->act_nfft_plan->x[2*j+1]);
930 
931  /* top */
932  ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+1)*X(exp2i)(J);
933 
934  if((r==rr)&&(J-rr!=rr))
935  RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
936 
937  if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
938  if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
939  nfft_trafo_direct(ths->act_nfft_plan);
940  else
941  short_nfft_trafo_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
942  else
943  nfft_trafo(ths->act_nfft_plan);
944 
945  if((r==rr)&&(J-rr!=rr))
946  RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
947 
948  for (j=0; j<M; j++)
949  ths->f[j] += ths->act_nfft_plan->f[j] *
950  cexp( + _Complex_I*temp*ths->act_nfft_plan->x[2*j+0]);
951 
952  /* left */
953  ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+2)*X(exp2i)(J);
954 
955  if(r<rr)
956  RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
957 
958  if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
959  if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
960  nfft_trafo_direct(ths->act_nfft_plan);
961  else
962  short_nfft_trafo_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
963  else
964  nfft_trafo(ths->act_nfft_plan);
965 
966  if(r<rr)
967  RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
968 
969  for (j=0; j<M; j++)
970  ths->f[j] += ths->act_nfft_plan->f[j] *
971  cexp( - I*temp*ths->act_nfft_plan->x[2*j+1]);
972 
973  /* bottom */
974  ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+3)*X(exp2i)(J);
975 
976  if((r==rr)&&(J-rr!=rr))
977  RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
978 
979  if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
980  if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
981  nfft_trafo_direct(ths->act_nfft_plan);
982  else
983  short_nfft_trafo_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
984  else
985  nfft_trafo(ths->act_nfft_plan);
986 
987  if((r==rr)&&(J-rr!=rr))
988  RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
989 
990  for (j=0; j<M; j++)
991  ths->f[j] += ths->act_nfft_plan->f[j] *
992  cexp( - I*temp*ths->act_nfft_plan->x[2*j+0]);
993  } /* for(rr) */
994 } /* void nsfft_trafo_2d */
995 
996 static void nsfft_adjoint_2d(nsfft_plan *ths)
997 {
998  int r,rr,j;
999  double temp;
1000 
1001  int M=ths->M_total;
1002  int J=ths->J;
1003 
1004  /* center */
1005  for (j=0; j<M; j++)
1006  ths->center_nfft_plan->f[j] = ths->f[j];
1007 
1008  ths->center_nfft_plan->f_hat=ths->f_hat+4*((J+1)/2+1)*X(exp2i)(J);
1009 
1010  if (ths->center_nfft_plan->N[0]<=ths->center_nfft_plan->m)
1011  nfft_adjoint_direct(ths->center_nfft_plan);
1012  else
1013  nfft_adjoint(ths->center_nfft_plan);
1014 
1015  for(rr=0;rr<=(J+1)/2;rr++)
1016  {
1017  r=MIN(rr,J-rr);
1018  ths->act_nfft_plan->my_fftw_plan2 = ths->set_fftw_plan2[r];
1019  ths->act_nfft_plan->N[0]=X(exp2i)(r); ths->act_nfft_plan->n[0]=ths->sigma*ths->act_nfft_plan->N[0];
1020  ths->act_nfft_plan->N[1]=X(exp2i)(J-r); ths->act_nfft_plan->n[1]=ths->sigma*ths->act_nfft_plan->N[1];
1021 
1022  /*printf("%d x %d\n",ths->act_nfft_plan->N[0],ths->act_nfft_plan->N[1]);*/
1023 
1024  temp=-3.0*KPI*X(exp2i)(J-rr);
1025 
1026  /* right */
1027  ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+0)*X(exp2i)(J);
1028 
1029  for (j=0; j<M; j++)
1030  ths->act_nfft_plan->f[j]= ths->f[j] *
1031  cexp( - I*temp*ths->act_nfft_plan->x[2*j+1]);
1032 
1033  if(r<rr)
1034  RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
1035 
1036  if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1037  if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1038  nfft_adjoint_direct(ths->act_nfft_plan);
1039  else
1040  short_nfft_adjoint_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1041  else
1042  nfft_adjoint(ths->act_nfft_plan);
1043 
1044  if(r<rr)
1045  RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
1046 
1047  /* top */
1048  ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+1)*X(exp2i)(J);
1049 
1050  for (j=0; j<M; j++)
1051  ths->act_nfft_plan->f[j]= ths->f[j] *
1052  cexp( - I*temp*ths->act_nfft_plan->x[2*j+0]);
1053 
1054  if((r==rr)&&(J-rr!=rr))
1055  RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
1056 
1057  if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1058  if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1059  nfft_adjoint_direct(ths->act_nfft_plan);
1060  else
1061  short_nfft_adjoint_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1062  else
1063  nfft_adjoint(ths->act_nfft_plan);
1064 
1065  if((r==rr)&&(J-rr!=rr))
1066  RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
1067 
1068  /* left */
1069  ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+2)*X(exp2i)(J);
1070 
1071  for (j=0; j<M; j++)
1072  ths->act_nfft_plan->f[j]= ths->f[j] *
1073  cexp( + _Complex_I*temp*ths->act_nfft_plan->x[2*j+1]);
1074 
1075  if(r<rr)
1076  RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
1077 
1078  if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1079  if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1080  nfft_adjoint_direct(ths->act_nfft_plan);
1081  else
1082  short_nfft_adjoint_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1083  else
1084  nfft_adjoint(ths->act_nfft_plan);
1085 
1086  if(r<rr)
1087  RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
1088 
1089  /* bottom */
1090  ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+3)*X(exp2i)(J);
1091 
1092  for (j=0; j<M; j++)
1093  ths->act_nfft_plan->f[j]= ths->f[j] *
1094  cexp( + _Complex_I*temp*ths->act_nfft_plan->x[2*j+0]);
1095 
1096  if((r==rr)&&(J-rr!=rr))
1097  RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
1098 
1099  if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1100  if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1101  nfft_adjoint_direct(ths->act_nfft_plan);
1102  else
1103  short_nfft_adjoint_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1104  else
1105  nfft_adjoint(ths->act_nfft_plan);
1106 
1107  if((r==rr)&&(J-rr!=rr))
1108  RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
1109  } /* for(rr) */
1110 } /* void nsfft_adjoint_2d */
1111 
1112 static void nsfft_trafo_3d(nsfft_plan *ths)
1113 {
1114  int r,rr,j;
1115  double temp;
1116  int sum_N_B_less_r,N_B_r,a,b;
1117 
1118  int M=ths->M_total;
1119  int J=ths->J;
1120 
1121  /* center */
1122  ths->center_nfft_plan->f_hat=ths->f_hat+6*X(exp2i)(J)*(X(exp2i)((J+1)/2+1)-1);
1123 
1124  if (ths->center_nfft_plan->N[0]<=ths->center_nfft_plan->m)
1125  nfft_trafo_direct(ths->center_nfft_plan);
1126  else
1127  nfft_trafo(ths->center_nfft_plan);
1128 
1129  for (j=0; j<M; j++)
1130  ths->f[j] = ths->center_nfft_plan->f[j];
1131 
1132  sum_N_B_less_r=0;
1133  for(rr=0;rr<=(J+1)/2;rr++)
1134  {
1135  a=X(exp2i)(J-rr);
1136  b=X(exp2i)(rr);
1137 
1138  N_B_r=a*b*b;
1139 
1140  r=MIN(rr,J-rr);
1141  ths->act_nfft_plan->my_fftw_plan1 = ths->set_fftw_plan1[rr];
1142 
1143  ths->act_nfft_plan->N[0]=X(exp2i)(r);
1144  if(a<b)
1145  ths->act_nfft_plan->N[1]=X(exp2i)(J-r);
1146  else
1147  ths->act_nfft_plan->N[1]=X(exp2i)(r);
1148  ths->act_nfft_plan->N[2]=X(exp2i)(J-r);
1149 
1150  /*printf("\n\n%d x %d x %d:\t",ths->act_nfft_plan->N[0],ths->act_nfft_plan->N[1],ths->act_nfft_plan->N[2]); fflush(stdout);*/
1151 
1152  ths->act_nfft_plan->N_total=ths->act_nfft_plan->N[0]*ths->act_nfft_plan->N[1]*ths->act_nfft_plan->N[2];
1153  ths->act_nfft_plan->n[0]=ths->sigma*ths->act_nfft_plan->N[0];
1154  ths->act_nfft_plan->n[1]=ths->sigma*ths->act_nfft_plan->N[1];
1155  ths->act_nfft_plan->n[2]=ths->sigma*ths->act_nfft_plan->N[2];
1156  ths->act_nfft_plan->n_total=ths->act_nfft_plan->n[0]*ths->act_nfft_plan->n[1]*ths->act_nfft_plan->n[2];
1157 
1158  /* only for right - rear - top */
1159  if((J==0)||((J==1)&&(rr==1)))
1160  temp=-2.0*KPI;
1161  else
1162  temp=-3.0*KPI*X(exp2i)(J-rr);
1163 
1164  /* right */
1165  ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*0;
1166 
1167  if(a>b)
1168  RSWAP(ths->act_nfft_plan->x,ths->x_120);
1169 
1170  if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1171  if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1172  if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1173  nfft_trafo_direct(ths->act_nfft_plan);
1174  else
1175  short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1176  else
1177  short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1178  else
1179  nfft_trafo(ths->act_nfft_plan);
1180 
1181  if(a>b)
1182  RSWAP(ths->act_nfft_plan->x,ths->x_120);
1183 
1184  for (j=0; j<M; j++)
1185  ths->f[j] += ths->act_nfft_plan->f[j] *
1186  cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+0]);
1187 
1188  /* rear */
1189  ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*1;
1190 
1191  if(a>b)
1192  RSWAP(ths->act_nfft_plan->x,ths->x_021);
1193  if(a<b)
1194  RSWAP(ths->act_nfft_plan->x,ths->x_102);
1195 
1196  if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1197  if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1198  if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1199  nfft_trafo_direct(ths->act_nfft_plan);
1200  else
1201  short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1202  else
1203  short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1204  else
1205  nfft_trafo(ths->act_nfft_plan);
1206 
1207  if(a>b)
1208  RSWAP(ths->act_nfft_plan->x,ths->x_021);
1209  if(a<b)
1210  RSWAP(ths->act_nfft_plan->x,ths->x_102);
1211 
1212  for (j=0; j<M; j++)
1213  ths->f[j] += ths->act_nfft_plan->f[j] *
1214  cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+1]);
1215 
1216  /* top */
1217  ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*2;
1218 
1219  if(a<b)
1220  RSWAP(ths->act_nfft_plan->x,ths->x_201);
1221 
1222  if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1223  if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1224  if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1225  nfft_trafo_direct(ths->act_nfft_plan);
1226  else
1227  short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1228  else
1229  short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1230  else
1231  nfft_trafo(ths->act_nfft_plan);
1232 
1233  if(a<b)
1234  RSWAP(ths->act_nfft_plan->x,ths->x_201);
1235 
1236  for (j=0; j<M; j++)
1237  ths->f[j] += ths->act_nfft_plan->f[j] *
1238  cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+2]);
1239 
1240  /* only for left - front - bottom */
1241  if((J==0)||((J==1)&&(rr==1)))
1242  temp=-4.0*KPI;
1243  else
1244  temp=-3.0*KPI*X(exp2i)(J-rr);
1245 
1246  /* left */
1247  ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*3;
1248 
1249  if(a>b)
1250  RSWAP(ths->act_nfft_plan->x,ths->x_120);
1251 
1252  if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1253  if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1254  if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1255  nfft_trafo_direct(ths->act_nfft_plan);
1256  else
1257  short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1258  else
1259  short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1260  else
1261  nfft_trafo(ths->act_nfft_plan);
1262 
1263  if(a>b)
1264  RSWAP(ths->act_nfft_plan->x,ths->x_120);
1265 
1266  for (j=0; j<M; j++)
1267  ths->f[j] += ths->act_nfft_plan->f[j] *
1268  cexp( - I*temp*ths->act_nfft_plan->x[3*j+0]);
1269 
1270  /* front */
1271  ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*4;
1272 
1273  if(a>b)
1274  RSWAP(ths->act_nfft_plan->x,ths->x_021);
1275  if(a<b)
1276  RSWAP(ths->act_nfft_plan->x,ths->x_102);
1277 
1278  if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1279  if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1280  if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1281  nfft_trafo_direct(ths->act_nfft_plan);
1282  else
1283  short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1284  else
1285  short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1286  else
1287  nfft_trafo(ths->act_nfft_plan);
1288 
1289  if(a>b)
1290  RSWAP(ths->act_nfft_plan->x,ths->x_021);
1291  if(a<b)
1292  RSWAP(ths->act_nfft_plan->x,ths->x_102);
1293 
1294  for (j=0; j<M; j++)
1295  ths->f[j] += ths->act_nfft_plan->f[j] *
1296  cexp( - I*temp*ths->act_nfft_plan->x[3*j+1]);
1297 
1298  /* bottom */
1299  ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*5;
1300 
1301  if(a<b)
1302  RSWAP(ths->act_nfft_plan->x,ths->x_201);
1303 
1304  if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1305  if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1306  if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1307  nfft_trafo_direct(ths->act_nfft_plan);
1308  else
1309  short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1310  else
1311  short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1312  else
1313  nfft_trafo(ths->act_nfft_plan);
1314 
1315  if(a<b)
1316  RSWAP(ths->act_nfft_plan->x,ths->x_201);
1317 
1318  for (j=0; j<M; j++)
1319  ths->f[j] += ths->act_nfft_plan->f[j] *
1320  cexp( - I*temp*ths->act_nfft_plan->x[3*j+2]);
1321 
1322  sum_N_B_less_r+=6*N_B_r;
1323  } /* for(rr) */
1324 } /* void nsfft_trafo_3d */
1325 
1326 static void nsfft_adjoint_3d(nsfft_plan *ths)
1327 {
1328  int r,rr,j;
1329  double temp;
1330  int sum_N_B_less_r,N_B_r,a,b;
1331 
1332  int M=ths->M_total;
1333  int J=ths->J;
1334 
1335  /* center */
1336  for (j=0; j<M; j++)
1337  ths->center_nfft_plan->f[j] = ths->f[j];
1338 
1339  ths->center_nfft_plan->f_hat=ths->f_hat+6*X(exp2i)(J)*(X(exp2i)((J+1)/2+1)-1);
1340 
1341  if (ths->center_nfft_plan->N[0]<=ths->center_nfft_plan->m)
1342  nfft_adjoint_direct(ths->center_nfft_plan);
1343  else
1344  nfft_adjoint(ths->center_nfft_plan);
1345 
1346  sum_N_B_less_r=0;
1347  for(rr=0;rr<=(J+1)/2;rr++)
1348  {
1349  a=X(exp2i)(J-rr);
1350  b=X(exp2i)(rr);
1351 
1352  N_B_r=a*b*b;
1353 
1354  r=MIN(rr,J-rr);
1355  ths->act_nfft_plan->my_fftw_plan1 = ths->set_fftw_plan1[rr];
1356  ths->act_nfft_plan->my_fftw_plan2 = ths->set_fftw_plan2[rr];
1357 
1358  ths->act_nfft_plan->N[0]=X(exp2i)(r);
1359  if(a<b)
1360  ths->act_nfft_plan->N[1]=X(exp2i)(J-r);
1361  else
1362  ths->act_nfft_plan->N[1]=X(exp2i)(r);
1363  ths->act_nfft_plan->N[2]=X(exp2i)(J-r);
1364 
1365  /*printf("\n\n%d x %d x %d:\t",ths->act_nfft_plan->N[0],ths->act_nfft_plan->N[1],ths->act_nfft_plan->N[2]); fflush(stdout);*/
1366 
1367  ths->act_nfft_plan->N_total=ths->act_nfft_plan->N[0]*ths->act_nfft_plan->N[1]*ths->act_nfft_plan->N[2];
1368  ths->act_nfft_plan->n[0]=ths->sigma*ths->act_nfft_plan->N[0];
1369  ths->act_nfft_plan->n[1]=ths->sigma*ths->act_nfft_plan->N[1];
1370  ths->act_nfft_plan->n[2]=ths->sigma*ths->act_nfft_plan->N[2];
1371  ths->act_nfft_plan->n_total=ths->act_nfft_plan->n[0]*ths->act_nfft_plan->n[1]*ths->act_nfft_plan->n[2];
1372 
1373  /* only for right - rear - top */
1374  if((J==0)||((J==1)&&(rr==1)))
1375  temp=-2.0*KPI;
1376  else
1377  temp=-3.0*KPI*X(exp2i)(J-rr);
1378 
1379  /* right */
1380  ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*0;
1381 
1382  for (j=0; j<M; j++)
1383  ths->act_nfft_plan->f[j]= ths->f[j] *
1384  cexp( - I*temp*ths->act_nfft_plan->x[3*j+0]);
1385 
1386  if(a>b)
1387  RSWAP(ths->act_nfft_plan->x,ths->x_120);
1388 
1389  if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1390  if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1391  if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1392  nfft_adjoint_direct(ths->act_nfft_plan);
1393  else
1394  short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1395  else
1396  short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1397  else
1398  nfft_adjoint(ths->act_nfft_plan);
1399 
1400  if(a>b)
1401  RSWAP(ths->act_nfft_plan->x,ths->x_120);
1402 
1403  /* rear */
1404  ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*1;
1405 
1406  for (j=0; j<M; j++)
1407  ths->act_nfft_plan->f[j]= ths->f[j] *
1408  cexp( - I*temp*ths->act_nfft_plan->x[3*j+1]);
1409 
1410  if(a>b)
1411  RSWAP(ths->act_nfft_plan->x,ths->x_021);
1412  if(a<b)
1413  RSWAP(ths->act_nfft_plan->x,ths->x_102);
1414 
1415  if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1416  if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1417  if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1418  nfft_adjoint_direct(ths->act_nfft_plan);
1419  else
1420  short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1421  else
1422  short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1423  else
1424  nfft_adjoint(ths->act_nfft_plan);
1425 
1426  if(a>b)
1427  RSWAP(ths->act_nfft_plan->x,ths->x_021);
1428  if(a<b)
1429  RSWAP(ths->act_nfft_plan->x,ths->x_102);
1430 
1431  /* top */
1432  ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*2;
1433 
1434  for (j=0; j<M; j++)
1435  ths->act_nfft_plan->f[j]= ths->f[j] *
1436  cexp( - I*temp*ths->act_nfft_plan->x[3*j+2]);
1437 
1438  if(a<b)
1439  RSWAP(ths->act_nfft_plan->x,ths->x_201);
1440 
1441  if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1442  if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1443  if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1444  nfft_adjoint_direct(ths->act_nfft_plan);
1445  else
1446  short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1447  else
1448  short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1449  else
1450  nfft_adjoint(ths->act_nfft_plan);
1451 
1452  if(a<b)
1453  RSWAP(ths->act_nfft_plan->x,ths->x_201);
1454 
1455  /* only for left - front - bottom */
1456  if((J==0)||((J==1)&&(rr==1)))
1457  temp=-4.0*KPI;
1458  else
1459  temp=-3.0*KPI*X(exp2i)(J-rr);
1460 
1461  /* left */
1462  ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*3;
1463 
1464  for (j=0; j<M; j++)
1465  ths->act_nfft_plan->f[j]= ths->f[j] *
1466  cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+0]);
1467 
1468  if(a>b)
1469  RSWAP(ths->act_nfft_plan->x,ths->x_120);
1470 
1471  if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1472  if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1473  if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1474  nfft_adjoint_direct(ths->act_nfft_plan);
1475  else
1476  short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1477  else
1478  short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1479  else
1480  nfft_adjoint(ths->act_nfft_plan);
1481 
1482  if(a>b)
1483  RSWAP(ths->act_nfft_plan->x,ths->x_120);
1484 
1485  /* front */
1486  ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*4;
1487 
1488  for (j=0; j<M; j++)
1489  ths->act_nfft_plan->f[j]= ths->f[j] *
1490  cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+1]);
1491 
1492  if(a>b)
1493  RSWAP(ths->act_nfft_plan->x,ths->x_021);
1494  if(a<b)
1495  RSWAP(ths->act_nfft_plan->x,ths->x_102);
1496 
1497  if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1498  if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1499  if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1500  nfft_adjoint_direct(ths->act_nfft_plan);
1501  else
1502  short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1503  else
1504  short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1505  else
1506  nfft_adjoint(ths->act_nfft_plan);
1507 
1508  if(a>b)
1509  RSWAP(ths->act_nfft_plan->x,ths->x_021);
1510  if(a<b)
1511  RSWAP(ths->act_nfft_plan->x,ths->x_102);
1512 
1513  /* bottom */
1514  ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*5;
1515 
1516  for (j=0; j<M; j++)
1517  ths->act_nfft_plan->f[j]= ths->f[j] *
1518  cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+2]);
1519 
1520  if(a<b)
1521  RSWAP(ths->act_nfft_plan->x,ths->x_201);
1522 
1523  if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1524  if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1525  if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1526  nfft_adjoint_direct(ths->act_nfft_plan);
1527  else
1528  short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1529  else
1530  short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1531  else
1532  nfft_adjoint(ths->act_nfft_plan);
1533 
1534  if(a<b)
1535  RSWAP(ths->act_nfft_plan->x,ths->x_201);
1536 
1537  sum_N_B_less_r+=6*N_B_r;
1538  } /* for(rr) */
1539 } /* void nsfft_adjoint_3d */
1540 
1541 void nsfft_trafo(nsfft_plan *ths)
1542 {
1543  if(ths->d==2)
1544  nsfft_trafo_2d(ths);
1545  else
1546  nsfft_trafo_3d(ths);
1547 }
1548 
1549 void nsfft_adjoint(nsfft_plan *ths)
1550 {
1551  if(ths->d==2)
1552  nsfft_adjoint_2d(ths);
1553  else
1554  nsfft_adjoint_3d(ths);
1555 }
1556 
1557 
1558 /*========================================================*/
1559 /* J >1, no precomputation at all!! */
1560 #ifdef GAUSSIAN
1561 static void nsfft_init_2d(nsfft_plan *ths, int J, int M, int m, unsigned snfft_flags)
1562 {
1563  int r;
1564  int N[2];
1565  int n[2];
1566 
1567  ths->flags=snfft_flags;
1568  ths->sigma=2;
1569  ths->J=J;
1570  ths->M_total=M;
1571  ths->N_total=(J+4)*X(exp2i)(J+1);
1572 
1573  /* memory allocation */
1574  ths->f = (double _Complex *)nfft_malloc(M*sizeof(double _Complex));
1575  ths->f_hat = (double _Complex *)nfft_malloc(ths->N_total*sizeof(double _Complex));
1576  ths->x_transposed= (double*)nfft_malloc(2*M*sizeof(double));
1577 
1578  ths->act_nfft_plan = (nfft_plan*)nfft_malloc(sizeof(nfft_plan));
1579  ths->center_nfft_plan = (nfft_plan*)nfft_malloc(sizeof(nfft_plan));
1580 
1581  ths->set_fftw_plan1=(fftw_plan*) nfft_malloc((J/2+1)*sizeof(fftw_plan));
1582  ths->set_fftw_plan2=(fftw_plan*) nfft_malloc((J/2+1)*sizeof(fftw_plan));
1583 
1584  ths->set_nfft_plan_1d = (nfft_plan*) nfft_malloc((X(log2i)(m)+1)*(sizeof(nfft_plan)));
1585 
1586  /* planning the small nffts */
1587  /* r=0 */
1588  N[0]=1; n[0]=ths->sigma*N[0];
1589  N[1]=X(exp2i)(J); n[1]=ths->sigma*N[1];
1590 
1591  nfft_init_guru(ths->act_nfft_plan,2,N,M,n,m, FG_PSI| MALLOC_X| MALLOC_F| FFTW_INIT, FFTW_MEASURE);
1592 
1593  if(ths->act_nfft_plan->flags & PRE_ONE_PSI)
1594  nfft_precompute_one_psi(ths->act_nfft_plan);
1595 
1596  ths->set_fftw_plan1[0]=ths->act_nfft_plan->my_fftw_plan1;
1597  ths->set_fftw_plan2[0]=ths->act_nfft_plan->my_fftw_plan2;
1598 
1599  for(r=1;r<=J/2;r++)
1600  {
1601  N[0]=X(exp2i)(r); n[0]=ths->sigma*N[0];
1602  N[1]=X(exp2i)(J-r); n[1]=ths->sigma*N[1];
1603  ths->set_fftw_plan1[r] =
1604  fftw_plan_dft(2, n, ths->act_nfft_plan->g1, ths->act_nfft_plan->g2,
1605  FFTW_FORWARD, ths->act_nfft_plan->fftw_flags);
1606 
1607  ths->set_fftw_plan2[r] =
1608  fftw_plan_dft(2, n, ths->act_nfft_plan->g2, ths->act_nfft_plan->g1,
1609  FFTW_BACKWARD, ths->act_nfft_plan->fftw_flags);
1610  }
1611 
1612  /* planning the 1d nffts */
1613  for(r=0;r<=X(log2i)(m);r++)
1614  {
1615  N[0]=X(exp2i)(J-r); n[0]=ths->sigma*N[0]; /* ==N[1] of the 2 dimensional plan */
1616 
1617  nfft_init_guru(&(ths->set_nfft_plan_1d[r]),1,N,M,n,m, MALLOC_X| MALLOC_F| FFTW_INIT, FFTW_MEASURE);
1618  ths->set_nfft_plan_1d[r].flags = ths->set_nfft_plan_1d[r].flags | FG_PSI;
1619  ths->set_nfft_plan_1d[r].K=ths->act_nfft_plan->K;
1620  ths->set_nfft_plan_1d[r].psi=ths->act_nfft_plan->psi;
1621  }
1622 
1623  /* center plan */
1624  /* J/2 == floor(((double)J) / 2.0) */
1625  N[0]=X(exp2i)(J/2+1); n[0]=ths->sigma*N[0];
1626  N[1]=X(exp2i)(J/2+1); n[1]=ths->sigma*N[1];
1627  nfft_init_guru(ths->center_nfft_plan,2,N,M,n, m, MALLOC_F| FFTW_INIT,
1628  FFTW_MEASURE);
1629  ths->center_nfft_plan->x= ths->act_nfft_plan->x;
1630  ths->center_nfft_plan->flags = ths->center_nfft_plan->flags|
1631  FG_PSI;
1632  ths->center_nfft_plan->K=ths->act_nfft_plan->K;
1633  ths->center_nfft_plan->psi=ths->act_nfft_plan->psi;
1634 
1635  if(ths->flags & NSDFT)
1636  {
1637  ths->index_sparse_to_full=(int*)nfft_malloc(ths->N_total*sizeof(int));
1638  init_index_sparse_to_full_2d(ths);
1639  }
1640 }
1641 #endif
1642 
1643 /*========================================================*/
1644 /* J >1, no precomputation at all!! */
1645 #ifdef GAUSSIAN
1646 static void nsfft_init_3d(nsfft_plan *ths, int J, int M, int m, unsigned snfft_flags)
1647 {
1648  int r,rr,a,b;
1649  int N[3];
1650  int n[3];
1651 
1652  ths->flags=snfft_flags;
1653  ths->sigma=2;
1654  ths->J=J;
1655  ths->M_total=M;
1656  ths->N_total=6*X(exp2i)(J)*(X(exp2i)((J+1)/2+1)-1)+X(exp2i)(3*(J/2+1));
1657 
1658  /* memory allocation */
1659  ths->f = (double _Complex *)nfft_malloc(M*sizeof(double _Complex));
1660  ths->f_hat = (double _Complex *)nfft_malloc(ths->N_total*sizeof(double _Complex));
1661 
1662  ths->x_102= (double*)nfft_malloc(3*M*sizeof(double));
1663  ths->x_201= (double*)nfft_malloc(3*M*sizeof(double));
1664  ths->x_120= (double*)nfft_malloc(3*M*sizeof(double));
1665  ths->x_021= (double*)nfft_malloc(3*M*sizeof(double));
1666 
1667  ths->act_nfft_plan = (nfft_plan*)nfft_malloc(sizeof(nfft_plan));
1668  ths->center_nfft_plan = (nfft_plan*)nfft_malloc(sizeof(nfft_plan));
1669 
1670  ths->set_fftw_plan1=(fftw_plan*) nfft_malloc(((J+1)/2+1)*sizeof(fftw_plan));
1671  ths->set_fftw_plan2=(fftw_plan*) nfft_malloc(((J+1)/2+1)*sizeof(fftw_plan));
1672 
1673  ths->set_nfft_plan_1d = (nfft_plan*) nfft_malloc((X(log2i)(m)+1)*(sizeof(nfft_plan)));
1674  ths->set_nfft_plan_2d = (nfft_plan*) nfft_malloc((X(log2i)(m)+1)*(sizeof(nfft_plan)));
1675 
1676  /* planning the small nffts */
1677  /* r=0 */
1678  N[0]=1; n[0]=ths->sigma*N[0];
1679  N[1]=1; n[1]=ths->sigma*N[1];
1680  N[2]=X(exp2i)(J); n[2]=ths->sigma*N[2];
1681 
1682  nfft_init_guru(ths->act_nfft_plan,3,N,M,n,m, FG_PSI| MALLOC_X| MALLOC_F, FFTW_MEASURE);
1683 
1684  if(ths->act_nfft_plan->flags & PRE_ONE_PSI)
1685  nfft_precompute_one_psi(ths->act_nfft_plan);
1686 
1687  /* malloc g1, g2 for maximal size */
1688  ths->act_nfft_plan->g1 = nfft_malloc(ths->sigma*ths->sigma*ths->sigma*X(exp2i)(J+(J+1)/2)*sizeof(double _Complex));
1689  ths->act_nfft_plan->g2 = nfft_malloc(ths->sigma*ths->sigma*ths->sigma*X(exp2i)(J+(J+1)/2)*sizeof(double _Complex));
1690 
1691  ths->act_nfft_plan->my_fftw_plan1 =
1692  fftw_plan_dft(3, n, ths->act_nfft_plan->g1, ths->act_nfft_plan->g2,
1693  FFTW_FORWARD, ths->act_nfft_plan->fftw_flags);
1694  ths->act_nfft_plan->my_fftw_plan2 =
1695  fftw_plan_dft(3, n, ths->act_nfft_plan->g2, ths->act_nfft_plan->g1,
1696  FFTW_BACKWARD, ths->act_nfft_plan->fftw_flags);
1697 
1698  ths->set_fftw_plan1[0]=ths->act_nfft_plan->my_fftw_plan1;
1699  ths->set_fftw_plan2[0]=ths->act_nfft_plan->my_fftw_plan2;
1700 
1701  for(rr=1;rr<=(J+1)/2;rr++)
1702  {
1703  a=X(exp2i)(J-rr);
1704  b=X(exp2i)(rr);
1705 
1706  r=MIN(rr,J-rr);
1707 
1708  n[0]=ths->sigma*X(exp2i)(r);
1709  if(a<b)
1710  n[1]=ths->sigma*X(exp2i)(J-r);
1711  else
1712  n[1]=ths->sigma*X(exp2i)(r);
1713  n[2]=ths->sigma*X(exp2i)(J-r);
1714 
1715  ths->set_fftw_plan1[rr] =
1716  fftw_plan_dft(3, n, ths->act_nfft_plan->g1, ths->act_nfft_plan->g2,
1717  FFTW_FORWARD, ths->act_nfft_plan->fftw_flags);
1718  ths->set_fftw_plan2[rr] =
1719  fftw_plan_dft(3, n, ths->act_nfft_plan->g2, ths->act_nfft_plan->g1,
1720  FFTW_BACKWARD, ths->act_nfft_plan->fftw_flags);
1721  }
1722 
1723  /* planning the 1d nffts */
1724  for(r=0;r<=X(log2i)(m);r++)
1725  {
1726  N[0]=X(exp2i)(J-r); n[0]=ths->sigma*N[0];
1727  N[1]=X(exp2i)(J-r); n[1]=ths->sigma*N[1];
1728 
1729  if(N[0]>m)
1730  {
1731  nfft_init_guru(&(ths->set_nfft_plan_1d[r]),1,N,M,n,m, MALLOC_X| MALLOC_F| FFTW_INIT, FFTW_MEASURE);
1732  ths->set_nfft_plan_1d[r].flags = ths->set_nfft_plan_1d[r].flags | FG_PSI;
1733  ths->set_nfft_plan_1d[r].K=ths->act_nfft_plan->K;
1734  ths->set_nfft_plan_1d[r].psi=ths->act_nfft_plan->psi;
1735  nfft_init_guru(&(ths->set_nfft_plan_2d[r]),2,N,M,n,m, MALLOC_X| MALLOC_F| FFTW_INIT, FFTW_MEASURE);
1736  ths->set_nfft_plan_2d[r].flags = ths->set_nfft_plan_2d[r].flags | FG_PSI;
1737  ths->set_nfft_plan_2d[r].K=ths->act_nfft_plan->K;
1738  ths->set_nfft_plan_2d[r].psi=ths->act_nfft_plan->psi;
1739  }
1740  }
1741 
1742  /* center plan */
1743  /* J/2 == floor(((double)J) / 2.0) */
1744  N[0]=X(exp2i)(J/2+1); n[0]=ths->sigma*N[0];
1745  N[1]=X(exp2i)(J/2+1); n[1]=ths->sigma*N[1];
1746  N[2]=X(exp2i)(J/2+1); n[2]=ths->sigma*N[2];
1747  nfft_init_guru(ths->center_nfft_plan,3,N,M,n, m, MALLOC_F| FFTW_INIT,
1748  FFTW_MEASURE);
1749  ths->center_nfft_plan->x= ths->act_nfft_plan->x;
1750  ths->center_nfft_plan->flags = ths->center_nfft_plan->flags|
1751  FG_PSI;
1752  ths->center_nfft_plan->K=ths->act_nfft_plan->K;
1753  ths->center_nfft_plan->psi=ths->act_nfft_plan->psi;
1754 
1755  if(ths->flags & NSDFT)
1756  {
1757  ths->index_sparse_to_full=(int*)nfft_malloc(ths->N_total*sizeof(int));
1758  init_index_sparse_to_full_3d(ths);
1759  }
1760 }
1761 #endif
1762 
1763 #ifdef GAUSSIAN
1764 void nsfft_init(nsfft_plan *ths, int d, int J, int M, int m, unsigned flags)
1765 {
1766  ths->d=d;
1767 
1768  if(ths->d==2)
1769  nsfft_init_2d(ths, J, M, m, flags);
1770  else
1771  nsfft_init_3d(ths, J, M, m, flags);
1772 
1773 
1774  ths->mv_trafo = (void (*) (void* ))nsfft_trafo;
1775  ths->mv_adjoint = (void (*) (void* ))nsfft_adjoint;
1776 }
1777 #else
1778 void nsfft_init(nsfft_plan *ths, int d, int J, int M, int m, unsigned flags)
1779 {
1780  UNUSED(ths);
1781  UNUSED(d);
1782  UNUSED(J);
1783  UNUSED(M);
1784  UNUSED(m);
1785  UNUSED(flags);
1786  fprintf(stderr,
1787  "\nError in kernel/nsfft_init: require GAUSSIAN window function\n");
1788 }
1789 #endif
1790 
1791 static void nsfft_finalize_2d(nsfft_plan *ths)
1792 {
1793  int r;
1794 
1795  if(ths->flags & NSDFT)
1796  nfft_free(ths->index_sparse_to_full);
1797 
1798  /* center plan */
1799  ths->center_nfft_plan->flags = ths->center_nfft_plan->flags ^ FG_PSI;
1800  nfft_finalize(ths->center_nfft_plan);
1801 
1802  /* the 1d nffts */
1803  for(r=0;r<=X(log2i)(ths->act_nfft_plan->m);r++)
1804  {
1805  ths->set_nfft_plan_1d[r].flags =
1806  ths->set_nfft_plan_1d[r].flags ^ FG_PSI;
1807  nfft_finalize(&(ths->set_nfft_plan_1d[r]));
1808  }
1809 
1810  /* finalize the small nffts */
1811  ths->act_nfft_plan->my_fftw_plan2=ths->set_fftw_plan2[0];
1812  ths->act_nfft_plan->my_fftw_plan1=ths->set_fftw_plan1[0];
1813 
1814  for(r=1;r<=ths->J/2;r++)
1815  {
1816  fftw_destroy_plan(ths->set_fftw_plan2[r]);
1817  fftw_destroy_plan(ths->set_fftw_plan1[r]);
1818  }
1819 
1820  /* r=0 */
1821  nfft_finalize(ths->act_nfft_plan);
1822 
1823  nfft_free(ths->set_nfft_plan_1d);
1824 
1825  nfft_free(ths->set_fftw_plan2);
1826  nfft_free(ths->set_fftw_plan1);
1827 
1828  nfft_free(ths->x_transposed);
1829 
1830  nfft_free(ths->f_hat);
1831  nfft_free(ths->f);
1832 }
1833 
1834 static void nsfft_finalize_3d(nsfft_plan *ths)
1835 {
1836  int r;
1837 
1838  if(ths->flags & NSDFT)
1839  nfft_free(ths->index_sparse_to_full);
1840 
1841  /* center plan */
1842  ths->center_nfft_plan->flags = ths->center_nfft_plan->flags ^ FG_PSI;
1843  nfft_finalize(ths->center_nfft_plan);
1844 
1845  /* the 1d and 2d nffts */
1846  for(r=0;r<=X(log2i)(ths->act_nfft_plan->m);r++)
1847  {
1848  if(X(exp2i)(ths->J-r)>ths->act_nfft_plan->m)
1849  {
1850  ths->set_nfft_plan_2d[r].flags = ths->set_nfft_plan_2d[r].flags ^ FG_PSI;
1851  nfft_finalize(&(ths->set_nfft_plan_2d[r]));
1852  ths->set_nfft_plan_1d[r].flags = ths->set_nfft_plan_1d[r].flags ^ FG_PSI;
1853  nfft_finalize(&(ths->set_nfft_plan_1d[r]));
1854  }
1855  }
1856 
1857  /* finalize the small nffts */
1858  ths->act_nfft_plan->my_fftw_plan2=ths->set_fftw_plan2[0];
1859  ths->act_nfft_plan->my_fftw_plan1=ths->set_fftw_plan1[0];
1860 
1861  for(r=1;r<=(ths->J+1)/2;r++)
1862  {
1863  fftw_destroy_plan(ths->set_fftw_plan2[r]);
1864  fftw_destroy_plan(ths->set_fftw_plan1[r]);
1865  }
1866 
1867  /* r=0 */
1868  nfft_finalize(ths->act_nfft_plan);
1869 
1870  nfft_free(ths->set_nfft_plan_1d);
1871  nfft_free(ths->set_nfft_plan_2d);
1872 
1873  nfft_free(ths->set_fftw_plan2);
1874  nfft_free(ths->set_fftw_plan1);
1875 
1876  nfft_free(ths->x_102);
1877  nfft_free(ths->x_201);
1878  nfft_free(ths->x_120);
1879  nfft_free(ths->x_021);
1880 
1881  nfft_free(ths->f_hat);
1882  nfft_free(ths->f);
1883 }
1884 
1885 void nsfft_finalize(nsfft_plan *ths)
1886 {
1887  if(ths->d==2)
1888  nsfft_finalize_2d(ths);
1889  else
1890  nsfft_finalize_3d(ths);
1891 }
nfft_plan::f_hat
fftw_complex * f_hat
Fourier coefficients.
Definition: nfft3.h:192
nfft_plan::d
NFFT_INT d
Dimension (rank).
Definition: nfft3.h:192
MALLOC_X
#define MALLOC_X
Definition: nfft3.h:199
nsfft_plan::d
int d
dimension, rank; d = 2, 3
Definition: nfft3.h:475
nsfft_plan::flags
unsigned flags
flags for precomputation, malloc
Definition: nfft3.h:475
nfft_plan::n
NFFT_INT * n
Length of FFTW transforms.
Definition: nfft3.h:192
FG_PSI
#define FG_PSI
Definition: nfft3.h:194
nsfft_plan::sigma
int sigma
oversampling-factor
Definition: nfft3.h:475
nsfft_plan::N_total
NFFT_INT N_total
Total number of Fourier coefficients.
Definition: nfft3.h:475
nsfft_plan::x_021
double * x_021
coordinate exchanged nodes, d=3
Definition: nfft3.h:475
nfft_vrand_shifted_unit_double
void nfft_vrand_shifted_unit_double(double *x, const NFFT_INT n)
Inits a vector of random double numbers in .
nsfft_plan::set_nfft_plan_2d
nfft_plan * set_nfft_plan_2d
nfft plans for short nffts
Definition: nfft3.h:475
nfft_plan::g1
fftw_complex * g1
Input of fftw.
Definition: nfft3.h:192
nfft_adjoint
void nfft_adjoint(nfft_plan *ths)
nfft_plan::f
fftw_complex * f
Samples.
Definition: nfft3.h:192
nsfft_init
void nsfft_init(nsfft_plan *ths, int d, int J, int M, int m, unsigned flags)
Definition: nsfft.c:1778
nsfft_plan::center_nfft_plan
nfft_plan * center_nfft_plan
central nfft block
Definition: nfft3.h:475
nsfft_plan::M_total
NFFT_INT M_total
Total number of samples.
Definition: nfft3.h:475
nfft_free
void nfft_free(void *p)
nsfft_plan::f
fftw_complex * f
Samples.
Definition: nfft3.h:475
nfft_plan
data structure for an NFFT (nonequispaced fast Fourier transform) plan with double precision ...
Definition: nfft3.h:192
nfft_vrand_unit_complex
void nfft_vrand_unit_complex(fftw_complex *x, const NFFT_INT n)
Inits a vector of random complex numbers in .
FFTW_INIT
#define FFTW_INIT
Definition: nfft3.h:203
nsfft_plan::set_nfft_plan_1d
nfft_plan * set_nfft_plan_1d
nfft plans for short nffts
Definition: nfft3.h:475
nfft_plan::N_total
NFFT_INT N_total
Total number of Fourier coefficients.
Definition: nfft3.h:192
nsfft_cp
void nsfft_cp(nsfft_plan *ths, nfft_plan *ths_full_plan)
Definition: nsfft.c:599
MALLOC_F
#define MALLOC_F
Definition: nfft3.h:201
nsfft_plan::index_sparse_to_full
int * index_sparse_to_full
index conversation, overflow for d=3, J=9!
Definition: nfft3.h:475
nfft_plan::M_total
NFFT_INT M_total
Total number of samples.
Definition: nfft3.h:192
X
#define X(name)
Include header for C99 complex datatype.
Definition: fastsum.h:57
RSWAP
#define RSWAP(x, y)
Swap two vectors.
Definition: infft.h:143
UNUSED
#define UNUSED(x)
Dummy use of unused parameters to silence compiler warnings.
Definition: infft.h:1365
nfft_trafo
void nfft_trafo(nfft_plan *ths)
nfft_plan::n_total
NFFT_INT n_total
Combined total length of FFTW transforms.
Definition: nfft3.h:192
nfft_plan::K
NFFT_INT K
Number of equispaced samples of window function.
Definition: nfft3.h:192
nfft_malloc
void * nfft_malloc(size_t n)
nsfft_plan::x_transposed
double * x_transposed
coordinate exchanged nodes, d = 2
Definition: nfft3.h:475
nfft_finalize
void nfft_finalize(nfft_plan *ths)
nfft_precompute_one_psi
void nfft_precompute_one_psi(nfft_plan *ths)
nsfft_plan::mv_adjoint
void(* mv_adjoint)(void *)
Adjoint transform.
Definition: nfft3.h:475
nsfft_plan::J
int J
problem size, i.e., d=2: N_total=(J+4) 2^(J+1) d=3: N_total=2^J 6(2^((J+1)/2+1)-1)+2^(3(J/2+1)) ...
Definition: nfft3.h:475
nsfft_init_random_nodes_coeffs
void nsfft_init_random_nodes_coeffs(nsfft_plan *ths)
Definition: nsfft.c:723
nfft_plan::psi
double * psi
Precomputed data for the sparse matrix , size depends on precomputation scheme.
Definition: nfft3.h:192
nsfft_plan::mv_trafo
void(* mv_trafo)(void *)
Transform.
Definition: nfft3.h:475
nsfft_adjoint
void nsfft_adjoint(nsfft_plan *ths)
Definition: nsfft.c:1549
nfft3.h
Header file for the nfft3 library.
nfft_plan::x
double * x
Nodes in time/spatial domain, size is doubles.
Definition: nfft3.h:192
nsfft_finalize
void nsfft_finalize(nsfft_plan *ths)
Definition: nsfft.c:1885
nfft_plan::N
NFFT_INT * N
Multi-bandwidth.
Definition: nfft3.h:192
PRE_ONE_PSI
#define PRE_ONE_PSI
Definition: nfft3.h:206
nfft_plan::flags
unsigned flags
Flags for precomputation, (de)allocation, and FFTW usage, default setting is PRE_PHI_HUT | PRE_PSI | ...
Definition: nfft3.h:192
nfft_plan::g2
fftw_complex * g2
Output of fftw.
Definition: nfft3.h:192
NSDFT
#define NSDFT
Definition: nfft3.h:476
nfft_plan::fftw_flags
unsigned fftw_flags
Flags for the FFTW, default is FFTW_ESTIMATE | FFTW_DESTROY_INPUT.
Definition: nfft3.h:192
nsfft_trafo
void nsfft_trafo(nsfft_plan *ths)
Definition: nsfft.c:1541
nsfft_plan
data structure for an NSFFT (nonequispaced sparse fast Fourier transform) plan with double precision ...
Definition: nfft3.h:475
nsfft_plan::act_nfft_plan
nfft_plan * act_nfft_plan
current nfft block
Definition: nfft3.h:475
nsfft_plan::f_hat
fftw_complex * f_hat
Fourier coefficients.
Definition: nfft3.h:475
nfft_plan::m
NFFT_INT m
Cut-off parameter for window function.
Definition: nfft3.h:192
Generated on Thu Apr 12 2018 21:00:49 for NFFT by   doxygen 1.8.6