46 #define X(name) NFST(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) SIN(x) 64 #define FOURIER_TRAFO FFTW_RODFT00 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] = 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] = K(0.5) * c_phi_inv_k[t] * MACRO_ ## which_phi; \ 319 #define MACRO_count_k_ks \ 324 while ((kg[i] == ths->N[i] - 1) && (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 \ 386 g[ths->psi_index_g[ix]] += ths->psi[ix] * (*fj); \ 389 #define MACRO_B_compute_A \ 391 (*fj) += phi_prod[ths->d] * g[ll_plain[ths->d]]; \ 394 #define MACRO_B_compute_T \ 396 g[ll_plain[ths->d]] += phi_prod[ths->d] * (*fj); \ 399 #define MACRO_init_uo_l_lj_t \ 401 for (t2 = 0; t2 < ths->d; t2++) \ 403 uo(ths, j, &u[t2], &o[t2], t2); \ 408 (u[(t2)] % (2 * NN(ths->n[(t2)]))) + (2 * NN(ths->n[(t2)])); \ 410 lg_offset[(t2)] = u[(t2)] % (2 * NN(ths->n[(t2)])); \ 411 if (lg_offset[(t2)] > NN(ths->n[(t2)])) \ 412 lg_offset[(t2)] = -(2 * NN(ths->n[(t2)]) - lg_offset[(t2)]); \ 414 if (lg_offset[t2] <= 0) \ 416 l[t2] = -lg_offset[t2]; \ 421 l[t2] = +lg_offset[t2]; \ 430 #define FOO_A ((R)count_lg[t]) 432 #define FOO_T ((R)count_lg[t]) 434 #define MACRO_update_phi_prod_ll_plain(which_one,which_psi) \ 436 for (t = t2; t < ths->d; t++) \ 438 if ((l[t] != 0) && (l[t] != NN(ths->n[t]))) \ 440 phi_prod[t+1] = (FOO_ ## which_one) * phi_prod[t] * (MACRO_ ## which_psi); \ 441 ll_plain[t+1] = ll_plain[t] * ths->n[t] + l[t] - 1; \ 445 phi_prod[t + 1] = K(0.0); \ 446 ll_plain[t+1] = ll_plain[t] * ths->n[t]; \ 451 #define MACRO_count_uo_l_lj_t \ 454 if ((l[(ths->d-1)] == 0) || (l[(ths->d-1)] == NN(ths->n[(ths->d-1)]))) \ 455 count_lg[(ths->d-1)] *= -1; \ 458 l[(ths->d-1)] += count_lg[(ths->d-1)]; \ 463 while ((lj[t2] == (2 * ths->m + 2)) && (t2 > 0)) \ 470 if ((l[(t2 - 1)] == 0) || (l[(t2 - 1)] == NN(ths->n[(t2 - 1)]))) \ 471 count_lg[(t2 - 1)] *= -1; \ 473 l[(t2 - 1)] += count_lg[(t2 - 1)]; \ 476 if (lg_offset[t2] <= 0) \ 478 l[t2] = -lg_offset[t2]; \ 483 l[t2] = +lg_offset[t2]; \ 491 #define MACRO_B(which_one) \ 492 static inline void B_ ## which_one (X(plan) *ths) \ 495 INT u[ths->d], o[ths->d]; \ 501 INT ll_plain[ths->d+1]; \ 502 R phi_prod[ths->d+1]; \ 506 R fg_psi[ths->d][2*ths->m+2]; \ 507 R fg_exp_l[ths->d][2*ths->m+2]; \ 509 R tmpEXP1, tmpEXP2, tmpEXP2sq, tmp1, tmp2, tmp3; \ 512 INT ip_s = ths->K/(ths->m+2); \ 513 INT lg_offset[ths->d]; \ 514 INT count_lg[ths->d]; \ 516 f = (R*)ths->f; g = (R*)ths->g; \ 518 MACRO_B_init_result_ ## which_one \ 520 if (ths->flags & PRE_FULL_PSI) \ 522 for (ix = 0, j = 0, fj = f; j < ths->M_total; j++, fj++) \ 524 for (l_L = 0; l_L < ths->psi_index_f[j]; l_L++, ix++) \ 526 MACRO_B_PRE_FULL_PSI_compute_ ## which_one; \ 532 phi_prod[0] = K(1.0); \ 535 for (t = 0, lprod = 1; t < ths->d; t++) \ 536 lprod *= (2 * ths->m + 2); \ 538 if (ths->flags & PRE_PSI) \ 540 for (j = 0, fj = f; j < ths->M_total; j++, fj++) \ 542 MACRO_init_uo_l_lj_t; \ 544 for (l_L = 0; l_L < lprod; l_L++) \ 546 MACRO_update_phi_prod_ll_plain(which_one, with_PRE_PSI); \ 548 MACRO_B_compute_ ## which_one; \ 550 MACRO_count_uo_l_lj_t; \ 556 if (ths->flags & PRE_FG_PSI) \ 558 for (t = 0; t < ths->d; t++) \ 560 tmpEXP2 = EXP(K(-1.0) / ths->b[t]); \ 561 tmpEXP2sq = tmpEXP2 * tmpEXP2; \ 564 fg_exp_l[t][0] = K(1.0); \ 566 for (lj_fg = 1; lj_fg <= (2 * ths->m + 2); lj_fg++) \ 568 tmp3 = tmp2 * tmpEXP2; \ 570 fg_exp_l[t][lj_fg] = fg_exp_l[t][lj_fg-1] * tmp3; \ 574 for (j = 0, fj = f; j < ths->M_total; j++, fj++) \ 576 MACRO_init_uo_l_lj_t; \ 578 for (t = 0; t < ths->d; t++) \ 580 fg_psi[t][0] = ths->psi[2 * (j * ths->d + t)]; \ 581 tmpEXP1 = ths->psi[2 * (j * ths->d + t) + 1]; \ 584 for (l_fg = u[t] + 1, lj_fg = 1; l_fg <= o[t]; l_fg++, lj_fg++) \ 587 fg_psi[t][lj_fg] = fg_psi[t][0] * tmp1 * fg_exp_l[t][lj_fg]; \ 591 for (l_L= 0; l_L < lprod; l_L++) \ 593 MACRO_update_phi_prod_ll_plain(which_one, with_FG_PSI); \ 595 MACRO_B_compute_ ## which_one; \ 597 MACRO_count_uo_l_lj_t; \ 603 if (ths->flags & FG_PSI) \ 605 for (t = 0; t < ths->d; t++) \ 607 tmpEXP2 = EXP(K(-1.0) / ths->b[t]); \ 608 tmpEXP2sq = tmpEXP2 * tmpEXP2; \ 611 fg_exp_l[t][0] = K(1.0); \ 612 for (lj_fg = 1; lj_fg <= (2 * ths->m + 2); lj_fg++) \ 614 tmp3 = tmp2 * tmpEXP2; \ 616 fg_exp_l[t][lj_fg] = fg_exp_l[t][lj_fg-1] * tmp3; \ 620 for (j = 0, fj = f; j < ths->M_total; j++, fj++) \ 622 MACRO_init_uo_l_lj_t; \ 624 for (t = 0; t < ths->d; t++) \ 626 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)));\ 628 tmpEXP1 = EXP(K(2.0) * ((2 * NN(ths->n[t])) * ths->x[j * ths->d + t] - u[t]) / ths->b[t]); \ 630 for (l_fg = u[t] + 1, lj_fg = 1; l_fg <= o[t]; l_fg++, lj_fg++) \ 633 fg_psi[t][lj_fg] = fg_psi[t][0] * tmp1 * fg_exp_l[t][lj_fg]; \ 637 for (l_L = 0; l_L < lprod; l_L++) \ 639 MACRO_update_phi_prod_ll_plain(which_one, with_FG_PSI); \ 641 MACRO_B_compute_ ## which_one; \ 643 MACRO_count_uo_l_lj_t; \ 649 if (ths->flags & PRE_LIN_PSI) \ 651 for (j = 0, fj = f; j < ths->M_total; j++, fj++) \ 653 MACRO_init_uo_l_lj_t; \ 655 for (t = 0; t < ths->d; t++) \ 657 y[t] = (((2 * NN(ths->n[t])) * ths->x[j * ths->d + t] - (R)u[t]) \ 658 * ((R)ths->K))/(ths->m + 2); \ 659 ip_u = LRINT(FLOOR(y[t])); \ 661 for (l_fg = u[t], lj_fg = 0; l_fg <= o[t]; l_fg++, lj_fg++) \ 663 fg_psi[t][lj_fg] = ths->psi[(ths->K+1)*t + ABS(ip_u-lj_fg*ip_s)] \ 664 * (1-ip_w) + ths->psi[(ths->K+1)*t + ABS(ip_u-lj_fg*ip_s+1)] \ 669 for (l_L = 0; l_L < lprod; l_L++) \ 671 MACRO_update_phi_prod_ll_plain(which_one, with_FG_PSI); \ 673 MACRO_B_compute_ ## which_one; \ 675 MACRO_count_uo_l_lj_t; \ 682 for (j = 0, fj = &f[0]; j < ths->M_total; j++, fj += 1) \ 684 MACRO_init_uo_l_lj_t; \ 686 for (l_L = 0; l_L < lprod; l_L++) \ 688 MACRO_update_phi_prod_ll_plain(which_one, without_PRE_PSI); \ 690 MACRO_B_compute_ ## which_one; \ 692 MACRO_count_uo_l_lj_t; \ 703 void X(trafo)(
X(plan) *ths)
710 ths->g_hat = ths->g1;
723 FFTW(execute)(ths->my_fftw_r2r_plan);
744 void X(adjoint)(
X(plan) *ths)
751 ths->g_hat = ths->g2;
767 FFTW(execute)(ths->my_fftw_r2r_plan);
781 static inline
void precompute_phi_hut(
X(plan) *ths)
786 ths->c_phi_inv = (R**) Y(malloc)((size_t)(ths->d) *
sizeof(R*));
788 for (t = 0; t < ths->d; t++)
790 ths->c_phi_inv[t] = (R*)Y(malloc)((size_t)(ths->N[t] - OFFSET) *
sizeof(R));
792 for (ks[t] = 0; ks[t] < ths->N[t] - OFFSET; ks[t]++)
794 ths->c_phi_inv[t][ks[t]] = (K(1.0) / (PHI_HUT((2 * NN(ths->n[t])), ks[t] + OFFSET, t)));
804 void X(precompute_lin_psi)(
X(plan) *ths)
810 for (t = 0; t < ths->d; t++)
812 step = ((R)(ths->m+2)) / (((R)ths->K) * (2 * NN(ths->n[t])));
814 for (j = 0; j <= ths->K; j++)
816 ths->psi[(ths->K + 1) * t + j] = PHI((2 * NN(ths->n[t])), (j * step), t);
821 void X(precompute_fg_psi)(
X(plan) *ths)
828 for (t = 0; t < ths->d; t++)
832 for (j = 0; j < ths->M_total; j++)
834 uo(ths, j, &u, &o, t);
836 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)));
837 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]);
843 void X(precompute_psi)(
X(plan) *ths)
851 for (t = 0; t < ths->d; t++)
855 for (j = 0; j < ths->M_total; j++)
857 uo(ths, j, &u, &o, t);
859 for(lj = 0; lj < (2 * ths->m + 2); lj++)
860 ths->psi[(j * ths->d + t) * (2 * ths->m + 2) + lj] =
861 (PHI((2 * NN(ths->n[t])), ((ths->x[(j) * ths->d + (t)]) - ((R)(lj + u)) / (K(2.0) * ((R)NN(ths->n[t])))), t));
866 void X(precompute_full_psi)(
X(plan) *ths)
878 INT ll_plain[ths->d+1];
880 INT u[ths->d], o[ths->d];
881 INT count_lg[ths->d];
882 INT lg_offset[ths->d];
884 R phi_prod[ths->d+1];
890 phi_prod[0] = K(1.0);
893 for (t = 0, lprod = 1; t < ths->d; t++)
894 lprod *= 2 * ths->m + 2;
896 for (j = 0, ix = 0, ix_old = 0; j < ths->M_total; j++)
898 MACRO_init_uo_l_lj_t;
900 for (l_L = 0; l_L < lprod; l_L++, ix++)
902 MACRO_update_phi_prod_ll_plain(A, without_PRE_PSI);
904 ths->psi_index_g[ix] = ll_plain[ths->d];
905 ths->psi[ix] = phi_prod[ths->d];
907 MACRO_count_uo_l_lj_t;
910 ths->psi_index_f[j] = ix - ix_old;
916 void X(precompute_one_psi)(
X(plan) *ths)
918 if(ths->flags & PRE_PSI)
919 X(precompute_psi)(ths);
920 if(ths->flags & PRE_FULL_PSI)
921 X(precompute_full_psi)(ths);
922 if(ths->flags & PRE_FG_PSI)
923 X(precompute_fg_psi)(ths);
924 if(ths->flags & PRE_LIN_PSI)
925 X(precompute_lin_psi)(ths);
928 static inline void init_help(
X(plan) *ths)
933 if (ths->flags & NFFT_OMP_BLOCKWISE_ADJOINT)
934 ths->flags |= NFFT_SORT_NODES;
936 ths->N_total = intprod(ths->N, OFFSET, ths->d);
937 ths->n_total = intprod(ths->n, 0, ths->d);
939 ths->sigma = (R*)Y(malloc)((size_t)(ths->d) *
sizeof(R));
941 for (t = 0; t < ths->d; t++)
942 ths->sigma[t] = ((R)NN(ths->n[t])) / ths->N[t];
945 ths->r2r_kind = (FFTW(r2r_kind)*)Y(malloc)((size_t)(ths->d) *
sizeof (FFTW(r2r_kind)));
946 for (t = 0; t < ths->d; t++)
947 ths->r2r_kind[t] = FOURIER_TRAFO;
951 if (ths->flags & MALLOC_X)
952 ths->x = (R*)Y(malloc)((size_t)(ths->d * ths->M_total) *
sizeof(R));
954 if (ths->flags & MALLOC_F_HAT)
955 ths->f_hat = (R*)Y(malloc)((size_t)(ths->N_total) *
sizeof(R));
957 if (ths->flags & MALLOC_F)
958 ths->f = (R*)Y(malloc)((size_t)(ths->M_total) *
sizeof(R));
960 if (ths->flags & PRE_PHI_HUT)
961 precompute_phi_hut(ths);
963 if(ths->flags & PRE_LIN_PSI)
965 ths->K = (1U<< 10) * (ths->m+2);
966 ths->psi = (R*) Y(malloc)((size_t)((ths->K + 1) * ths->d) *
sizeof(R));
969 if(ths->flags & PRE_FG_PSI)
970 ths->psi = (R*) Y(malloc)((size_t)(ths->M_total * ths->d * 2) *
sizeof(R));
972 if (ths->flags & PRE_PSI)
973 ths->psi = (R*) Y(malloc)((size_t)(ths->M_total * ths->d * (2 * ths->m + 2 )) *
sizeof(R));
975 if(ths->flags & PRE_FULL_PSI)
977 for (t = 0, lprod = 1; t < ths->d; t++)
978 lprod *= 2 * ths->m + 2;
980 ths->psi = (R*) Y(malloc)((size_t)(ths->M_total * lprod) *
sizeof(R));
982 ths->psi_index_f = (INT*) Y(malloc)((size_t)(ths->M_total) *
sizeof(INT));
983 ths->psi_index_g = (INT*) Y(malloc)((size_t)(ths->M_total * lprod) *
sizeof(INT));
986 if (ths->flags & FFTW_INIT)
988 ths->g1 = (R*)Y(malloc)((size_t)(ths->n_total) *
sizeof(R));
990 if (ths->flags & FFT_OUT_OF_PLACE)
991 ths->g2 = (R*) Y(malloc)((size_t)(ths->n_total) *
sizeof(R));
996 int *_n = Y(malloc)((size_t)(ths->d) *
sizeof(int));
998 for (t = 0; t < ths->d; t++)
999 _n[t] = (
int)(ths->n[t]);
1001 ths->my_fftw_r2r_plan = FFTW(plan_r2r)((int)ths->d, _n, ths->g1, ths->g2, ths->r2r_kind, ths->fftw_flags);
1011 ths->mv_trafo = (void (*) (
void* ))
X(trafo);
1012 ths->mv_adjoint = (void (*) (
void* ))
X(adjoint);
1015 void X(init)(
X(plan) *ths,
int d,
int *N,
int M_total)
1021 ths->N = (INT*) Y(malloc)((size_t)(d) *
sizeof(INT));
1023 for (t = 0; t < d; t++)
1024 ths->N[t] = (INT)N[t];
1026 ths->M_total = (INT)M_total;
1028 ths->n = (INT*) Y(malloc)((size_t)(d) *
sizeof(INT));
1030 for (t = 0; t < d; t++)
1031 ths->n[t] = 2 * (Y(next_power_of_2)(ths->N[t]) - 1) + OFFSET;
1033 ths->m = WINDOW_HELP_ESTIMATE_m;
1042 ths->flags = PRE_PHI_HUT | PRE_PSI | MALLOC_X| MALLOC_F_HAT | MALLOC_F |
1043 FFTW_INIT | FFT_OUT_OF_PLACE | NFFT_SORT_NODES;
1047 ths->flags = PRE_PHI_HUT | PRE_PSI | MALLOC_X| MALLOC_F_HAT | MALLOC_F |
1048 FFTW_INIT | FFT_OUT_OF_PLACE;
1050 ths->fftw_flags = FFTW_ESTIMATE | FFTW_DESTROY_INPUT;
1055 void X(init_guru)(
X(plan) *ths,
int d,
int *N,
int M_total,
int *n,
int m,
1056 unsigned flags,
unsigned fftw_flags)
1061 ths->M_total = (INT)M_total;
1062 ths->N = (INT*)Y(malloc)((size_t)(ths->d) *
sizeof(INT));
1064 for (t = 0; t < d; t++)
1065 ths->N[t] = (INT)N[t];
1067 ths->n = (INT*)Y(malloc)((size_t)(ths->d) *
sizeof(INT));
1069 for (t = 0; t < d; t++)
1070 ths->n[t] = (INT)n[t];
1075 ths->fftw_flags = fftw_flags;
1080 void X(init_1d)(
X(plan) *ths,
int N1,
int M_total)
1086 X(init)(ths, 1, N, M_total);
1089 void X(init_2d)(
X(plan) *ths,
int N1,
int N2,
int M_total)
1096 X(init)(ths, 2, N, M_total);
1099 void X(init_3d)(
X(plan) *ths,
int N1,
int N2,
int N3,
int M_total)
1107 X(init)(ths, 3, N, M_total);
1110 const char*
X(check)(
X(plan) *ths)
1115 return "Member f not initialized.";
1118 return "Member x not initialized.";
1121 return "Member f_hat not initialized.";
1123 for (j = 0; j < ths->M_total * ths->d; j++)
1125 if ((ths->x[j] < K(0.0)) || (ths->x[j] >= K(0.5)))
1127 return "ths->x out of range [0.0,0.5)";
1131 for (j = 0; j < ths->d; j++)
1133 if (ths->sigma[j] <= 1)
1134 return "Oversampling factor too small";
1136 if(ths->N[j] - 1 <= ths->m)
1137 return "Polynomial degree N is smaller than cut-off m";
1139 if(ths->N[j]%2 == 1)
1140 return "polynomial degree N has to be even";
1145 void X(finalize)(
X(plan) *ths)
1152 if (ths->flags & FFTW_INIT)
1155 #pragma omp critical (nfft_omp_critical_fftw_plan) 1157 FFTW(destroy_plan)(ths->my_fftw_r2r_plan);
1159 if (ths->flags & FFT_OUT_OF_PLACE)
1165 if(ths->flags & PRE_FULL_PSI)
1167 Y(free)(ths->psi_index_g);
1168 Y(free)(ths->psi_index_f);
1172 if (ths->flags & PRE_PSI)
1175 if(ths->flags & PRE_FG_PSI)
1178 if(ths->flags & PRE_LIN_PSI)
1181 if (ths->flags & PRE_PHI_HUT)
1183 for (t = 0; t < ths->d; t++)
1184 Y(free)(ths->c_phi_inv[t]);
1185 Y(free)(ths->c_phi_inv);
1188 if (ths->flags & MALLOC_F)
1191 if(ths->flags & MALLOC_F_HAT)
1192 Y(free)(ths->f_hat);
1194 if (ths->flags & MALLOC_X)
1197 WINDOW_HELP_FINALIZE;
1201 Y(free)(ths->sigma);
1203 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.