44 #define X(name) NFCT(name)
47 static inline INT intprod(
const INT *vec,
const INT a,
const INT d)
52 for (t = 0; t < d; t++)
59 #define BASE(x) COS(x)
62 #define FOURIER_TRAFO FFTW_REDFT00
63 #define FFTW_DEFAULT_FLAGS FFTW_ESTIMATE | FFTW_DESTROY_INPUT
65 #define NODE(p,r) (ths->x[(p) * ths->d + (r)])
67 #define MACRO_with_FG_PSI fg_psi[t][lj[t]]
68 #define MACRO_with_PRE_PSI ths->psi[(j * ths->d + t) * (2 * ths->m + 2) + lj[t]]
69 #define MACRO_without_PRE_PSI PHI((2 * NN(ths->n[t])), ((ths->x[(j) * ths->d + t]) \
70 - ((R)(lj[t] + u[t])) / (K(2.0) * ((R)NN(ths->n[t])))), t)
71 #define MACRO_compute_PSI PHI((2 * NN(ths->n[t])), (NODE(j,t) - ((R)(lj[t] + u[t])) / (K(2.0) * ((R)NN(ths->n[t])))), t)
88 void X(trafo_direct)(
const X(plan) *ths)
90 R *f_hat = (R*)ths->f_hat, *f = (R*)ths->f;
92 memset(f, 0, (
size_t)(ths->M_total) *
sizeof(R));
99 #pragma omp parallel for default(shared) private(j)
101 for (j = 0; j < ths->M_total; j++)
104 for (k_L = 0; k_L < ths->N_total; k_L++)
106 R omega = K2PI * ((R)(k_L + OFFSET)) * ths->x[j];
107 f[j] += f_hat[k_L] * BASE(omega);
116 #pragma omp parallel for default(shared) private(j)
118 for (j = 0; j < ths->M_total; j++)
120 R x[ths->d], omega, Omega[ths->d + 1];
121 INT t, t2, k_L, k[ths->d];
123 for (t = 0; t < ths->d; t++)
126 x[t] = K2PI * ths->x[j * ths->d + t];
127 Omega[t+1] = BASE(((R)(k[t])) * x[t]) * Omega[t];
129 omega = Omega[ths->d];
131 for (k_L = 0; k_L < ths->N_total; k_L++)
133 f[j] += f_hat[k_L] * omega;
135 for (t = ths->d - 1; (t >= 1) && (k[t] == (ths->N[t] - 1)); t--)
140 for (t2 = t; t2 < ths->d; t2++)
141 Omega[t2+1] = BASE(((R)(k[t2])) * x[t2]) * Omega[t2];
143 omega = Omega[ths->d];
150 void X(adjoint_direct)(
const X(plan) *ths)
152 R *f_hat = (R*)ths->f_hat, *f = (R*)ths->f;
154 memset(f_hat, 0, (
size_t)(ths->N_total) *
sizeof(R));
161 #pragma omp parallel for default(shared) private(k_L)
162 for (k_L = 0; k_L < ths->N_total; k_L++)
165 for (j = 0; j < ths->M_total; j++)
167 R omega = K2PI * ((R)(k_L + OFFSET)) * ths->x[j];
168 f_hat[k_L] += f[j] * BASE(omega);
173 for (j = 0; j < ths->M_total; j++)
176 for (k_L = 0; k_L < ths->N_total; k_L++)
178 R omega = K2PI * ((R)(k_L + OFFSET)) * ths->x[j];
179 f_hat[k_L] += f[j] * BASE(omega);
189 #pragma omp parallel for default(shared) private(j, k_L)
190 for (k_L = 0; k_L < ths->N_total; k_L++)
192 INT k[ths->d], k_temp, t;
196 for (t = ths->d - 1; t >= 0; t--)
198 k[t] = k_temp % ths->N[t];
202 for (j = 0; j < ths->M_total; j++)
205 for (t = 0; t < ths->d; t++)
206 omega *= BASE(K2PI * (k[t] + OFFSET) * ths->x[j * ths->d + t]);
207 f_hat[k_L] += f[j] * omega;
211 for (j = 0; j < ths->M_total; j++)
213 R x[ths->d], omega, Omega[ths->d+1];
214 INT t, t2, k[ths->d];
216 for (t = 0; t < ths->d; t++)
219 x[t] = K2PI * ths->x[j * ths->d + t];
220 Omega[t+1] = BASE(((R)(k[t])) * x[t]) * Omega[t];
222 omega = Omega[ths->d];
223 for (k_L = 0; k_L < ths->N_total; k_L++)
225 f_hat[k_L] += f[j] * omega;
227 for (t = ths->d-1; (t >= 1) && (k[t] == ths->N[t] - 1); t--)
232 for (t2 = t; t2 < ths->d; t2++)
233 Omega[t2+1] = BASE(((R)(k[t2])) * x[t2]) * Omega[t2];
235 omega = Omega[ths->d];
261 static inline void uo(
const X(plan) *ths,
const INT j, INT *up, INT *op,
264 const R xj = ths->x[j * ths->d + act_dim];
265 INT c = LRINT(xj * (2 * NN(ths->n[(act_dim)])));
267 (*up) = c - (ths->m);
268 (*op) = c + 1 + (ths->m);
271 #define MACRO_D_compute_A \
273 g_hat[kg_plain[ths->d]] = f_hat[k_L] * c_phi_inv_k[ths->d]; \
276 #define MACRO_D_compute_T \
278 f_hat[k_L] = g_hat[kg_plain[ths->d]] * c_phi_inv_k[ths->d]; \
281 #define MACRO_D_init_result_A memset(g_hat, 0, (size_t)(ths->n_total) * sizeof(R));
283 #define MACRO_D_init_result_T memset(f_hat, 0, (size_t)(ths->N_total) * sizeof(R));
285 #define MACRO_with_PRE_PHI_HUT ths->c_phi_inv[t][kg[t]]
287 #define MACRO_compute_PHI_HUT_INV (K(1.0) / (PHI_HUT((2 * NN(ths->n[t])), kg[t] + OFFSET, t)))
289 #define MACRO_init_k_ks \
291 for (t = 0; t < ths->d; t++) \
298 #define MACRO_update_c_phi_inv_k(what_kind, which_phi) \
300 for (t = i; t < ths->d; t++) \
302 MACRO_update_c_phi_inv_k_ ## what_kind(which_phi); \
303 kg_plain[t+1] = kg_plain[t] * ths->n[t] + kg[t]; \
307 #define MACRO_update_c_phi_inv_k_A(which_phi) \
309 c_phi_inv_k[t+1] = (kg[t] == 0 ? K(1.0) : K(0.5)) * c_phi_inv_k[t] * MACRO_ ## which_phi; \
312 #define MACRO_update_c_phi_inv_k_T(which_phi) \
314 c_phi_inv_k[t+1] = c_phi_inv_k[t] * MACRO_ ## which_phi; \
317 #define MACRO_count_k_ks \
322 while ((kg[i] == ths->N[i]) && (i > 0)) \
331 #define MACRO_D(which_one) \
332 static inline void D_ ## which_one (X(plan) *ths) \
335 R c_phi_inv_k[ths->d+1]; \
340 INT kg_plain[ths->d+1]; \
342 f_hat = (R*)ths->f_hat; g_hat = (R*)ths->g_hat; \
343 MACRO_D_init_result_ ## which_one; \
345 c_phi_inv_k[0] = K(1.0); \
350 if (ths->flags & PRE_PHI_HUT) \
352 for (k_L = 0; k_L < ths->N_total; k_L++) \
354 MACRO_update_c_phi_inv_k(which_one, with_PRE_PHI_HUT); \
355 MACRO_D_compute_ ## which_one; \
361 for (k_L = 0; k_L < ths->N_total; k_L++) \
363 MACRO_update_c_phi_inv_k(which_one,compute_PHI_HUT_INV); \
364 MACRO_D_compute_ ## which_one; \
374 #define MACRO_B_init_result_A memset(f, 0, (size_t)(ths->M_total) * sizeof(R));
375 #define MACRO_B_init_result_T memset(g, 0, (size_t)(ths->n_total) * sizeof(R));
377 #define MACRO_B_PRE_FULL_PSI_compute_A \
379 (*fj) += ths->psi[ix] * g[ths->psi_index_g[ix]]; \
382 #define MACRO_B_PRE_FULL_PSI_compute_T \
385 INT d = ths->psi_index_g[ix]; \
386 for (t = ths->d - 1; t >= 0; t--) \
388 INT m = d % ths->n[t]; \
389 if (m != 0 && m != ths->n[t] - 1) \
393 g[ths->psi_index_g[ix]] += factor * ths->psi[ix] * (*fj); \
396 #define MACRO_B_compute_A \
398 (*fj) += phi_prod[ths->d] * g[ll_plain[ths->d]]; \
401 #define MACRO_B_compute_T \
403 g[ll_plain[ths->d]] += phi_prod[ths->d] * (*fj); \
406 #define MACRO_init_uo_l_lj_t \
408 for (t2 = 0; t2 < ths->d; t2++) \
410 uo(ths, j, &u[t2], &o[t2], t2); \
415 (u[(t2)] % (2 * NN(ths->n[(t2)]))) + (2 * NN(ths->n[(t2)])); \
417 lg_offset[(t2)] = u[(t2)] % (2 * NN(ths->n[(t2)])); \
418 if (lg_offset[(t2)] > NN(ths->n[(t2)])) \
419 lg_offset[(t2)] = -(2 * NN(ths->n[(t2)]) - lg_offset[(t2)]); \
421 if (lg_offset[t2] <= 0) \
423 l[t2] = -lg_offset[t2]; \
428 l[t2] = +lg_offset[t2]; \
439 #define FOO_T ((l[t] == 0 || l[t] == ths->n[t] - 1) ? K(1.0) : K(0.5))
441 #define MACRO_update_phi_prod_ll_plain(which_one,which_psi) \
443 for (t = t2; t < ths->d; t++) \
445 phi_prod[t+1] = (FOO_ ## which_one) * phi_prod[t] * (MACRO_ ## which_psi); \
446 ll_plain[t+1] = ll_plain[t] * ths->n[t] + l[t]; \
450 #define MACRO_count_uo_l_lj_t \
453 if ((l[(ths->d-1)] == 0) || (l[(ths->d-1)] == NN(ths->n[(ths->d-1)]))) \
454 count_lg[(ths->d-1)] *= -1; \
457 l[(ths->d-1)] += count_lg[(ths->d-1)]; \
462 while ((lj[t2] == (2 * ths->m + 2)) && (t2 > 0)) \
469 if ((l[(t2 - 1)] == 0) || (l[(t2 - 1)] == NN(ths->n[(t2 - 1)]))) \
470 count_lg[(t2 - 1)] *= -1; \
472 l[(t2 - 1)] += count_lg[(t2 - 1)]; \
475 if (lg_offset[t2] <= 0) \
477 l[t2] = -lg_offset[t2]; \
482 l[t2] = +lg_offset[t2]; \
490 #define MACRO_B(which_one) \
491 static inline void B_ ## which_one (X(plan) *ths) \
494 INT u[ths->d], o[ths->d]; \
500 INT ll_plain[ths->d+1]; \
501 R phi_prod[ths->d+1]; \
505 R fg_psi[ths->d][2*ths->m+2]; \
506 R fg_exp_l[ths->d][2*ths->m+2]; \
508 R tmpEXP1, tmpEXP2, tmpEXP2sq, tmp1, tmp2, tmp3; \
511 INT ip_s = ths->K/(ths->m+2); \
512 INT lg_offset[ths->d]; \
513 INT count_lg[ths->d]; \
515 f = (R*)ths->f; g = (R*)ths->g; \
517 MACRO_B_init_result_ ## which_one \
519 if (ths->flags & PRE_FULL_PSI) \
521 for (ix = 0, j = 0, fj = f; j < ths->M_total; j++, fj++) \
523 for (l_L = 0; l_L < ths->psi_index_f[j]; l_L++, ix++) \
525 MACRO_B_PRE_FULL_PSI_compute_ ## which_one; \
531 phi_prod[0] = K(1.0); \
534 for (t = 0, lprod = 1; t < ths->d; t++) \
535 lprod *= (2 * ths->m + 2); \
537 if (ths->flags & PRE_PSI) \
539 for (j = 0, fj = f; j < ths->M_total; j++, fj++) \
541 MACRO_init_uo_l_lj_t; \
543 for (l_L = 0; l_L < lprod; l_L++) \
545 MACRO_update_phi_prod_ll_plain(which_one, with_PRE_PSI); \
547 MACRO_B_compute_ ## which_one; \
549 MACRO_count_uo_l_lj_t; \
555 if (ths->flags & PRE_FG_PSI) \
557 for (t = 0; t < ths->d; t++) \
559 tmpEXP2 = EXP(K(-1.0) / ths->b[t]); \
560 tmpEXP2sq = tmpEXP2 * tmpEXP2; \
563 fg_exp_l[t][0] = K(1.0); \
565 for (lj_fg = 1; lj_fg <= (2 * ths->m + 2); lj_fg++) \
567 tmp3 = tmp2 * tmpEXP2; \
569 fg_exp_l[t][lj_fg] = fg_exp_l[t][lj_fg-1] * tmp3; \
573 for (j = 0, fj = f; j < ths->M_total; j++, fj++) \
575 MACRO_init_uo_l_lj_t; \
577 for (t = 0; t < ths->d; t++) \
579 fg_psi[t][0] = ths->psi[2 * (j * ths->d + t)]; \
580 tmpEXP1 = ths->psi[2 * (j * ths->d + t) + 1]; \
583 for (l_fg = u[t] + 1, lj_fg = 1; l_fg <= o[t]; l_fg++, lj_fg++) \
586 fg_psi[t][lj_fg] = fg_psi[t][0] * tmp1 * fg_exp_l[t][lj_fg]; \
590 for (l_L= 0; l_L < lprod; l_L++) \
592 MACRO_update_phi_prod_ll_plain(which_one, with_FG_PSI); \
594 MACRO_B_compute_ ## which_one; \
596 MACRO_count_uo_l_lj_t; \
602 if (ths->flags & FG_PSI) \
604 for (t = 0; t < ths->d; t++) \
606 tmpEXP2 = EXP(K(-1.0) / ths->b[t]); \
607 tmpEXP2sq = tmpEXP2 * tmpEXP2; \
610 fg_exp_l[t][0] = K(1.0); \
611 for (lj_fg = 1; lj_fg <= (2 * ths->m + 2); lj_fg++) \
613 tmp3 = tmp2 * tmpEXP2; \
615 fg_exp_l[t][lj_fg] = fg_exp_l[t][lj_fg-1] * tmp3; \
619 for (j = 0, fj = f; j < ths->M_total; j++, fj++) \
621 MACRO_init_uo_l_lj_t; \
623 for (t = 0; t < ths->d; t++) \
625 fg_psi[t][0] = (PHI((2 * NN(ths->n[t])), (ths->x[j*ths->d+t] - ((R)u[t])/(2 * NN(ths->n[t]))),(t)));\
627 tmpEXP1 = EXP(K(2.0) * ((2 * NN(ths->n[t])) * ths->x[j * ths->d + t] - u[t]) / ths->b[t]); \
629 for (l_fg = u[t] + 1, lj_fg = 1; l_fg <= o[t]; l_fg++, lj_fg++) \
632 fg_psi[t][lj_fg] = fg_psi[t][0] * tmp1 * fg_exp_l[t][lj_fg]; \
636 for (l_L = 0; l_L < lprod; l_L++) \
638 MACRO_update_phi_prod_ll_plain(which_one, with_FG_PSI); \
640 MACRO_B_compute_ ## which_one; \
642 MACRO_count_uo_l_lj_t; \
648 if (ths->flags & PRE_LIN_PSI) \
650 for (j = 0, fj = f; j < ths->M_total; j++, fj++) \
652 MACRO_init_uo_l_lj_t; \
654 for (t = 0; t < ths->d; t++) \
656 y[t] = (((2 * NN(ths->n[t])) * ths->x[j * ths->d + t] - (R)u[t]) \
657 * ((R)ths->K))/(ths->m + 2); \
658 ip_u = LRINT(FLOOR(y[t])); \
660 for (l_fg = u[t], lj_fg = 0; l_fg <= o[t]; l_fg++, lj_fg++) \
662 fg_psi[t][lj_fg] = ths->psi[(ths->K+1)*t + ABS(ip_u-lj_fg*ip_s)] \
663 * (1-ip_w) + ths->psi[(ths->K+1)*t + ABS(ip_u-lj_fg*ip_s+1)] \
668 for (l_L = 0; l_L < lprod; l_L++) \
670 MACRO_update_phi_prod_ll_plain(which_one, with_FG_PSI); \
672 MACRO_B_compute_ ## which_one; \
674 MACRO_count_uo_l_lj_t; \
681 for (j = 0, fj = &f[0]; j < ths->M_total; j++, fj += 1) \
683 MACRO_init_uo_l_lj_t; \
685 for (l_L = 0; l_L < lprod; l_L++) \
687 MACRO_update_phi_prod_ll_plain(which_one, without_PRE_PSI); \
689 MACRO_B_compute_ ## which_one; \
691 MACRO_count_uo_l_lj_t; \
702 void X(trafo)(
X(plan) *ths)
709 ths->g_hat = ths->g1;
722 FFTW(execute)(ths->my_fftw_r2r_plan);
743 void X(adjoint)(
X(plan) *ths)
750 ths->g_hat = ths->g2;
766 FFTW(execute)(ths->my_fftw_r2r_plan);
780 static inline
void precompute_phi_hut(
X(plan) *ths)
785 ths->c_phi_inv = (R**) Y(malloc)((size_t)(ths->d) *
sizeof(R*));
787 for (t = 0; t < ths->d; t++)
789 ths->c_phi_inv[t] = (R*)Y(malloc)((size_t)(ths->N[t] - OFFSET) *
sizeof(R));
791 for (ks[t] = 0; ks[t] < ths->N[t] - OFFSET; ks[t]++)
793 ths->c_phi_inv[t][ks[t]] = (K(1.0) / (PHI_HUT((2 * NN(ths->n[t])), ks[t] + OFFSET, t)));
803 void X(precompute_lin_psi)(
X(plan) *ths)
809 for (t = 0; t < ths->d; t++)
811 step = ((R)(ths->m+2)) / (((R)ths->K) * (2 * NN(ths->n[t])));
813 for (j = 0; j <= ths->K; j++)
815 ths->psi[(ths->K + 1) * t + j] = PHI((2 * NN(ths->n[t])), (j * step), t);
820 void X(precompute_fg_psi)(
X(plan) *ths)
827 for (t = 0; t < ths->d; t++)
831 for (j = 0; j < ths->M_total; j++)
833 uo(ths, j, &u, &o, t);
835 ths->psi[2 * (j*ths->d + t)] = (PHI((2 * NN(ths->n[t])),(ths->x[j * ths->d + t] - ((R)u) / (2 * NN(ths->n[t]))),(t)));
836 ths->psi[2 * (j*ths->d + t) + 1] = EXP(K(2.0) * ( (2 * NN(ths->n[t])) * ths->x[j * ths->d + t] - u) / ths->b[t]);
842 void X(precompute_psi)(
X(plan) *ths)
850 for (t = 0; t < ths->d; t++)
854 for (j = 0; j < ths->M_total; j++)
856 uo(ths, j, &u, &o, t);
858 for(lj = 0; lj < (2 * ths->m + 2); lj++)
859 ths->psi[(j * ths->d + t) * (2 * ths->m + 2) + lj] =
860 (PHI((2 * NN(ths->n[t])), ((ths->x[(j) * ths->d + (t)]) - ((R)(lj + u)) / (K(2.0) * ((R)NN(ths->n[t])))), t));
865 void X(precompute_full_psi)(
X(plan) *ths)
877 INT ll_plain[ths->d+1];
879 INT u[ths->d], o[ths->d];
880 INT count_lg[ths->d];
881 INT lg_offset[ths->d];
883 R phi_prod[ths->d+1];
889 phi_prod[0] = K(1.0);
892 for (t = 0, lprod = 1; t < ths->d; t++)
893 lprod *= 2 * ths->m + 2;
895 for (j = 0, ix = 0, ix_old = 0; j < ths->M_total; j++)
897 MACRO_init_uo_l_lj_t;
899 for (l_L = 0; l_L < lprod; l_L++, ix++)
901 MACRO_update_phi_prod_ll_plain(A, without_PRE_PSI);
903 ths->psi_index_g[ix] = ll_plain[ths->d];
904 ths->psi[ix] = phi_prod[ths->d];
906 MACRO_count_uo_l_lj_t;
909 ths->psi_index_f[j] = ix - ix_old;
915 void X(precompute_one_psi)(
X(plan) *ths)
918 X(precompute_psi)(ths);
920 X(precompute_full_psi)(ths);
922 X(precompute_fg_psi)(ths);
924 X(precompute_lin_psi)(ths);
927 static inline void init_help(
X(plan) *ths)
932 if (ths->flags & NFFT_OMP_BLOCKWISE_ADJOINT)
933 ths->flags |= NFFT_SORT_NODES;
935 ths->N_total = intprod(ths->N, OFFSET, ths->d);
936 ths->n_total = intprod(ths->n, 0, ths->d);
938 ths->sigma = (R*)Y(malloc)((size_t)(ths->d) *
sizeof(R));
940 for (t = 0; t < ths->d; t++)
941 ths->sigma[t] = ((R)NN(ths->n[t])) / ths->N[t];
944 ths->r2r_kind = (FFTW(r2r_kind)*)Y(malloc)((size_t)(ths->d) *
sizeof (FFTW(r2r_kind)));
945 for (t = 0; t < ths->d; t++)
946 ths->r2r_kind[t] = FOURIER_TRAFO;
951 ths->x = (R*)Y(malloc)((size_t)(ths->d * ths->M_total) *
sizeof(R));
954 ths->f_hat = (R*)Y(malloc)((size_t)(ths->N_total) *
sizeof(R));
957 ths->f = (R*)Y(malloc)((size_t)(ths->M_total) *
sizeof(R));
960 precompute_phi_hut(ths);
964 ths->K = (1U<< 10) * (ths->m+2);
965 ths->psi = (R*) Y(malloc)((size_t)((ths->K + 1) * ths->d) *
sizeof(R));
969 ths->psi = (R*) Y(malloc)((size_t)(ths->M_total * ths->d * 2) *
sizeof(R));
972 ths->psi = (R*) Y(malloc)((size_t)(ths->M_total * ths->d * (2 * ths->m + 2 )) *
sizeof(R));
976 for (t = 0, lprod = 1; t < ths->d; t++)
977 lprod *= 2 * ths->m + 2;
979 ths->psi = (R*) Y(malloc)((size_t)(ths->M_total * lprod) *
sizeof(R));
981 ths->psi_index_f = (INT*) Y(malloc)((size_t)(ths->M_total) *
sizeof(INT));
982 ths->psi_index_g = (INT*) Y(malloc)((size_t)(ths->M_total * lprod) *
sizeof(INT));
987 ths->g1 = (R*)Y(malloc)((size_t)(ths->n_total) *
sizeof(R));
990 ths->g2 = (R*) Y(malloc)((size_t)(ths->n_total) *
sizeof(R));
995 int *_n = Y(malloc)((size_t)(ths->d) *
sizeof(int));
997 for (t = 0; t < ths->d; t++)
998 _n[t] = (
int)(ths->n[t]);
1000 ths->my_fftw_r2r_plan = FFTW(plan_r2r)((int)ths->d, _n, ths->g1, ths->g2, ths->r2r_kind, ths->fftw_flags);
1010 ths->mv_trafo = (void (*) (
void* ))
X(trafo);
1011 ths->mv_adjoint = (void (*) (
void* ))
X(adjoint);
1014 void X(init)(
X(plan) *ths,
int d,
int *N,
int M_total)
1020 ths->N = (INT*) Y(malloc)((size_t)(d) *
sizeof(INT));
1022 for (t = 0; t < d; t++)
1023 ths->N[t] = (INT)N[t];
1025 ths->M_total = (INT)M_total;
1027 ths->n = (INT*) Y(malloc)((size_t)(d) *
sizeof(INT));
1029 for (t = 0; t < d; t++)
1030 ths->n[t] = 2 * (Y(next_power_of_2)(ths->N[t]) - 1) + OFFSET;
1032 ths->m = WINDOW_HELP_ESTIMATE_m;
1049 ths->fftw_flags = FFTW_ESTIMATE | FFTW_DESTROY_INPUT;
1054 void X(init_guru)(
X(plan) *ths,
int d,
int *N,
int M_total,
int *n,
int m,
1055 unsigned flags,
unsigned fftw_flags)
1060 ths->M_total = (INT)M_total;
1061 ths->N = (INT*)Y(malloc)((size_t)(ths->d) *
sizeof(INT));
1063 for (t = 0; t < d; t++)
1064 ths->N[t] = (INT)N[t];
1066 ths->n = (INT*)Y(malloc)((size_t)(ths->d) *
sizeof(INT));
1068 for (t = 0; t < d; t++)
1069 ths->n[t] = (INT)n[t];
1074 ths->fftw_flags = fftw_flags;
1079 void X(init_1d)(
X(plan) *ths,
int N1,
int M_total)
1085 X(init)(ths, 1, N, M_total);
1088 void X(init_2d)(
X(plan) *ths,
int N1,
int N2,
int M_total)
1095 X(init)(ths, 2, N, M_total);
1098 void X(init_3d)(
X(plan) *ths,
int N1,
int N2,
int N3,
int M_total)
1106 X(init)(ths, 3, N, M_total);
1109 const char*
X(check)(
X(plan) *ths)
1114 return "Member f not initialized.";
1117 return "Member x not initialized.";
1120 return "Member f_hat not initialized.";
1122 for (j = 0; j < ths->M_total * ths->d; j++)
1124 if ((ths->x[j] < K(0.0)) || (ths->x[j] >= K(0.5)))
1126 return "ths->x out of range [0.0,0.5)";
1130 for (j = 0; j < ths->d; j++)
1132 if (ths->sigma[j] <= 1)
1133 return "Oversampling factor too small";
1135 if(ths->N[j] - 1 <= ths->m)
1136 return "Polynomial degree N is smaller than cut-off m";
1138 if(ths->N[j]%2 == 1)
1139 return "polynomial degree N has to be even";
1144 void X(finalize)(
X(plan) *ths)
1154 #pragma omp critical (nfft_omp_critical_fftw_plan)
1156 FFTW(destroy_plan)(ths->my_fftw_r2r_plan);
1166 Y(free)(ths->psi_index_g);
1167 Y(free)(ths->psi_index_f);
1182 for (t = 0; t < ths->d; t++)
1183 Y(free)(ths->c_phi_inv[t]);
1184 Y(free)(ths->c_phi_inv);
1191 Y(free)(ths->f_hat);
1196 WINDOW_HELP_FINALIZE;
1200 Y(free)(ths->sigma);
1202 Y(free)(ths->r2r_kind);
#define TIC(a)
Timing, method works since the inaccurate timer is updated mostly in the measured function...
#define X(name)
Include header for C99 complex datatype.
Header file for the nfft3 library.