30 #define NSFTT_DISABLE_TEST
40 for(j=0;j<ths->M_total;j++)
43 plan_1d->x[j] = ths->x[ths->d * j + 1];
46 for(k0=0;k0<ths->N[0];k0++)
48 plan_1d->f_hat = ths->f_hat + k0*ths->N[1];
52 for(j=0;j<ths->M_total;j++)
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);
65 for(j=0;j<ths->M_total;j++)
66 plan_1d->x[j] = ths->x[ths->d * j + 1];
68 for(k0=0;k0<ths->N[0];k0++)
70 for(j=0;j<ths->M_total;j++)
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);
76 plan_1d->f_hat = ths->f_hat + k0*ths->N[1];
78 nfft_adjoint(plan_1d);
90 for(j=0;j<ths->M_total;j++)
93 plan_1d->x[j] = ths->x[ths->d * j + 2];
96 for(k0=0;k0<ths->N[0];k0++)
97 for(k1=0;k1<ths->N[1];k1++)
99 plan_1d->f_hat = ths->f_hat + (k0*ths->N[1]+k1)*ths->N[2];
103 for(j=0;j<ths->M_total;j++)
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);
117 for(j=0;j<ths->M_total;j++)
118 plan_1d->x[j] = ths->x[ths->d * j + 2];
120 for(k0=0;k0<ths->N[0];k0++)
121 for(k1=0;k1<ths->N[1];k1++)
123 for(j=0;j<ths->M_total;j++)
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);
130 plan_1d->f_hat = ths->f_hat + (k0*ths->N[1]+k1)*ths->N[2];
132 nfft_adjoint(plan_1d);
144 for(j=0;j<ths->M_total;j++)
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];
151 for(k0=0;k0<ths->N[0];k0++)
153 plan_2d->f_hat = ths->f_hat + k0*ths->N[1]*ths->N[2];
157 for(j=0;j<ths->M_total;j++)
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);
170 for(j=0;j<ths->M_total;j++)
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];
176 for(k0=0;k0<ths->N[0];k0++)
178 for(j=0;j<ths->M_total;j++)
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);
184 plan_2d->f_hat = ths->f_hat + k0*ths->N[1]*ths->N[2];
186 nfft_adjoint(plan_2d);
193 static int index_sparse_to_full_direct_2d(
int J,
int k)
201 int i, o, a, b,s,l,m1,m2;
204 if (k>=(J+4)*
X(exp2i)(J+1))
213 i=k-4*((J+1)/2+1)*N_B;
314 static inline int index_sparse_to_full_2d(
nsfft_plan *ths,
int k)
317 if( k < ths->N_total)
318 return ths->index_sparse_to_full[k];
323 #ifndef NSFTT_DISABLE_TEST
324 static int index_full_to_sparse_2d(
int J,
int k)
336 if ( (k1>=-(a/2)) && (k1<a/2) && (k2>=(-a/2)) && (k2<a/2) )
338 return(4*((J+1)/2+1)*N_B+(k1+a/2)*a+(k2+a/2));
341 for (r=0; r<=(J+1)/2; r++)
345 if ( (k1>=-(b/2)) && (k1<(b+1)/2) && (k2>=a) && (k2<2*a) )
348 return((4*r+0)*N_B+(k1+b/2)*a+(k2-a));
350 return((4*r+0)*N_B+(k2-a)*b+(k1+b/2));
352 else if ( (k1>=a) && (k1<2*a) && (k2>=-(b/2)) && (k2<(b+1)/2) )
355 return((4*r+1)*N_B+(k2+b/2)*a+(k1-a));
357 return((4*r+1)*N_B+(k1-a)*b+(k2+b/2));
359 else if ( (k1>=-(b/2)) && (k1<(b+1)/2) && (k2>=-2*a) && (k2<-a) )
362 return((4*r+2)*N_B+(k1+b/2)*a+(k2+2*a));
364 return((4*r+2)*N_B+(k2+2*a)*b+(k1+b/2));
366 else if ( (k1>=-2*a) && (k1<-a) && (k2>=-(b/2)) && (k2<(b+1)/2) )
369 return((4*r+3)*N_B+(k2+b/2)*a+(k1+2*a));
371 return((4*r+3)*N_B+(k1+2*a)*b+(k2+b/2));
380 static void init_index_sparse_to_full_2d(
nsfft_plan *ths)
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);
390 static inline int index_sparse_to_full_3d(
nsfft_plan *ths,
int k)
393 if( k < ths->N_total)
394 return ths->index_sparse_to_full[k];
400 #ifndef NSFTT_DISABLE_TEST
401 static int index_full_to_sparse_3d(
int J,
int k)
410 int k2=((k/N)%N)-N/2;
415 if((k1>=-(a/2)) && (k1<a/2) && (k2>=(-a/2)) && (k2<a/2) && (k3>=(-a/2)) &&
418 return(6*
X(exp2i)(J)*(
X(exp2i)((J+1)/2+1)-1)+((k1+a/2)*a+(k2+a/2))*a+
423 for (r=0; r<=(J+1)/2; r++)
431 if ((k1>=a) && (k1<2*a) && (k2>=-(b/2)) && (k2<(b+1)/2) &&
432 (k3>=-(b/2)) && (k3<(b+1)/2))
435 return sum_N_B_less_r+N_B_r*0 + ((k2+b/2)*b+k3+b/2)*a + (k1-a);
437 return sum_N_B_less_r+N_B_r*0 + ((k1-a)*b+(k2+b/2))*b + (k3+b/2);
439 else if ((k2>=a) && (k2<2*a) && (k1>=-(b/2)) && (k1<(b+1)/2) &&
440 (k3>=-(b/2)) && (k3<(b+1)/2))
443 return sum_N_B_less_r+N_B_r*1 + ((k1+b/2)*b+k3+b/2)*a + (k2-a);
445 return sum_N_B_less_r+N_B_r*1 + ((k1+b/2)*b+(k2-a))*a + (k3+b/2);
447 return sum_N_B_less_r+N_B_r*1 + ((k2-a)*b+(k1+b/2))*b + (k3+b/2);
449 else if ((k3>=a) && (k3<2*a) && (k1>=-(b/2)) && (k1<(b+1)/2) &&
450 (k2>=-(b/2)) && (k2<(b+1)/2))
453 return sum_N_B_less_r+N_B_r*2 + ((k1+b/2)*b+k2+b/2)*a + (k3-a);
455 return sum_N_B_less_r+N_B_r*2 + ((k3-a)*b+(k1+b/2))*b + (k2+b/2);
458 else if ((k1>=-2*a) && (k1<-a) && (k2>=-(b/2)) && (k2<(b+1)/2) &&
459 (k3>=-(b/2)) && (k3<(b+1)/2))
462 return sum_N_B_less_r+N_B_r*3 + ((k2+b/2)*b+k3+b/2)*a + (k1+2*a);
464 return sum_N_B_less_r+N_B_r*3 + ((k1+2*a)*b+(k2+b/2))*b + (k3+b/2);
466 else if ((k2>=-2*a) && (k2<-a) && (k1>=-(b/2)) && (k1<(b+1)/2) &&
467 (k3>=-(b/2)) && (k3<(b+1)/2))
470 return sum_N_B_less_r+N_B_r*4 + ((k1+b/2)*b+k3+b/2)*a + (k2+2*a);
472 return sum_N_B_less_r+N_B_r*4 + ((k1+b/2)*b+(k2+2*a))*a + (k3+b/2);
474 return sum_N_B_less_r+N_B_r*4 + ((k2+2*a)*b+(k1+b/2))*b + (k3+b/2);
476 else if ((k3>=-2*a) && (k3<-a) && (k1>=-(b/2)) && (k1<(b+1)/2) &&
477 (k2>=-(b/2)) && (k2<(b+1)/2))
480 return sum_N_B_less_r+N_B_r*5 + ((k1+b/2)*b+k2+b/2)*a + (k3+2*a);
482 return sum_N_B_less_r+N_B_r*5 + ((k3+2*a)*b+(k1+b/2))*b + (k2+b/2);
485 sum_N_B_less_r+=6*N_B_r;
493 static void init_index_sparse_to_full_3d(
nsfft_plan *ths)
497 int N=
X(exp2i)(ths->J+2);
498 int Nc=ths->center_nfft_plan->N[0];
500 for (k_s=0, r=0; r<=(ths->J+1)/2; r++)
502 a=
X(exp2i)(ths->J-r);
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;
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;
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;
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;
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;
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;
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;
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;
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;
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;
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;
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;
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;
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;
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;
604 memset(ths_full_plan->f_hat, 0, ths_full_plan->N_total*
sizeof(
double _Complex));
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];
611 memcpy(ths_full_plan->x,ths->act_nfft_plan->x,ths->M_total*ths->d*
sizeof(
double));
614 #ifndef NSFTT_DISABLE_TEST
622 const int N=ths_full_plan->N[0];
623 const int N_B=
X(exp2i)(J);
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);
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]));
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]));
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]));
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]));
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]));
683 #ifndef NSFTT_DISABLE_TEST
684 static void test_sparse_to_full_2d(
nsfft_plan* ths)
687 int N=
X(exp2i)(ths->J+2);
689 printf(
"N=%d\n\n",N);
694 k_S=index_full_to_sparse_2d(ths->J, k1*N+k2);
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]);
702 #ifndef NSFTT_DISABLE_TEST
703 static void test_sparse_to_full_3d(
nsfft_plan* ths)
706 int N=
X(exp2i)(ths->J+2);
708 printf(
"N=%d\n\n",N);
714 k_S=index_full_to_sparse_3d(ths->J, (k1*N+k2)*N+k3);
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]);
728 nfft_vrand_unit_complex(ths->f_hat, ths->N_total);
731 nfft_vrand_shifted_unit_double(ths->act_nfft_plan->x, ths->d * ths->M_total);
734 for(j=0;j<ths->M_total;j++)
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];
740 for(j=0;j<ths->M_total;j++)
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];
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];
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];
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];
764 int N=
X(exp2i)(ths->J+2);
766 memset(ths->f,0,ths->M_total*
sizeof(
double _Complex));
768 for(k_S=0;k_S<ths->N_total;k_S++)
770 k_L=ths->index_sparse_to_full[k_S];
774 for(j=0;j<ths->M_total;j++)
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);
788 int N=
X(exp2i)(ths->J+2);
791 memset(ths->f,0,ths->M_total*
sizeof(
double _Complex));
793 for(k_S=0;k_S<ths->N_total;k_S++)
795 k_L=ths->index_sparse_to_full[k_S];
801 for(j=0;j<ths->M_total;j++)
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);
824 int N=
X(exp2i)(ths->J+2);
826 memset(ths->f_hat,0,ths->N_total*
sizeof(
double _Complex));
828 for(k_S=0;k_S<ths->N_total;k_S++)
830 k_L=ths->index_sparse_to_full[k_S];
834 for(j=0;j<ths->M_total;j++)
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);
848 int N=
X(exp2i)(ths->J+2);
851 memset(ths->f_hat,0,ths->N_total*
sizeof(
double _Complex));
853 for(k_S=0;k_S<ths->N_total;k_S++)
855 k_L=ths->index_sparse_to_full[k_S];
861 for(j=0;j<ths->M_total;j++)
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);
875 nsdft_adjoint_2d(ths);
877 nsdft_adjoint_3d(ths);
889 ths->center_nfft_plan->f_hat=ths->f_hat+4*((J+1)/2+1)*
X(exp2i)(J);
891 if (ths->center_nfft_plan->N[0]<=ths->center_nfft_plan->m)
892 nfft_trafo_direct(ths->center_nfft_plan);
894 nfft_trafo(ths->center_nfft_plan);
897 ths->f[j] = ths->center_nfft_plan->f[j];
899 for(rr=0;rr<=(J+1)/2;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];
908 temp=-3.0*KPI*
X(exp2i)(J-rr);
911 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+0)*
X(exp2i)(J);
914 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
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);
920 short_nfft_trafo_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
922 nfft_trafo(ths->act_nfft_plan);
925 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
928 ths->f[j] += ths->act_nfft_plan->f[j] *
929 cexp( + _Complex_I*temp*ths->act_nfft_plan->x[2*j+1]);
932 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+1)*
X(exp2i)(J);
934 if((r==rr)&&(J-rr!=rr))
935 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
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);
941 short_nfft_trafo_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
943 nfft_trafo(ths->act_nfft_plan);
945 if((r==rr)&&(J-rr!=rr))
946 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
949 ths->f[j] += ths->act_nfft_plan->f[j] *
950 cexp( + _Complex_I*temp*ths->act_nfft_plan->x[2*j+0]);
953 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+2)*
X(exp2i)(J);
956 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
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);
962 short_nfft_trafo_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
964 nfft_trafo(ths->act_nfft_plan);
967 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
970 ths->f[j] += ths->act_nfft_plan->f[j] *
971 cexp( - I*temp*ths->act_nfft_plan->x[2*j+1]);
974 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+3)*
X(exp2i)(J);
976 if((r==rr)&&(J-rr!=rr))
977 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
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);
983 short_nfft_trafo_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
985 nfft_trafo(ths->act_nfft_plan);
987 if((r==rr)&&(J-rr!=rr))
988 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
991 ths->f[j] += ths->act_nfft_plan->f[j] *
992 cexp( - I*temp*ths->act_nfft_plan->x[2*j+0]);
1006 ths->center_nfft_plan->f[j] = ths->f[j];
1008 ths->center_nfft_plan->f_hat=ths->f_hat+4*((J+1)/2+1)*
X(exp2i)(J);
1010 if (ths->center_nfft_plan->N[0]<=ths->center_nfft_plan->m)
1011 nfft_adjoint_direct(ths->center_nfft_plan);
1013 nfft_adjoint(ths->center_nfft_plan);
1015 for(rr=0;rr<=(J+1)/2;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];
1024 temp=-3.0*KPI*
X(exp2i)(J-rr);
1027 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+0)*
X(exp2i)(J);
1030 ths->act_nfft_plan->f[j]= ths->f[j] *
1031 cexp( - I*temp*ths->act_nfft_plan->x[2*j+1]);
1034 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
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);
1040 short_nfft_adjoint_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1042 nfft_adjoint(ths->act_nfft_plan);
1045 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
1048 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+1)*
X(exp2i)(J);
1051 ths->act_nfft_plan->f[j]= ths->f[j] *
1052 cexp( - I*temp*ths->act_nfft_plan->x[2*j+0]);
1054 if((r==rr)&&(J-rr!=rr))
1055 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
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);
1061 short_nfft_adjoint_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1063 nfft_adjoint(ths->act_nfft_plan);
1065 if((r==rr)&&(J-rr!=rr))
1066 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
1069 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+2)*
X(exp2i)(J);
1072 ths->act_nfft_plan->f[j]= ths->f[j] *
1073 cexp( + _Complex_I*temp*ths->act_nfft_plan->x[2*j+1]);
1076 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
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);
1082 short_nfft_adjoint_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1084 nfft_adjoint(ths->act_nfft_plan);
1087 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
1090 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+3)*
X(exp2i)(J);
1093 ths->act_nfft_plan->f[j]= ths->f[j] *
1094 cexp( + _Complex_I*temp*ths->act_nfft_plan->x[2*j+0]);
1096 if((r==rr)&&(J-rr!=rr))
1097 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
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);
1103 short_nfft_adjoint_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1105 nfft_adjoint(ths->act_nfft_plan);
1107 if((r==rr)&&(J-rr!=rr))
1108 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
1116 int sum_N_B_less_r,N_B_r,a,b;
1122 ths->center_nfft_plan->f_hat=ths->f_hat+6*
X(exp2i)(J)*(
X(exp2i)((J+1)/2+1)-1);
1124 if (ths->center_nfft_plan->N[0]<=ths->center_nfft_plan->m)
1125 nfft_trafo_direct(ths->center_nfft_plan);
1127 nfft_trafo(ths->center_nfft_plan);
1130 ths->f[j] = ths->center_nfft_plan->f[j];
1133 for(rr=0;rr<=(J+1)/2;rr++)
1141 ths->act_nfft_plan->my_fftw_plan1 = ths->set_fftw_plan1[rr];
1143 ths->act_nfft_plan->N[0]=
X(exp2i)(r);
1145 ths->act_nfft_plan->N[1]=
X(exp2i)(J-r);
1147 ths->act_nfft_plan->N[1]=
X(exp2i)(r);
1148 ths->act_nfft_plan->N[2]=
X(exp2i)(J-r);
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];
1159 if((J==0)||((J==1)&&(rr==1)))
1162 temp=-3.0*KPI*
X(exp2i)(J-rr);
1165 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*0;
1168 RSWAP(ths->act_nfft_plan->x,ths->x_120);
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);
1175 short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1177 short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1179 nfft_trafo(ths->act_nfft_plan);
1182 RSWAP(ths->act_nfft_plan->x,ths->x_120);
1185 ths->f[j] += ths->act_nfft_plan->f[j] *
1186 cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+0]);
1189 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*1;
1192 RSWAP(ths->act_nfft_plan->x,ths->x_021);
1194 RSWAP(ths->act_nfft_plan->x,ths->x_102);
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);
1201 short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1203 short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1205 nfft_trafo(ths->act_nfft_plan);
1208 RSWAP(ths->act_nfft_plan->x,ths->x_021);
1210 RSWAP(ths->act_nfft_plan->x,ths->x_102);
1213 ths->f[j] += ths->act_nfft_plan->f[j] *
1214 cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+1]);
1217 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*2;
1220 RSWAP(ths->act_nfft_plan->x,ths->x_201);
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);
1227 short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1229 short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1231 nfft_trafo(ths->act_nfft_plan);
1234 RSWAP(ths->act_nfft_plan->x,ths->x_201);
1237 ths->f[j] += ths->act_nfft_plan->f[j] *
1238 cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+2]);
1241 if((J==0)||((J==1)&&(rr==1)))
1244 temp=-3.0*KPI*
X(exp2i)(J-rr);
1247 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*3;
1250 RSWAP(ths->act_nfft_plan->x,ths->x_120);
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);
1257 short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1259 short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1261 nfft_trafo(ths->act_nfft_plan);
1264 RSWAP(ths->act_nfft_plan->x,ths->x_120);
1267 ths->f[j] += ths->act_nfft_plan->f[j] *
1268 cexp( - I*temp*ths->act_nfft_plan->x[3*j+0]);
1271 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*4;
1274 RSWAP(ths->act_nfft_plan->x,ths->x_021);
1276 RSWAP(ths->act_nfft_plan->x,ths->x_102);
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);
1283 short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1285 short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1287 nfft_trafo(ths->act_nfft_plan);
1290 RSWAP(ths->act_nfft_plan->x,ths->x_021);
1292 RSWAP(ths->act_nfft_plan->x,ths->x_102);
1295 ths->f[j] += ths->act_nfft_plan->f[j] *
1296 cexp( - I*temp*ths->act_nfft_plan->x[3*j+1]);
1299 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*5;
1302 RSWAP(ths->act_nfft_plan->x,ths->x_201);
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);
1309 short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1311 short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1313 nfft_trafo(ths->act_nfft_plan);
1316 RSWAP(ths->act_nfft_plan->x,ths->x_201);
1319 ths->f[j] += ths->act_nfft_plan->f[j] *
1320 cexp( - I*temp*ths->act_nfft_plan->x[3*j+2]);
1322 sum_N_B_less_r+=6*N_B_r;
1326 static void nsfft_adjoint_3d(
nsfft_plan *ths)
1330 int sum_N_B_less_r,N_B_r,a,b;
1337 ths->center_nfft_plan->f[j] = ths->f[j];
1339 ths->center_nfft_plan->f_hat=ths->f_hat+6*
X(exp2i)(J)*(
X(exp2i)((J+1)/2+1)-1);
1341 if (ths->center_nfft_plan->N[0]<=ths->center_nfft_plan->m)
1342 nfft_adjoint_direct(ths->center_nfft_plan);
1344 nfft_adjoint(ths->center_nfft_plan);
1347 for(rr=0;rr<=(J+1)/2;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];
1358 ths->act_nfft_plan->N[0]=
X(exp2i)(r);
1360 ths->act_nfft_plan->N[1]=
X(exp2i)(J-r);
1362 ths->act_nfft_plan->N[1]=
X(exp2i)(r);
1363 ths->act_nfft_plan->N[2]=
X(exp2i)(J-r);
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];
1374 if((J==0)||((J==1)&&(rr==1)))
1377 temp=-3.0*KPI*
X(exp2i)(J-rr);
1380 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*0;
1383 ths->act_nfft_plan->f[j]= ths->f[j] *
1384 cexp( - I*temp*ths->act_nfft_plan->x[3*j+0]);
1387 RSWAP(ths->act_nfft_plan->x,ths->x_120);
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);
1394 short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1396 short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1398 nfft_adjoint(ths->act_nfft_plan);
1401 RSWAP(ths->act_nfft_plan->x,ths->x_120);
1404 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*1;
1407 ths->act_nfft_plan->f[j]= ths->f[j] *
1408 cexp( - I*temp*ths->act_nfft_plan->x[3*j+1]);
1411 RSWAP(ths->act_nfft_plan->x,ths->x_021);
1413 RSWAP(ths->act_nfft_plan->x,ths->x_102);
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);
1420 short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1422 short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1424 nfft_adjoint(ths->act_nfft_plan);
1427 RSWAP(ths->act_nfft_plan->x,ths->x_021);
1429 RSWAP(ths->act_nfft_plan->x,ths->x_102);
1432 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*2;
1435 ths->act_nfft_plan->f[j]= ths->f[j] *
1436 cexp( - I*temp*ths->act_nfft_plan->x[3*j+2]);
1439 RSWAP(ths->act_nfft_plan->x,ths->x_201);
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);
1446 short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1448 short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1450 nfft_adjoint(ths->act_nfft_plan);
1453 RSWAP(ths->act_nfft_plan->x,ths->x_201);
1456 if((J==0)||((J==1)&&(rr==1)))
1459 temp=-3.0*KPI*
X(exp2i)(J-rr);
1462 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*3;
1465 ths->act_nfft_plan->f[j]= ths->f[j] *
1466 cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+0]);
1469 RSWAP(ths->act_nfft_plan->x,ths->x_120);
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);
1476 short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1478 short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1480 nfft_adjoint(ths->act_nfft_plan);
1483 RSWAP(ths->act_nfft_plan->x,ths->x_120);
1486 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*4;
1489 ths->act_nfft_plan->f[j]= ths->f[j] *
1490 cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+1]);
1493 RSWAP(ths->act_nfft_plan->x,ths->x_021);
1495 RSWAP(ths->act_nfft_plan->x,ths->x_102);
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);
1502 short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1504 short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1506 nfft_adjoint(ths->act_nfft_plan);
1509 RSWAP(ths->act_nfft_plan->x,ths->x_021);
1511 RSWAP(ths->act_nfft_plan->x,ths->x_102);
1514 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*5;
1517 ths->act_nfft_plan->f[j]= ths->f[j] *
1518 cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+2]);
1521 RSWAP(ths->act_nfft_plan->x,ths->x_201);
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);
1528 short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1530 short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1532 nfft_adjoint(ths->act_nfft_plan);
1535 RSWAP(ths->act_nfft_plan->x,ths->x_201);
1537 sum_N_B_less_r+=6*N_B_r;
1544 nsfft_trafo_2d(ths);
1546 nsfft_trafo_3d(ths);
1552 nsfft_adjoint_2d(ths);
1554 nsfft_adjoint_3d(ths);
1561 static void nsfft_init_2d(
nsfft_plan *ths,
int J,
int M,
int m,
unsigned snfft_flags)
1567 ths->flags=snfft_flags;
1571 ths->N_total=(J+4)*
X(exp2i)(J+1);
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));
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));
1588 N[0]=1; n[0]=ths->sigma*N[0];
1589 N[1]=
X(exp2i)(J); n[1]=ths->sigma*N[1];
1594 nfft_precompute_one_psi(ths->act_nfft_plan);
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;
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);
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);
1613 for(r=0;r<=
X(log2i)(m);r++)
1615 N[0]=
X(exp2i)(J-r); n[0]=ths->sigma*N[0];
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;
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];
1629 ths->center_nfft_plan->x= ths->act_nfft_plan->x;
1630 ths->center_nfft_plan->flags = ths->center_nfft_plan->flags|
1632 ths->center_nfft_plan->K=ths->act_nfft_plan->K;
1633 ths->center_nfft_plan->psi=ths->act_nfft_plan->psi;
1635 if(ths->flags &
NSDFT)
1637 ths->index_sparse_to_full=(
int*)nfft_malloc(ths->N_total*
sizeof(
int));
1638 init_index_sparse_to_full_2d(ths);
1646 static void nsfft_init_3d(
nsfft_plan *ths,
int J,
int M,
int m,
unsigned snfft_flags)
1652 ths->flags=snfft_flags;
1656 ths->N_total=6*
X(exp2i)(J)*(
X(exp2i)((J+1)/2+1)-1)+
X(exp2i)(3*(J/2+1));
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));
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));
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));
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];
1685 nfft_precompute_one_psi(ths->act_nfft_plan);
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));
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);
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;
1701 for(rr=1;rr<=(J+1)/2;rr++)
1708 n[0]=ths->sigma*
X(exp2i)(r);
1710 n[1]=ths->sigma*
X(exp2i)(J-r);
1712 n[1]=ths->sigma*
X(exp2i)(r);
1713 n[2]=ths->sigma*
X(exp2i)(J-r);
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);
1724 for(r=0;r<=
X(log2i)(m);r++)
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];
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;
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;
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];
1749 ths->center_nfft_plan->x= ths->act_nfft_plan->x;
1750 ths->center_nfft_plan->flags = ths->center_nfft_plan->flags|
1752 ths->center_nfft_plan->K=ths->act_nfft_plan->K;
1753 ths->center_nfft_plan->psi=ths->act_nfft_plan->psi;
1755 if(ths->flags &
NSDFT)
1757 ths->index_sparse_to_full=(
int*)nfft_malloc(ths->N_total*
sizeof(
int));
1758 init_index_sparse_to_full_3d(ths);
1769 nsfft_init_2d(ths, J, M, m, flags);
1771 nsfft_init_3d(ths, J, M, m, flags);
1787 "\nError in kernel/nsfft_init: require GAUSSIAN window function\n");
1791 static void nsfft_finalize_2d(
nsfft_plan *ths)
1795 if(ths->flags &
NSDFT)
1796 nfft_free(ths->index_sparse_to_full);
1799 ths->center_nfft_plan->flags = ths->center_nfft_plan->flags ^
FG_PSI;
1800 nfft_finalize(ths->center_nfft_plan);
1803 for(r=0;r<=
X(log2i)(ths->act_nfft_plan->m);r++)
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]));
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];
1814 for(r=1;r<=ths->J/2;r++)
1816 fftw_destroy_plan(ths->set_fftw_plan2[r]);
1817 fftw_destroy_plan(ths->set_fftw_plan1[r]);
1821 nfft_finalize(ths->act_nfft_plan);
1823 nfft_free(ths->set_nfft_plan_1d);
1825 nfft_free(ths->set_fftw_plan2);
1826 nfft_free(ths->set_fftw_plan1);
1828 nfft_free(ths->x_transposed);
1830 nfft_free(ths->f_hat);
1834 static void nsfft_finalize_3d(
nsfft_plan *ths)
1838 if(ths->flags &
NSDFT)
1839 nfft_free(ths->index_sparse_to_full);
1842 ths->center_nfft_plan->flags = ths->center_nfft_plan->flags ^
FG_PSI;
1843 nfft_finalize(ths->center_nfft_plan);
1846 for(r=0;r<=
X(log2i)(ths->act_nfft_plan->m);r++)
1848 if(
X(exp2i)(ths->J-r)>ths->act_nfft_plan->m)
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]));
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];
1861 for(r=1;r<=(ths->J+1)/2;r++)
1863 fftw_destroy_plan(ths->set_fftw_plan2[r]);
1864 fftw_destroy_plan(ths->set_fftw_plan1[r]);
1868 nfft_finalize(ths->act_nfft_plan);
1870 nfft_free(ths->set_nfft_plan_1d);
1871 nfft_free(ths->set_nfft_plan_2d);
1873 nfft_free(ths->set_fftw_plan2);
1874 nfft_free(ths->set_fftw_plan1);
1876 nfft_free(ths->x_102);
1877 nfft_free(ths->x_201);
1878 nfft_free(ths->x_120);
1879 nfft_free(ths->x_021);
1881 nfft_free(ths->f_hat);
1888 nsfft_finalize_2d(ths);
1890 nsfft_finalize_3d(ths);
*We expand this macro for each supported precision * X
#define RSWAP(x, y)
Swap two vectors.
#define UNUSED(x)
Dummy use of unused parameters to silence compiler warnings.
void nsfft_init(nsfft_plan *ths, int d, int J, int M, int m, unsigned flags)
void nsfft_cp(nsfft_plan *ths, nfft_plan *ths_full_plan)
void nsfft_init_random_nodes_coeffs(nsfft_plan *ths)
void nsfft_adjoint(nsfft_plan *ths)
void nsfft_trafo(nsfft_plan *ths)
void nsfft_finalize(nsfft_plan *ths)
Internal header file for auxiliary definitions and functions.
Header file for the nfft3 library.