44 #define X(name) NFST(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) SIN(x)
62 #define FOURIER_TRAFO FFTW_RODFT00
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] = 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] = K(0.5) * c_phi_inv_k[t] * MACRO_ ## which_phi; \
317 #define MACRO_count_k_ks \
322 while ((kg[i] == ths->N[i] - 1) && (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 \
384 g[ths->psi_index_g[ix]] += ths->psi[ix] * (*fj); \
387 #define MACRO_B_compute_A \
389 (*fj) += phi_prod[ths->d] * g[ll_plain[ths->d]]; \
392 #define MACRO_B_compute_T \
394 g[ll_plain[ths->d]] += phi_prod[ths->d] * (*fj); \
397 #define MACRO_init_uo_l_lj_t \
399 for (t2 = 0; t2 < ths->d; t2++) \
401 uo(ths, j, &u[t2], &o[t2], t2); \
406 (u[(t2)] % (2 * NN(ths->n[(t2)]))) + (2 * NN(ths->n[(t2)])); \
408 lg_offset[(t2)] = u[(t2)] % (2 * NN(ths->n[(t2)])); \
409 if (lg_offset[(t2)] > NN(ths->n[(t2)])) \
410 lg_offset[(t2)] = -(2 * NN(ths->n[(t2)]) - lg_offset[(t2)]); \
412 if (lg_offset[t2] <= 0) \
414 l[t2] = -lg_offset[t2]; \
419 l[t2] = +lg_offset[t2]; \
428 #define FOO_A ((R)count_lg[t])
430 #define FOO_T ((R)count_lg[t])
432 #define MACRO_update_phi_prod_ll_plain(which_one,which_psi) \
434 for (t = t2; t < ths->d; t++) \
436 if ((l[t] != 0) && (l[t] != NN(ths->n[t]))) \
438 phi_prod[t+1] = (FOO_ ## which_one) * phi_prod[t] * (MACRO_ ## which_psi); \
439 ll_plain[t+1] = ll_plain[t] * ths->n[t] + l[t] - 1; \
443 phi_prod[t + 1] = K(0.0); \
444 ll_plain[t+1] = ll_plain[t] * ths->n[t]; \
449 #define MACRO_count_uo_l_lj_t \
452 if ((l[(ths->d-1)] == 0) || (l[(ths->d-1)] == NN(ths->n[(ths->d-1)]))) \
453 count_lg[(ths->d-1)] *= -1; \
456 l[(ths->d-1)] += count_lg[(ths->d-1)]; \
461 while ((lj[t2] == (2 * ths->m + 2)) && (t2 > 0)) \
468 if ((l[(t2 - 1)] == 0) || (l[(t2 - 1)] == NN(ths->n[(t2 - 1)]))) \
469 count_lg[(t2 - 1)] *= -1; \
471 l[(t2 - 1)] += count_lg[(t2 - 1)]; \
474 if (lg_offset[t2] <= 0) \
476 l[t2] = -lg_offset[t2]; \
481 l[t2] = +lg_offset[t2]; \
489 #define MACRO_B(which_one) \
490 static inline void B_ ## which_one (X(plan) *ths) \
493 INT u[ths->d], o[ths->d]; \
499 INT ll_plain[ths->d+1]; \
500 R phi_prod[ths->d+1]; \
504 R fg_psi[ths->d][2*ths->m+2]; \
505 R fg_exp_l[ths->d][2*ths->m+2]; \
507 R tmpEXP1, tmpEXP2, tmpEXP2sq, tmp1, tmp2, tmp3; \
510 INT ip_s = ths->K/(ths->m+2); \
511 INT lg_offset[ths->d]; \
512 INT count_lg[ths->d]; \
514 f = (R*)ths->f; g = (R*)ths->g; \
516 MACRO_B_init_result_ ## which_one \
518 if (ths->flags & PRE_FULL_PSI) \
520 for (ix = 0, j = 0, fj = f; j < ths->M_total; j++, fj++) \
522 for (l_L = 0; l_L < ths->psi_index_f[j]; l_L++, ix++) \
524 MACRO_B_PRE_FULL_PSI_compute_ ## which_one; \
530 phi_prod[0] = K(1.0); \
533 for (t = 0, lprod = 1; t < ths->d; t++) \
534 lprod *= (2 * ths->m + 2); \
536 if (ths->flags & PRE_PSI) \
538 for (j = 0, fj = f; j < ths->M_total; j++, fj++) \
540 MACRO_init_uo_l_lj_t; \
542 for (l_L = 0; l_L < lprod; l_L++) \
544 MACRO_update_phi_prod_ll_plain(which_one, with_PRE_PSI); \
546 MACRO_B_compute_ ## which_one; \
548 MACRO_count_uo_l_lj_t; \
554 if (ths->flags & PRE_FG_PSI) \
556 for (t = 0; t < ths->d; t++) \
558 tmpEXP2 = EXP(K(-1.0) / ths->b[t]); \
559 tmpEXP2sq = tmpEXP2 * tmpEXP2; \
562 fg_exp_l[t][0] = K(1.0); \
564 for (lj_fg = 1; lj_fg <= (2 * ths->m + 2); lj_fg++) \
566 tmp3 = tmp2 * tmpEXP2; \
568 fg_exp_l[t][lj_fg] = fg_exp_l[t][lj_fg-1] * tmp3; \
572 for (j = 0, fj = f; j < ths->M_total; j++, fj++) \
574 MACRO_init_uo_l_lj_t; \
576 for (t = 0; t < ths->d; t++) \
578 fg_psi[t][0] = ths->psi[2 * (j * ths->d + t)]; \
579 tmpEXP1 = ths->psi[2 * (j * ths->d + t) + 1]; \
582 for (l_fg = u[t] + 1, lj_fg = 1; l_fg <= o[t]; l_fg++, lj_fg++) \
585 fg_psi[t][lj_fg] = fg_psi[t][0] * tmp1 * fg_exp_l[t][lj_fg]; \
589 for (l_L= 0; l_L < lprod; l_L++) \
591 MACRO_update_phi_prod_ll_plain(which_one, with_FG_PSI); \
593 MACRO_B_compute_ ## which_one; \
595 MACRO_count_uo_l_lj_t; \
601 if (ths->flags & FG_PSI) \
603 for (t = 0; t < ths->d; t++) \
605 tmpEXP2 = EXP(K(-1.0) / ths->b[t]); \
606 tmpEXP2sq = tmpEXP2 * tmpEXP2; \
609 fg_exp_l[t][0] = K(1.0); \
610 for (lj_fg = 1; lj_fg <= (2 * ths->m + 2); lj_fg++) \
612 tmp3 = tmp2 * tmpEXP2; \
614 fg_exp_l[t][lj_fg] = fg_exp_l[t][lj_fg-1] * tmp3; \
618 for (j = 0, fj = f; j < ths->M_total; j++, fj++) \
620 MACRO_init_uo_l_lj_t; \
622 for (t = 0; t < ths->d; t++) \
624 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)));\
626 tmpEXP1 = EXP(K(2.0) * ((2 * NN(ths->n[t])) * ths->x[j * ths->d + t] - u[t]) / ths->b[t]); \
628 for (l_fg = u[t] + 1, lj_fg = 1; l_fg <= o[t]; l_fg++, lj_fg++) \
631 fg_psi[t][lj_fg] = fg_psi[t][0] * tmp1 * fg_exp_l[t][lj_fg]; \
635 for (l_L = 0; l_L < lprod; l_L++) \
637 MACRO_update_phi_prod_ll_plain(which_one, with_FG_PSI); \
639 MACRO_B_compute_ ## which_one; \
641 MACRO_count_uo_l_lj_t; \
647 if (ths->flags & PRE_LIN_PSI) \
649 for (j = 0, fj = f; j < ths->M_total; j++, fj++) \
651 MACRO_init_uo_l_lj_t; \
653 for (t = 0; t < ths->d; t++) \
655 y[t] = (((2 * NN(ths->n[t])) * ths->x[j * ths->d + t] - (R)u[t]) \
656 * ((R)ths->K))/(ths->m + 2); \
657 ip_u = LRINT(FLOOR(y[t])); \
659 for (l_fg = u[t], lj_fg = 0; l_fg <= o[t]; l_fg++, lj_fg++) \
661 fg_psi[t][lj_fg] = ths->psi[(ths->K+1)*t + ABS(ip_u-lj_fg*ip_s)] \
662 * (1-ip_w) + ths->psi[(ths->K+1)*t + ABS(ip_u-lj_fg*ip_s+1)] \
667 for (l_L = 0; l_L < lprod; l_L++) \
669 MACRO_update_phi_prod_ll_plain(which_one, with_FG_PSI); \
671 MACRO_B_compute_ ## which_one; \
673 MACRO_count_uo_l_lj_t; \
680 for (j = 0, fj = &f[0]; j < ths->M_total; j++, fj += 1) \
682 MACRO_init_uo_l_lj_t; \
684 for (l_L = 0; l_L < lprod; l_L++) \
686 MACRO_update_phi_prod_ll_plain(which_one, without_PRE_PSI); \
688 MACRO_B_compute_ ## which_one; \
690 MACRO_count_uo_l_lj_t; \
701 void X(trafo)(
X(plan) *ths)
708 ths->g_hat = ths->g1;
721 FFTW(execute)(ths->my_fftw_r2r_plan);
742 void X(adjoint)(
X(plan) *ths)
749 ths->g_hat = ths->g2;
765 FFTW(execute)(ths->my_fftw_r2r_plan);
779 static inline
void precompute_phi_hut(
X(plan) *ths)
784 ths->c_phi_inv = (R**) Y(malloc)((size_t)(ths->d) *
sizeof(R*));
786 for (t = 0; t < ths->d; t++)
788 ths->c_phi_inv[t] = (R*)Y(malloc)((size_t)(ths->N[t] - OFFSET) *
sizeof(R));
790 for (ks[t] = 0; ks[t] < ths->N[t] - OFFSET; ks[t]++)
792 ths->c_phi_inv[t][ks[t]] = (K(1.0) / (PHI_HUT((2 * NN(ths->n[t])), ks[t] + OFFSET, t)));
802 void X(precompute_lin_psi)(
X(plan) *ths)
808 for (t = 0; t < ths->d; t++)
810 step = ((R)(ths->m+2)) / (((R)ths->K) * (2 * NN(ths->n[t])));
812 for (j = 0; j <= ths->K; j++)
814 ths->psi[(ths->K + 1) * t + j] = PHI((2 * NN(ths->n[t])), (j * step), t);
819 void X(precompute_fg_psi)(
X(plan) *ths)
826 for (t = 0; t < ths->d; t++)
830 for (j = 0; j < ths->M_total; j++)
832 uo(ths, j, &u, &o, t);
834 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)));
835 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]);
841 void X(precompute_psi)(
X(plan) *ths)
849 for (t = 0; t < ths->d; t++)
853 for (j = 0; j < ths->M_total; j++)
855 uo(ths, j, &u, &o, t);
857 for(lj = 0; lj < (2 * ths->m + 2); lj++)
858 ths->psi[(j * ths->d + t) * (2 * ths->m + 2) + lj] =
859 (PHI((2 * NN(ths->n[t])), ((ths->x[(j) * ths->d + (t)]) - ((R)(lj + u)) / (K(2.0) * ((R)NN(ths->n[t])))), t));
864 void X(precompute_full_psi)(
X(plan) *ths)
876 INT ll_plain[ths->d+1];
878 INT u[ths->d], o[ths->d];
879 INT count_lg[ths->d];
880 INT lg_offset[ths->d];
882 R phi_prod[ths->d+1];
888 phi_prod[0] = K(1.0);
891 for (t = 0, lprod = 1; t < ths->d; t++)
892 lprod *= 2 * ths->m + 2;
894 for (j = 0, ix = 0, ix_old = 0; j < ths->M_total; j++)
896 MACRO_init_uo_l_lj_t;
898 for (l_L = 0; l_L < lprod; l_L++, ix++)
900 MACRO_update_phi_prod_ll_plain(A, without_PRE_PSI);
902 ths->psi_index_g[ix] = ll_plain[ths->d];
903 ths->psi[ix] = phi_prod[ths->d];
905 MACRO_count_uo_l_lj_t;
908 ths->psi_index_f[j] = ix - ix_old;
914 void X(precompute_one_psi)(
X(plan) *ths)
917 X(precompute_psi)(ths);
919 X(precompute_full_psi)(ths);
921 X(precompute_fg_psi)(ths);
923 X(precompute_lin_psi)(ths);
926 static inline void init_help(
X(plan) *ths)
931 if (ths->flags & NFFT_OMP_BLOCKWISE_ADJOINT)
932 ths->flags |= NFFT_SORT_NODES;
934 ths->N_total = intprod(ths->N, OFFSET, ths->d);
935 ths->n_total = intprod(ths->n, 0, ths->d);
937 ths->sigma = (R*)Y(malloc)((size_t)(ths->d) *
sizeof(R));
939 for (t = 0; t < ths->d; t++)
940 ths->sigma[t] = ((R)NN(ths->n[t])) / ths->N[t];
943 ths->r2r_kind = (FFTW(r2r_kind)*)Y(malloc)((size_t)(ths->d) *
sizeof (FFTW(r2r_kind)));
944 for (t = 0; t < ths->d; t++)
945 ths->r2r_kind[t] = FOURIER_TRAFO;
950 ths->x = (R*)Y(malloc)((size_t)(ths->d * ths->M_total) *
sizeof(R));
953 ths->f_hat = (R*)Y(malloc)((size_t)(ths->N_total) *
sizeof(R));
956 ths->f = (R*)Y(malloc)((size_t)(ths->M_total) *
sizeof(R));
959 precompute_phi_hut(ths);
963 ths->K = (1U<< 10) * (ths->m+2);
964 ths->psi = (R*) Y(malloc)((size_t)((ths->K + 1) * ths->d) *
sizeof(R));
968 ths->psi = (R*) Y(malloc)((size_t)(ths->M_total * ths->d * 2) *
sizeof(R));
971 ths->psi = (R*) Y(malloc)((size_t)(ths->M_total * ths->d * (2 * ths->m + 2 )) *
sizeof(R));
975 for (t = 0, lprod = 1; t < ths->d; t++)
976 lprod *= 2 * ths->m + 2;
978 ths->psi = (R*) Y(malloc)((size_t)(ths->M_total * lprod) *
sizeof(R));
980 ths->psi_index_f = (INT*) Y(malloc)((size_t)(ths->M_total) *
sizeof(INT));
981 ths->psi_index_g = (INT*) Y(malloc)((size_t)(ths->M_total * lprod) *
sizeof(INT));
986 ths->g1 = (R*)Y(malloc)((size_t)(ths->n_total) *
sizeof(R));
989 ths->g2 = (R*) Y(malloc)((size_t)(ths->n_total) *
sizeof(R));
994 int *_n = Y(malloc)((size_t)(ths->d) *
sizeof(int));
996 for (t = 0; t < ths->d; t++)
997 _n[t] = (
int)(ths->n[t]);
999 ths->my_fftw_r2r_plan = FFTW(plan_r2r)((int)ths->d, _n, ths->g1, ths->g2, ths->r2r_kind, ths->fftw_flags);
1009 ths->mv_trafo = (void (*) (
void* ))
X(trafo);
1010 ths->mv_adjoint = (void (*) (
void* ))
X(adjoint);
1013 void X(init)(
X(plan) *ths,
int d,
int *N,
int M_total)
1019 ths->N = (INT*) Y(malloc)((size_t)(d) *
sizeof(INT));
1021 for (t = 0; t < d; t++)
1022 ths->N[t] = (INT)N[t];
1024 ths->M_total = (INT)M_total;
1026 ths->n = (INT*) Y(malloc)((size_t)(d) *
sizeof(INT));
1028 for (t = 0; t < d; t++)
1029 ths->n[t] = 2 * (Y(next_power_of_2)(ths->N[t]) - 1) + OFFSET;
1031 ths->m = WINDOW_HELP_ESTIMATE_m;
1048 ths->fftw_flags = FFTW_ESTIMATE | FFTW_DESTROY_INPUT;
1053 void X(init_guru)(
X(plan) *ths,
int d,
int *N,
int M_total,
int *n,
int m,
1054 unsigned flags,
unsigned fftw_flags)
1059 ths->M_total = (INT)M_total;
1060 ths->N = (INT*)Y(malloc)((size_t)(ths->d) *
sizeof(INT));
1062 for (t = 0; t < d; t++)
1063 ths->N[t] = (INT)N[t];
1065 ths->n = (INT*)Y(malloc)((size_t)(ths->d) *
sizeof(INT));
1067 for (t = 0; t < d; t++)
1068 ths->n[t] = (INT)n[t];
1073 ths->fftw_flags = fftw_flags;
1078 void X(init_1d)(
X(plan) *ths,
int N1,
int M_total)
1084 X(init)(ths, 1, N, M_total);
1087 void X(init_2d)(
X(plan) *ths,
int N1,
int N2,
int M_total)
1094 X(init)(ths, 2, N, M_total);
1097 void X(init_3d)(
X(plan) *ths,
int N1,
int N2,
int N3,
int M_total)
1105 X(init)(ths, 3, N, M_total);
1108 const char*
X(check)(
X(plan) *ths)
1113 return "Member f not initialized.";
1116 return "Member x not initialized.";
1119 return "Member f_hat not initialized.";
1121 for (j = 0; j < ths->M_total * ths->d; j++)
1123 if ((ths->x[j] < K(0.0)) || (ths->x[j] >= K(0.5)))
1125 return "ths->x out of range [0.0,0.5)";
1129 for (j = 0; j < ths->d; j++)
1131 if (ths->sigma[j] <= 1)
1132 return "Oversampling factor too small";
1134 if(ths->N[j] - 1 <= ths->m)
1135 return "Polynomial degree N is smaller than cut-off m";
1137 if(ths->N[j]%2 == 1)
1138 return "polynomial degree N has to be even";
1143 void X(finalize)(
X(plan) *ths)
1153 #pragma omp critical (nfft_omp_critical_fftw_plan)
1155 FFTW(destroy_plan)(ths->my_fftw_r2r_plan);
1165 Y(free)(ths->psi_index_g);
1166 Y(free)(ths->psi_index_f);
1181 for (t = 0; t < ths->d; t++)
1182 Y(free)(ths->c_phi_inv[t]);
1183 Y(free)(ths->c_phi_inv);
1190 Y(free)(ths->f_hat);
1195 WINDOW_HELP_FINALIZE;
1199 Y(free)(ths->sigma);
1201 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.