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.