46 #define X(name) NFCT(name) 49 static inline INT intprod(
const INT *vec,
const INT a,
const INT d)
54 for (t = 0; t < d; t++)
61 #define BASE(x) COS(x) 64 #define FOURIER_TRAFO FFTW_REDFT00 65 #define FFTW_DEFAULT_FLAGS FFTW_ESTIMATE | FFTW_DESTROY_INPUT 67 #define NODE(p,r) (ths->x[(p) * ths->d + (r)]) 69 #define MACRO_with_FG_PSI fg_psi[t][lj[t]] 70 #define MACRO_with_PRE_PSI ths->psi[(j * ths->d + t) * (2 * ths->m + 2) + lj[t]] 71 #define MACRO_without_PRE_PSI PHI((2 * NN(ths->n[t])), ((ths->x[(j) * ths->d + t]) \ 72 - ((R)(lj[t] + u[t])) / (K(2.0) * ((R)NN(ths->n[t])))), t) 73 #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) 90 void X(trafo_direct)(
const X(plan) *ths)
92 R *f_hat = (R*)ths->f_hat, *f = (R*)ths->f;
94 memset(f, 0, (
size_t)(ths->M_total) *
sizeof(R));
101 #pragma omp parallel for default(shared) private(j) 103 for (j = 0; j < ths->M_total; j++)
106 for (k_L = 0; k_L < ths->N_total; k_L++)
108 R omega = K2PI * ((R)(k_L + OFFSET)) * ths->x[j];
109 f[j] += f_hat[k_L] * BASE(omega);
118 #pragma omp parallel for default(shared) private(j) 120 for (j = 0; j < ths->M_total; j++)
122 R x[ths->d], omega, Omega[ths->d + 1];
123 INT t, t2, k_L, k[ths->d];
125 for (t = 0; t < ths->d; t++)
128 x[t] = K2PI * ths->x[j * ths->d + t];
129 Omega[t+1] = BASE(((R)(k[t])) * x[t]) * Omega[t];
131 omega = Omega[ths->d];
133 for (k_L = 0; k_L < ths->N_total; k_L++)
135 f[j] += f_hat[k_L] * omega;
137 for (t = ths->d - 1; (t >= 1) && (k[t] == (ths->N[t] - 1)); t--)
142 for (t2 = t; t2 < ths->d; t2++)
143 Omega[t2+1] = BASE(((R)(k[t2])) * x[t2]) * Omega[t2];
145 omega = Omega[ths->d];
152 void X(adjoint_direct)(
const X(plan) *ths)
154 R *f_hat = (R*)ths->f_hat, *f = (R*)ths->f;
156 memset(f_hat, 0, (
size_t)(ths->N_total) *
sizeof(R));
163 #pragma omp parallel for default(shared) private(k_L) 164 for (k_L = 0; k_L < ths->N_total; k_L++)
167 for (j = 0; j < ths->M_total; j++)
169 R omega = K2PI * ((R)(k_L + OFFSET)) * ths->x[j];
170 f_hat[k_L] += f[j] * BASE(omega);
175 for (j = 0; j < ths->M_total; j++)
178 for (k_L = 0; k_L < ths->N_total; k_L++)
180 R omega = K2PI * ((R)(k_L + OFFSET)) * ths->x[j];
181 f_hat[k_L] += f[j] * BASE(omega);
191 #pragma omp parallel for default(shared) private(j, k_L) 192 for (k_L = 0; k_L < ths->N_total; k_L++)
194 INT k[ths->d], k_temp, t;
198 for (t = ths->d - 1; t >= 0; t--)
200 k[t] = k_temp % ths->N[t];
204 for (j = 0; j < ths->M_total; j++)
207 for (t = 0; t < ths->d; t++)
208 omega *= BASE(K2PI * (k[t] + OFFSET) * ths->x[j * ths->d + t]);
209 f_hat[k_L] += f[j] * omega;
213 for (j = 0; j < ths->M_total; j++)
215 R x[ths->d], omega, Omega[ths->d+1];
216 INT t, t2, k[ths->d];
218 for (t = 0; t < ths->d; t++)
221 x[t] = K2PI * ths->x[j * ths->d + t];
222 Omega[t+1] = BASE(((R)(k[t])) * x[t]) * Omega[t];
224 omega = Omega[ths->d];
225 for (k_L = 0; k_L < ths->N_total; k_L++)
227 f_hat[k_L] += f[j] * omega;
229 for (t = ths->d-1; (t >= 1) && (k[t] == ths->N[t] - 1); t--)
234 for (t2 = t; t2 < ths->d; t2++)
235 Omega[t2+1] = BASE(((R)(k[t2])) * x[t2]) * Omega[t2];
237 omega = Omega[ths->d];
263 static inline void uo(
const X(plan) *ths,
const INT j, INT *up, INT *op,
266 const R xj = ths->x[j * ths->d + act_dim];
267 INT c = LRINT(xj * (2 * NN(ths->n[(act_dim)])));
269 (*up) = c - (ths->m);
270 (*op) = c + 1 + (ths->m);
273 #define MACRO_D_compute_A \ 275 g_hat[kg_plain[ths->d]] = f_hat[k_L] * c_phi_inv_k[ths->d]; \ 278 #define MACRO_D_compute_T \ 280 f_hat[k_L] = g_hat[kg_plain[ths->d]] * c_phi_inv_k[ths->d]; \ 283 #define MACRO_D_init_result_A memset(g_hat, 0, (size_t)(ths->n_total) * sizeof(R)); 285 #define MACRO_D_init_result_T memset(f_hat, 0, (size_t)(ths->N_total) * sizeof(R)); 287 #define MACRO_with_PRE_PHI_HUT ths->c_phi_inv[t][kg[t]] 289 #define MACRO_compute_PHI_HUT_INV (K(1.0) / (PHI_HUT((2 * NN(ths->n[t])), kg[t] + OFFSET, t))) 291 #define MACRO_init_k_ks \ 293 for (t = 0; t < ths->d; t++) \ 300 #define MACRO_update_c_phi_inv_k(what_kind, which_phi) \ 302 for (t = i; t < ths->d; t++) \ 304 MACRO_update_c_phi_inv_k_ ## what_kind(which_phi); \ 305 kg_plain[t+1] = kg_plain[t] * ths->n[t] + kg[t]; \ 309 #define MACRO_update_c_phi_inv_k_A(which_phi) \ 311 c_phi_inv_k[t+1] = (kg[t] == 0 ? K(1.0) : K(0.5)) * c_phi_inv_k[t] * MACRO_ ## which_phi; \ 314 #define MACRO_update_c_phi_inv_k_T(which_phi) \ 316 c_phi_inv_k[t+1] = c_phi_inv_k[t] * MACRO_ ## which_phi; \ 319 #define MACRO_count_k_ks \ 324 while ((kg[i] == ths->N[i]) && (i > 0)) \ 333 #define MACRO_D(which_one) \ 334 static inline void D_ ## which_one (X(plan) *ths) \ 337 R c_phi_inv_k[ths->d+1]; \ 342 INT kg_plain[ths->d+1]; \ 344 f_hat = (R*)ths->f_hat; g_hat = (R*)ths->g_hat; \ 345 MACRO_D_init_result_ ## which_one; \ 347 c_phi_inv_k[0] = K(1.0); \ 352 if (ths->flags & PRE_PHI_HUT) \ 354 for (k_L = 0; k_L < ths->N_total; k_L++) \ 356 MACRO_update_c_phi_inv_k(which_one, with_PRE_PHI_HUT); \ 357 MACRO_D_compute_ ## which_one; \ 363 for (k_L = 0; k_L < ths->N_total; k_L++) \ 365 MACRO_update_c_phi_inv_k(which_one,compute_PHI_HUT_INV); \ 366 MACRO_D_compute_ ## which_one; \ 376 #define MACRO_B_init_result_A memset(f, 0, (size_t)(ths->M_total) * sizeof(R)); 377 #define MACRO_B_init_result_T memset(g, 0, (size_t)(ths->n_total) * sizeof(R)); 379 #define MACRO_B_PRE_FULL_PSI_compute_A \ 381 (*fj) += ths->psi[ix] * g[ths->psi_index_g[ix]]; \ 384 #define MACRO_B_PRE_FULL_PSI_compute_T \ 387 INT d = ths->psi_index_g[ix]; \ 388 for (t = ths->d - 1; t >= 0; t--) \ 390 INT m = d % ths->n[t]; \ 391 if (m != 0 && m != ths->n[t] - 1) \ 395 g[ths->psi_index_g[ix]] += factor * ths->psi[ix] * (*fj); \ 398 #define MACRO_B_compute_A \ 400 (*fj) += phi_prod[ths->d] * g[ll_plain[ths->d]]; \ 403 #define MACRO_B_compute_T \ 405 g[ll_plain[ths->d]] += phi_prod[ths->d] * (*fj); \ 408 #define MACRO_init_uo_l_lj_t \ 410 for (t2 = 0; t2 < ths->d; t2++) \ 412 uo(ths, j, &u[t2], &o[t2], t2); \ 417 (u[(t2)] % (2 * NN(ths->n[(t2)]))) + (2 * NN(ths->n[(t2)])); \ 419 lg_offset[(t2)] = u[(t2)] % (2 * NN(ths->n[(t2)])); \ 420 if (lg_offset[(t2)] > NN(ths->n[(t2)])) \ 421 lg_offset[(t2)] = -(2 * NN(ths->n[(t2)]) - lg_offset[(t2)]); \ 423 if (lg_offset[t2] <= 0) \ 425 l[t2] = -lg_offset[t2]; \ 430 l[t2] = +lg_offset[t2]; \ 441 #define FOO_T ((l[t] == 0 || l[t] == ths->n[t] - 1) ? K(1.0) : K(0.5)) 443 #define MACRO_update_phi_prod_ll_plain(which_one,which_psi) \ 445 for (t = t2; t < ths->d; t++) \ 447 phi_prod[t+1] = (FOO_ ## which_one) * phi_prod[t] * (MACRO_ ## which_psi); \ 448 ll_plain[t+1] = ll_plain[t] * ths->n[t] + l[t]; \ 452 #define MACRO_count_uo_l_lj_t \ 455 if ((l[(ths->d-1)] == 0) || (l[(ths->d-1)] == NN(ths->n[(ths->d-1)]))) \ 456 count_lg[(ths->d-1)] *= -1; \ 459 l[(ths->d-1)] += count_lg[(ths->d-1)]; \ 464 while ((lj[t2] == (2 * ths->m + 2)) && (t2 > 0)) \ 471 if ((l[(t2 - 1)] == 0) || (l[(t2 - 1)] == NN(ths->n[(t2 - 1)]))) \ 472 count_lg[(t2 - 1)] *= -1; \ 474 l[(t2 - 1)] += count_lg[(t2 - 1)]; \ 477 if (lg_offset[t2] <= 0) \ 479 l[t2] = -lg_offset[t2]; \ 484 l[t2] = +lg_offset[t2]; \ 492 #define MACRO_B(which_one) \ 493 static inline void B_ ## which_one (X(plan) *ths) \ 496 INT u[ths->d], o[ths->d]; \ 502 INT ll_plain[ths->d+1]; \ 503 R phi_prod[ths->d+1]; \ 507 R fg_psi[ths->d][2*ths->m+2]; \ 508 R fg_exp_l[ths->d][2*ths->m+2]; \ 510 R tmpEXP1, tmpEXP2, tmpEXP2sq, tmp1, tmp2, tmp3; \ 513 INT ip_s = ths->K/(ths->m+2); \ 514 INT lg_offset[ths->d]; \ 515 INT count_lg[ths->d]; \ 517 f = (R*)ths->f; g = (R*)ths->g; \ 519 MACRO_B_init_result_ ## which_one \ 521 if (ths->flags & PRE_FULL_PSI) \ 523 for (ix = 0, j = 0, fj = f; j < ths->M_total; j++, fj++) \ 525 for (l_L = 0; l_L < ths->psi_index_f[j]; l_L++, ix++) \ 527 MACRO_B_PRE_FULL_PSI_compute_ ## which_one; \ 533 phi_prod[0] = K(1.0); \ 536 for (t = 0, lprod = 1; t < ths->d; t++) \ 537 lprod *= (2 * ths->m + 2); \ 539 if (ths->flags & PRE_PSI) \ 541 for (j = 0, fj = f; j < ths->M_total; j++, fj++) \ 543 MACRO_init_uo_l_lj_t; \ 545 for (l_L = 0; l_L < lprod; l_L++) \ 547 MACRO_update_phi_prod_ll_plain(which_one, with_PRE_PSI); \ 549 MACRO_B_compute_ ## which_one; \ 551 MACRO_count_uo_l_lj_t; \ 557 if (ths->flags & PRE_FG_PSI) \ 559 for (t = 0; t < ths->d; t++) \ 561 tmpEXP2 = EXP(K(-1.0) / ths->b[t]); \ 562 tmpEXP2sq = tmpEXP2 * tmpEXP2; \ 565 fg_exp_l[t][0] = K(1.0); \ 567 for (lj_fg = 1; lj_fg <= (2 * ths->m + 2); lj_fg++) \ 569 tmp3 = tmp2 * tmpEXP2; \ 571 fg_exp_l[t][lj_fg] = fg_exp_l[t][lj_fg-1] * tmp3; \ 575 for (j = 0, fj = f; j < ths->M_total; j++, fj++) \ 577 MACRO_init_uo_l_lj_t; \ 579 for (t = 0; t < ths->d; t++) \ 581 fg_psi[t][0] = ths->psi[2 * (j * ths->d + t)]; \ 582 tmpEXP1 = ths->psi[2 * (j * ths->d + t) + 1]; \ 585 for (l_fg = u[t] + 1, lj_fg = 1; l_fg <= o[t]; l_fg++, lj_fg++) \ 588 fg_psi[t][lj_fg] = fg_psi[t][0] * tmp1 * fg_exp_l[t][lj_fg]; \ 592 for (l_L= 0; l_L < lprod; l_L++) \ 594 MACRO_update_phi_prod_ll_plain(which_one, with_FG_PSI); \ 596 MACRO_B_compute_ ## which_one; \ 598 MACRO_count_uo_l_lj_t; \ 604 if (ths->flags & FG_PSI) \ 606 for (t = 0; t < ths->d; t++) \ 608 tmpEXP2 = EXP(K(-1.0) / ths->b[t]); \ 609 tmpEXP2sq = tmpEXP2 * tmpEXP2; \ 612 fg_exp_l[t][0] = K(1.0); \ 613 for (lj_fg = 1; lj_fg <= (2 * ths->m + 2); lj_fg++) \ 615 tmp3 = tmp2 * tmpEXP2; \ 617 fg_exp_l[t][lj_fg] = fg_exp_l[t][lj_fg-1] * tmp3; \ 621 for (j = 0, fj = f; j < ths->M_total; j++, fj++) \ 623 MACRO_init_uo_l_lj_t; \ 625 for (t = 0; t < ths->d; t++) \ 627 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)));\ 629 tmpEXP1 = EXP(K(2.0) * ((2 * NN(ths->n[t])) * ths->x[j * ths->d + t] - u[t]) / ths->b[t]); \ 631 for (l_fg = u[t] + 1, lj_fg = 1; l_fg <= o[t]; l_fg++, lj_fg++) \ 634 fg_psi[t][lj_fg] = fg_psi[t][0] * tmp1 * fg_exp_l[t][lj_fg]; \ 638 for (l_L = 0; l_L < lprod; l_L++) \ 640 MACRO_update_phi_prod_ll_plain(which_one, with_FG_PSI); \ 642 MACRO_B_compute_ ## which_one; \ 644 MACRO_count_uo_l_lj_t; \ 650 if (ths->flags & PRE_LIN_PSI) \ 652 for (j = 0, fj = f; j < ths->M_total; j++, fj++) \ 654 MACRO_init_uo_l_lj_t; \ 656 for (t = 0; t < ths->d; t++) \ 658 y[t] = (((2 * NN(ths->n[t])) * ths->x[j * ths->d + t] - (R)u[t]) \ 659 * ((R)ths->K))/(ths->m + 2); \ 660 ip_u = LRINT(FLOOR(y[t])); \ 662 for (l_fg = u[t], lj_fg = 0; l_fg <= o[t]; l_fg++, lj_fg++) \ 664 fg_psi[t][lj_fg] = ths->psi[(ths->K+1)*t + ABS(ip_u-lj_fg*ip_s)] \ 665 * (1-ip_w) + ths->psi[(ths->K+1)*t + ABS(ip_u-lj_fg*ip_s+1)] \ 670 for (l_L = 0; l_L < lprod; l_L++) \ 672 MACRO_update_phi_prod_ll_plain(which_one, with_FG_PSI); \ 674 MACRO_B_compute_ ## which_one; \ 676 MACRO_count_uo_l_lj_t; \ 683 for (j = 0, fj = &f[0]; j < ths->M_total; j++, fj += 1) \ 685 MACRO_init_uo_l_lj_t; \ 687 for (l_L = 0; l_L < lprod; l_L++) \ 689 MACRO_update_phi_prod_ll_plain(which_one, without_PRE_PSI); \ 691 MACRO_B_compute_ ## which_one; \ 693 MACRO_count_uo_l_lj_t; \ 704 void X(trafo)(
X(plan) *ths)
711 ths->g_hat = ths->g1;
724 FFTW(execute)(ths->my_fftw_r2r_plan);
745 void X(adjoint)(
X(plan) *ths)
752 ths->g_hat = ths->g2;
768 FFTW(execute)(ths->my_fftw_r2r_plan);
782 static inline
void precompute_phi_hut(
X(plan) *ths)
787 ths->c_phi_inv = (R**) Y(malloc)((size_t)(ths->d) *
sizeof(R*));
789 for (t = 0; t < ths->d; t++)
791 ths->c_phi_inv[t] = (R*)Y(malloc)((size_t)(ths->N[t] - OFFSET) *
sizeof(R));
793 for (ks[t] = 0; ks[t] < ths->N[t] - OFFSET; ks[t]++)
795 ths->c_phi_inv[t][ks[t]] = (K(1.0) / (PHI_HUT((2 * NN(ths->n[t])), ks[t] + OFFSET, t)));
805 void X(precompute_lin_psi)(
X(plan) *ths)
811 for (t = 0; t < ths->d; t++)
813 step = ((R)(ths->m+2)) / (((R)ths->K) * (2 * NN(ths->n[t])));
815 for (j = 0; j <= ths->K; j++)
817 ths->psi[(ths->K + 1) * t + j] = PHI((2 * NN(ths->n[t])), (j * step), t);
822 void X(precompute_fg_psi)(
X(plan) *ths)
829 for (t = 0; t < ths->d; t++)
833 for (j = 0; j < ths->M_total; j++)
835 uo(ths, j, &u, &o, t);
837 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)));
838 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]);
844 void X(precompute_psi)(
X(plan) *ths)
852 for (t = 0; t < ths->d; t++)
856 for (j = 0; j < ths->M_total; j++)
858 uo(ths, j, &u, &o, t);
860 for(lj = 0; lj < (2 * ths->m + 2); lj++)
861 ths->psi[(j * ths->d + t) * (2 * ths->m + 2) + lj] =
862 (PHI((2 * NN(ths->n[t])), ((ths->x[(j) * ths->d + (t)]) - ((R)(lj + u)) / (K(2.0) * ((R)NN(ths->n[t])))), t));
867 void X(precompute_full_psi)(
X(plan) *ths)
879 INT ll_plain[ths->d+1];
881 INT u[ths->d], o[ths->d];
882 INT count_lg[ths->d];
883 INT lg_offset[ths->d];
885 R phi_prod[ths->d+1];
891 phi_prod[0] = K(1.0);
894 for (t = 0, lprod = 1; t < ths->d; t++)
895 lprod *= 2 * ths->m + 2;
897 for (j = 0, ix = 0, ix_old = 0; j < ths->M_total; j++)
899 MACRO_init_uo_l_lj_t;
901 for (l_L = 0; l_L < lprod; l_L++, ix++)
903 MACRO_update_phi_prod_ll_plain(A, without_PRE_PSI);
905 ths->psi_index_g[ix] = ll_plain[ths->d];
906 ths->psi[ix] = phi_prod[ths->d];
908 MACRO_count_uo_l_lj_t;
911 ths->psi_index_f[j] = ix - ix_old;
917 void X(precompute_one_psi)(
X(plan) *ths)
919 if(ths->flags & PRE_PSI)
920 X(precompute_psi)(ths);
921 if(ths->flags & PRE_FULL_PSI)
922 X(precompute_full_psi)(ths);
923 if(ths->flags & PRE_FG_PSI)
924 X(precompute_fg_psi)(ths);
925 if(ths->flags & PRE_LIN_PSI)
926 X(precompute_lin_psi)(ths);
929 static inline void init_help(
X(plan) *ths)
934 if (ths->flags & NFFT_OMP_BLOCKWISE_ADJOINT)
935 ths->flags |= NFFT_SORT_NODES;
937 ths->N_total = intprod(ths->N, OFFSET, ths->d);
938 ths->n_total = intprod(ths->n, 0, ths->d);
940 ths->sigma = (R*)Y(malloc)((size_t)(ths->d) *
sizeof(R));
942 for (t = 0; t < ths->d; t++)
943 ths->sigma[t] = ((R)NN(ths->n[t])) / ths->N[t];
946 ths->r2r_kind = (FFTW(r2r_kind)*)Y(malloc)((size_t)(ths->d) *
sizeof (FFTW(r2r_kind)));
947 for (t = 0; t < ths->d; t++)
948 ths->r2r_kind[t] = FOURIER_TRAFO;
952 if (ths->flags & MALLOC_X)
953 ths->x = (R*)Y(malloc)((size_t)(ths->d * ths->M_total) *
sizeof(R));
955 if (ths->flags & MALLOC_F_HAT)
956 ths->f_hat = (R*)Y(malloc)((size_t)(ths->N_total) *
sizeof(R));
958 if (ths->flags & MALLOC_F)
959 ths->f = (R*)Y(malloc)((size_t)(ths->M_total) *
sizeof(R));
961 if (ths->flags & PRE_PHI_HUT)
962 precompute_phi_hut(ths);
964 if(ths->flags & PRE_LIN_PSI)
966 ths->K = (1U<< 10) * (ths->m+2);
967 ths->psi = (R*) Y(malloc)((size_t)((ths->K + 1) * ths->d) *
sizeof(R));
970 if(ths->flags & PRE_FG_PSI)
971 ths->psi = (R*) Y(malloc)((size_t)(ths->M_total * ths->d * 2) *
sizeof(R));
973 if (ths->flags & PRE_PSI)
974 ths->psi = (R*) Y(malloc)((size_t)(ths->M_total * ths->d * (2 * ths->m + 2 )) *
sizeof(R));
976 if(ths->flags & PRE_FULL_PSI)
978 for (t = 0, lprod = 1; t < ths->d; t++)
979 lprod *= 2 * ths->m + 2;
981 ths->psi = (R*) Y(malloc)((size_t)(ths->M_total * lprod) *
sizeof(R));
983 ths->psi_index_f = (INT*) Y(malloc)((size_t)(ths->M_total) *
sizeof(INT));
984 ths->psi_index_g = (INT*) Y(malloc)((size_t)(ths->M_total * lprod) *
sizeof(INT));
987 if (ths->flags & FFTW_INIT)
989 ths->g1 = (R*)Y(malloc)((size_t)(ths->n_total) *
sizeof(R));
991 if (ths->flags & FFT_OUT_OF_PLACE)
992 ths->g2 = (R*) Y(malloc)((size_t)(ths->n_total) *
sizeof(R));
997 int *_n = Y(malloc)((size_t)(ths->d) *
sizeof(int));
999 for (t = 0; t < ths->d; t++)
1000 _n[t] = (
int)(ths->n[t]);
1002 ths->my_fftw_r2r_plan = FFTW(plan_r2r)((int)ths->d, _n, ths->g1, ths->g2, ths->r2r_kind, ths->fftw_flags);
1012 ths->mv_trafo = (void (*) (
void* ))
X(trafo);
1013 ths->mv_adjoint = (void (*) (
void* ))
X(adjoint);
1016 void X(init)(
X(plan) *ths,
int d,
int *N,
int M_total)
1022 ths->N = (INT*) Y(malloc)((size_t)(d) *
sizeof(INT));
1024 for (t = 0; t < d; t++)
1025 ths->N[t] = (INT)N[t];
1027 ths->M_total = (INT)M_total;
1029 ths->n = (INT*) Y(malloc)((size_t)(d) *
sizeof(INT));
1031 for (t = 0; t < d; t++)
1032 ths->n[t] = 2 * (Y(next_power_of_2)(ths->N[t]) - 1) + OFFSET;
1034 ths->m = WINDOW_HELP_ESTIMATE_m;
1043 ths->flags = PRE_PHI_HUT | PRE_PSI | MALLOC_X| MALLOC_F_HAT | MALLOC_F |
1044 FFTW_INIT | FFT_OUT_OF_PLACE | NFFT_SORT_NODES;
1048 ths->flags = PRE_PHI_HUT | PRE_PSI | MALLOC_X| MALLOC_F_HAT | MALLOC_F |
1049 FFTW_INIT | FFT_OUT_OF_PLACE;
1051 ths->fftw_flags = FFTW_ESTIMATE | FFTW_DESTROY_INPUT;
1056 void X(init_guru)(
X(plan) *ths,
int d,
int *N,
int M_total,
int *n,
int m,
1057 unsigned flags,
unsigned fftw_flags)
1062 ths->M_total = (INT)M_total;
1063 ths->N = (INT*)Y(malloc)((size_t)(ths->d) *
sizeof(INT));
1065 for (t = 0; t < d; t++)
1066 ths->N[t] = (INT)N[t];
1068 ths->n = (INT*)Y(malloc)((size_t)(ths->d) *
sizeof(INT));
1070 for (t = 0; t < d; t++)
1071 ths->n[t] = (INT)n[t];
1076 ths->fftw_flags = fftw_flags;
1081 void X(init_1d)(
X(plan) *ths,
int N1,
int M_total)
1087 X(init)(ths, 1, N, M_total);
1090 void X(init_2d)(
X(plan) *ths,
int N1,
int N2,
int M_total)
1097 X(init)(ths, 2, N, M_total);
1100 void X(init_3d)(
X(plan) *ths,
int N1,
int N2,
int N3,
int M_total)
1108 X(init)(ths, 3, N, M_total);
1111 const char*
X(check)(
X(plan) *ths)
1116 return "Member f not initialized.";
1119 return "Member x not initialized.";
1122 return "Member f_hat not initialized.";
1124 for (j = 0; j < ths->M_total * ths->d; j++)
1126 if ((ths->x[j] < K(0.0)) || (ths->x[j] >= K(0.5)))
1128 return "ths->x out of range [0.0,0.5)";
1132 for (j = 0; j < ths->d; j++)
1134 if (ths->sigma[j] <= 1)
1135 return "Oversampling factor too small";
1137 if(ths->N[j] - 1 <= ths->m)
1138 return "Polynomial degree N is smaller than cut-off m";
1140 if(ths->N[j]%2 == 1)
1141 return "polynomial degree N has to be even";
1146 void X(finalize)(
X(plan) *ths)
1153 if (ths->flags & FFTW_INIT)
1156 #pragma omp critical (nfft_omp_critical_fftw_plan) 1158 FFTW(destroy_plan)(ths->my_fftw_r2r_plan);
1160 if (ths->flags & FFT_OUT_OF_PLACE)
1166 if(ths->flags & PRE_FULL_PSI)
1168 Y(free)(ths->psi_index_g);
1169 Y(free)(ths->psi_index_f);
1173 if (ths->flags & PRE_PSI)
1176 if(ths->flags & PRE_FG_PSI)
1179 if(ths->flags & PRE_LIN_PSI)
1182 if (ths->flags & PRE_PHI_HUT)
1184 for (t = 0; t < ths->d; t++)
1185 Y(free)(ths->c_phi_inv[t]);
1186 Y(free)(ths->c_phi_inv);
1189 if (ths->flags & MALLOC_F)
1192 if(ths->flags & MALLOC_F_HAT)
1193 Y(free)(ths->f_hat);
1195 if (ths->flags & MALLOC_X)
1198 WINDOW_HELP_FINALIZE;
1202 Y(free)(ths->sigma);
1204 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.