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.