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