00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #include <stdio.h>
00025 #include <math.h>
00026 #include <string.h>
00027 #include <stdlib.h>
00028
00029 #include "nfft3util.h"
00030 #include "nfft3.h"
00031 #include "infft.h"
00032
00033 #undef X
00034 #if defined(NFFT_SINGLE)
00035 #define X(name) CONCAT(nfctf_, name)
00036 #elif defined(NFFT_LDOUBLE)
00037 #define X(name) CONCAT(nfctl_, name)
00038 #else
00039 #define X(name) CONCAT(nfct_, name)
00040 #endif
00041
00042 static inline int fftw_2N(int n)
00043 {
00044 return 2 * (n - 1);
00045 }
00046
00047 static inline int fftw_2N_rev(int n)
00048 {
00049 return (LRINT(K(0.5) * n) + 1);
00050 }
00051
00052 static inline int prod_int(int *vec, int d)
00053 {
00054 int t, prod = 1;
00055
00056 for (t = 0; t < d; t++)
00057 prod *= vec[t];
00058
00059 return prod;
00060 }
00061
00062
00063 #define NFCT_DEFAULT_FLAGS PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | \
00064 MALLOC_F | FFTW_INIT | FFT_OUT_OF_PLACE
00065
00066 #define FFTW_DEFAULT_FLAGS FFTW_ESTIMATE | FFTW_DESTROY_INPUT
00067
00068 #define NFCT_SUMMANDS (2 * ths->m + 2)
00069 #define NODE(p,r) (ths->x[(p) * ths->d + (r)])
00070
00071 #define MACRO_ndct_init_result_trafo \
00072 memset(f, 0, ths->M_total * sizeof(R));
00073 #define MACRO_ndct_init_result_adjoint \
00074 memset(f_hat, 0, ths->N_total * sizeof(R));
00075
00076 #define MACRO_nfct_D_init_result_A \
00077 memset(g_hat, 0, prod_int(ths->n, ths->d) * sizeof(R));
00078 #define MACRO_nfct_D_init_result_T \
00079 memset(f_hat, 0, ths->N_total * sizeof(R));
00080
00081 #define MACRO_nfct_B_init_result_A \
00082 memset(f, 0, ths->M_total * sizeof(R));
00083 #define MACRO_nfct_B_init_result_T \
00084 memset(g, 0, prod_int(ths->n, ths->d) * sizeof(R));
00085
00086 #define NFCT_PRE_WINFUN(d) ths->N[d] = 2 * ths->N[d]; \
00087 ths->n[d] = fftw_2N(ths->n[d]);
00088
00089 #define NFCT_POST_WINFUN(d) ths->N[d] = LRINT(K(0.5) * ths->N[d]); \
00090 ths->n[d] = fftw_2N_rev(ths->n[d]);
00091
00092 #define NFCT_WINDOW_HELP_INIT WINDOW_HELP_INIT
00093
00094 R X(phi_hut)(X(plan) *ths, int k, int d)
00095 {
00096 NFCT_PRE_WINFUN(d);
00097 R phi_hut_tmp = PHI_HUT(k, d);
00098 NFCT_POST_WINFUN(d);
00099
00100 return phi_hut_tmp;
00101 }
00102
00103 R X(phi)(X(plan) *ths, R x, int d)
00104 {
00105 NFCT_PRE_WINFUN(d);
00106 R phi_tmp = PHI(x, d);
00107 NFCT_POST_WINFUN(d);
00108
00109 return phi_tmp;
00110 }
00111
00112 #define MACRO_with_cos_vec cos_vec[t][ka[t]]
00113 #define MACRO_without_cos_vec COS(K(2.0) * KPI * ka[t] * NODE(j,t))
00114
00115 #define MACRO_with_PRE_PHI_HUT ths->c_phi_inv[t][kg[t]];
00116 #define MACRO_compute_PHI_HUT_INV (K(1.0) / (X(phi_hut)(ths, kg[t], t)))
00117
00118 #define MACRO_with_PRE_PSI ths->psi[(j * ths->d + t) * NFCT_SUMMANDS + lc[t]]
00119 #define MACRO_compute_PSI X(phi)(ths, NODE(j,t) - ((R)(lc[t] + lb[t])) / (K(2.0)*((R)(ths->n[t])-K(1.0))), t)
00120
00136 #define MACRO_ndct_malloc__cos_vec \
00137 R **cos_vec; \
00138 cos_vec = (R**)Y(malloc)(ths->d * sizeof(R*)); \
00139 for (t = 0; t < ths->d; t++) \
00140 cos_vec[t] = (R*)Y(malloc)(ths->N[t] * sizeof(R));
00141
00142 #define MACRO_ndct_free__cos_vec \
00143 { \
00144 \
00145 for (t = 0; t < ths->d; t++) \
00146 Y(free)(cos_vec[t]); \
00147 Y(free)(cos_vec); \
00148 }
00149
00150 #define MACRO_ndct_init__cos_vec \
00151 { \
00152 for(t = 0; t < ths->d; t++) \
00153 { \
00154 cos_vec[t][0] = K(1.0); \
00155 cos_vec[t][1] = COS(K(2.0) * KPI * NODE(j,t)); \
00156 for (k = 2; k < ths->N[t]; k++) \
00157 cos_vec[t][k] = K(2.0) * cos_vec[t][1] * cos_vec[t][k-1] - \
00158 cos_vec[t][k-2]; \
00159 } \
00160 }
00161
00162 #define MACRO_ndct_init__k__cos_k(which_one) \
00163 { \
00164 cos_k[0] = K(1.0); \
00165 for (t = 0; t < ths->d; t++) \
00166 ka[t] = 0; \
00167 \
00168 for (t = 0; t < ths->d; t++) \
00169 { \
00170 cos_k[t+1] = cos_k[t] * MACRO_ ##which_one; \
00171 } \
00172 }
00173
00174 #define MACRO_ndct_count__k__cos_k(which_one) \
00175 { \
00176 ka[ths->d-1]++; \
00177 i = ths->d - 1; \
00178 while ((ka[i] == ths->N[i]) && (i > 0)) \
00179 { \
00180 ka[i - 1]++; \
00181 ka[i] = 0; \
00182 i--; \
00183 } \
00184 for (t = i; t < ths->d; t++) \
00185 cos_k[t+1] = cos_k[t] * MACRO_ ##which_one; \
00186 }
00187
00188 #define MACRO_ndct_compute__trafo \
00189 { \
00190 f[j] += f_hat[k] * cos_k[ths->d]; \
00191 }
00192
00193 #define MACRO_ndct_compute__adjoint \
00194 { \
00195 f_hat[k] += f[j] * cos_k[ths->d]; \
00196 }
00197
00198
00199 #define MACRO_ndct(which_one) \
00200 void X(which_one ## _direct) (X(plan) *ths) \
00201 { \
00202 int j, k, t, i; \
00203 int ka[ths->d]; \
00204 R cos_k[ths->d+1]; \
00205 R *f = ths->f; \
00206 R *f_hat = ths->f_hat; \
00207 \
00208 MACRO_ndct_init_result_ ## which_one; \
00209 if (ths->d == 1) \
00210 for (j = 0; j < ths->M_total; j++) \
00211 { \
00212 for (k = 0; k < ths->N[0]; k++) \
00213 { \
00214 cos_k[ths->d] = COS(K(2.0) * KPI * k * NODE(j,0)); \
00215 MACRO_ndct_compute__ ## which_one; \
00216 } \
00217 } \
00218 else \
00219 { \
00220 \
00221 MACRO_ndct_malloc__cos_vec; \
00222 \
00223 for (j = 0; j < ths->M_total; j++) \
00224 { \
00225 MACRO_ndct_init__cos_vec; \
00226 MACRO_ndct_init__k__cos_k(with_cos_vec); \
00227 \
00228 for (k = 0; k < ths->N_total; k++) \
00229 { \
00230 MACRO_ndct_compute__ ## which_one; \
00231 MACRO_ndct_count__k__cos_k(with_cos_vec); \
00232 } \
00233 } \
00234 MACRO_ndct_free__cos_vec; \
00235 } \
00236 }
00237
00238 MACRO_ndct(trafo)
00239 MACRO_ndct(adjoint)
00240
00255 #define MACRO_nfct__lower_boundary(j,act_dim) \
00256 { \
00257 lb[(act_dim)] = \
00258 LRINT(NODE((j),(act_dim)) * fftw_2N(ths->n[(act_dim)])) - ths->m; \
00259 }
00260
00261 #define MACRO_nfct_D_compute_A \
00262 { \
00263 g_hat[kg_plain[ths->d]] = f_hat[k_L] * c_phi_inv_k[ths->d]; \
00264 }
00265
00266 #define MACRO_nfct_D_compute_T \
00267 { \
00268 f_hat[k_L] = g_hat[kg_plain[ths->d]] * c_phi_inv_k[ths->d]; \
00269 }
00270
00271 #define MACRO_init__kg \
00272 { \
00273 for (t = 0; t < ths->d; t++) \
00274 kg[t] = 0; \
00275 i = 0; \
00276 }
00277
00278 #define MACRO_count__kg \
00279 { \
00280 \
00281 kg[ths->d - 1]++; \
00282 i = ths->d - 1; \
00283 while ((kg[i] == ths->N[i]) && (i > 0)) \
00284 { \
00285 kg[i - 1]++; \
00286 kg[i] = 0; \
00287 i--; \
00288 } \
00289 }
00290
00291 #define MACRO_update__phi_inv_k__kg_plain(what_kind, which_phi) \
00292 { \
00293 for (t = i; t < ths->d; t++) \
00294 { \
00295 MACRO__c_phi_inv_k__ ## what_kind(which_phi); \
00296 kg_plain[t+1] = kg_plain[t] * ths->n[t] + kg[t]; \
00297 } \
00298 }
00299
00300 #define MACRO__c_phi_inv_k__A(which_phi) \
00301 { \
00302 if (kg[t] == 0) \
00303 { \
00304 c_phi_inv_k[t+1] = c_phi_inv_k[t] * MACRO_ ## which_phi; \
00305 } \
00306 else \
00307 { \
00308 c_phi_inv_k[t+1] = K(0.5) * c_phi_inv_k[t] * MACRO_ ## which_phi; \
00309 } \
00310 }
00311
00312 #define MACRO__c_phi_inv_k__T(which_phi) \
00313 { \
00314 c_phi_inv_k[t+1] = c_phi_inv_k[t] * MACRO_ ## which_phi; \
00315 }
00316
00317 #define MACRO_nfct_D(which_one) \
00318 static inline void D_ ## which_one (X(plan) *ths) \
00319 { \
00320 int k_L; \
00321 int i, t; \
00322 int kg[ths->d]; \
00323 R c_phi_inv_k[ths->d+1]; \
00324 int kg_plain[ths->d+1]; \
00325 R *g_hat, *f_hat; \
00326 \
00327 g_hat = ths->g_hat; \
00328 f_hat = ths->f_hat; \
00329 \
00330 MACRO_nfct_D_init_result_ ## which_one \
00331 \
00332 c_phi_inv_k[0] = K(1.0); \
00333 kg_plain[0] = 0; \
00334 \
00335 MACRO_init__kg; \
00336 \
00337 if (ths->nfct_flags & PRE_PHI_HUT) \
00338 for (k_L = 0; k_L < ths->N_total; k_L++) \
00339 { \
00340 MACRO_update__phi_inv_k__kg_plain(which_one, with_PRE_PHI_HUT); \
00341 MACRO_nfct_D_compute_ ## which_one; \
00342 MACRO_count__kg; \
00343 } \
00344 else \
00345 for (k_L = 0; k_L < ths->N_total; k_L++) \
00346 { \
00347 MACRO_update__phi_inv_k__kg_plain(which_one,compute_PHI_HUT_INV); \
00348 MACRO_nfct_D_compute_ ## which_one; \
00349 MACRO_count__kg; \
00350 } \
00351 }
00352
00353 MACRO_nfct_D(A)
00354 MACRO_nfct_D(T)
00355
00359 #define MACRO_nfct_B_PRE_FULL_PSI_compute_A \
00360 { \
00361 (*fj) += ths->psi[ix] * g[ths->psi_index_g[ix]]; \
00362 }
00363
00364 #define MACRO_nfct_B_PRE_FULL_PSI_compute_T \
00365 { \
00366 g[ths->psi_index_g[ix]] += ths->psi[ix] * (*fj); \
00367 }
00368
00369 #define MACRO_nfct_B_compute_A \
00370 { \
00371 (*fj) += phi_tilde[ths->d] * g[lg_plain[ths->d]]; \
00372 }
00373
00374 #define MACRO_nfct_B_compute_T \
00375 { \
00376 g[lg_plain[ths->d]] += phi_tilde[ths->d] * (*fj); \
00377 }
00378
00379 #define MACRO_compute_lg_offset__count_lg(i0) \
00380 { \
00381 \
00382 if (lb[(i0)] < 0) \
00383 lg_offset[(i0)] = \
00384 (lb[(i0)] % fftw_2N(ths->n[(i0)])) + fftw_2N(ths->n[(i0)]); \
00385 else \
00386 lg_offset[(i0)] = lb[(i0)] % (fftw_2N(ths->n[(i0)])); \
00387 if (lg_offset[(i0)] >= ths->n[(i0)]) \
00388 lg_offset[(i0)] = -(fftw_2N(ths->n[(i0)]) - lg_offset[(i0)]); \
00389 }
00390
00391 #define MACRO_set__lg__to__lg_offset \
00392 { \
00393 if (lg_offset[i] <= 0) \
00394 { \
00395 lg[i] = -lg_offset[i]; \
00396 count_lg[i] = -1; \
00397 } \
00398 else \
00399 { \
00400 lg[i] = +lg_offset[i]; \
00401 count_lg[i] = +1; \
00402 } \
00403 }
00404
00405 #define MACRO_count__lg(dim) \
00406 { \
00407 \
00408 if ((lg[(dim)] == 0) || (lg[(dim)] == ths->n[(dim)]-1)) \
00409 count_lg[(dim)] *= -1; \
00410 \
00411 lg[(dim)] += count_lg[(dim)]; \
00412 }
00413
00414 #define MACRO_init_lb_lg_lc \
00415 { \
00416 for (i = 0; i < ths->d; i++) \
00417 { \
00418 MACRO_nfct__lower_boundary(j, i); \
00419 MACRO_compute_lg_offset__count_lg(i); \
00420 MACRO_set__lg__to__lg_offset; \
00421 \
00422 lc[i] = 0; \
00423 } \
00424 i = 0; \
00425 }
00426
00427 #define MACRO_count__lg_lc \
00428 { \
00429 MACRO_count__lg(ths->d-1); \
00430 lc[ths->d - 1]++; \
00431 i = ths->d - 1; \
00432 while ((lc[i] == NFCT_SUMMANDS) && (i > 0)) \
00433 { \
00434 lc[i - 1]++; \
00435 lc[i] = 0; \
00436 \
00437 MACRO_count__lg(i - 1); \
00438 \
00439 MACRO_set__lg__to__lg_offset; \
00440 i--; \
00441 } \
00442 }
00443
00444 #define MACRO_update_phi_tilde_lg_plain(which_one, which_psi) \
00445 { \
00446 for (t = i; t < ths->d; t++) \
00447 { \
00448 MACRO__phi_tilde__ ## which_one(which_psi); \
00449 lg_plain[t+1] = lg_plain[t] * ths->n[t] + lg[t]; \
00450 } \
00451 }
00452
00453 #define MACRO__phi_tilde__A(which_psi) \
00454 { \
00455 phi_tilde[t+1] = phi_tilde[t] * MACRO_ ## which_psi; \
00456 }
00457
00458 #define MACRO__phi_tilde__T(which_psi) \
00459 { \
00460 if(lg[t] == 0 || lg[t] == ths->n[t] - 1) \
00461 { \
00462 phi_tilde[t+1] = phi_tilde[t] * MACRO_ ## which_psi; \
00463 } \
00464 else \
00465 { \
00466 phi_tilde[t+1] = K(0.5) * phi_tilde[t] * MACRO_ ## which_psi; \
00467 } \
00468 }
00469
00470 #define MACRO_nfct_B(which_one) \
00471 static inline void B_ ## which_one (nfct_plan *ths) \
00472 { \
00473 int lb[ths->d]; \
00474 int j, t, i; \
00475 int lprod, l_L, ix; \
00476 int lc[ths->d]; \
00477 int lg[ths->d]; \
00478 int lg_offset[ths->d]; \
00479 int count_lg[ths->d]; \
00480 int lg_plain[ths->d+1]; \
00481 R *f, *g; \
00482 R phi_tilde[ths->d+1]; \
00483 R *fj; \
00484 \
00485 f = ths->f; g = ths->g; \
00486 \
00487 MACRO_nfct_B_init_result_ ## which_one \
00488 \
00489 \
00490 if ((ths->nfct_flags & PRE_PSI)&&(ths->nfct_flags & PRE_FULL_PSI)) \
00491 { \
00492 for (ix = 0, j = 0, fj = &f[0]; j < ths->M_total; j++, fj += 1) \
00493 for (l_L = 0; l_L < ths->psi_index_f[j]; l_L++, ix++) \
00494 { \
00495 MACRO_nfct_B_PRE_FULL_PSI_compute_ ## which_one; \
00496 } \
00497 } \
00498 else \
00499 { \
00500 phi_tilde[0] = K(1.0); \
00501 lg_plain[0] = 0; \
00502 \
00503 for (t = 0, lprod = 1; t < ths->d; t++) \
00504 lprod *= NFCT_SUMMANDS; \
00505 \
00506 \
00507 if (ths->nfct_flags & PRE_PSI) \
00508 for (j = 0, fj = &f[0]; j < ths->M_total; j++, fj += 1) \
00509 { \
00510 MACRO_init_lb_lg_lc; \
00511 for (l_L = 0; l_L < lprod; l_L++) \
00512 { \
00513 MACRO_update_phi_tilde_lg_plain(which_one, with_PRE_PSI); \
00514 MACRO_nfct_B_compute_ ## which_one; \
00515 MACRO_count__lg_lc; \
00516 } \
00517 } \
00518 \
00519 \
00520 else \
00521 for (j = 0, fj = &f[0]; j < ths->M_total; j++, fj += 1) \
00522 { \
00523 MACRO_init_lb_lg_lc; \
00524 for (l_L = 0; l_L < lprod; l_L++) \
00525 { \
00526 MACRO_update_phi_tilde_lg_plain(which_one,compute_PSI); \
00527 MACRO_nfct_B_compute_ ## which_one; \
00528 MACRO_count__lg_lc; \
00529 } \
00530 } \
00531 } \
00532 }
00533
00534 MACRO_nfct_B(A)
00535 MACRO_nfct_B(T)
00536
00537
00538 #define MACRO_nfct_full_psi(which_one) \
00539 static inline void full_psi__ ## which_one(nfct_plan *ths) \
00540 { \
00541 int t, i; \
00542 int j; \
00543 int l_L; \
00544 int lc[ths->d]; \
00545 int lg_plain[ths->d+1]; \
00546 int count_lg[ths->d]; \
00547 int lg_offset[ths->d]; \
00548 int lg[ths->d]; \
00549 int lprod; \
00550 int lb[ths->d]; \
00551 \
00552 R phi_tilde[ths->d+1]; \
00553 R eps = ths->nfct_full_psi_eps; \
00554 \
00555 int *index_g, *index_f; \
00556 R *new_psi; \
00557 int ix, ix_old, size_psi; \
00558 \
00559 phi_tilde[0] = K(1.0); \
00560 lg_plain[0] = 0; \
00561 \
00562 if (ths->nfct_flags & PRE_PSI) \
00563 { \
00564 size_psi = ths->M_total; \
00565 index_f = (int*)Y(malloc)(ths->M_total * sizeof(int)); \
00566 index_g = (int*)Y(malloc)(size_psi * sizeof(int)); \
00567 new_psi = (R*)Y(malloc)(size_psi * sizeof(R)); \
00568 \
00569 for (t = 0,lprod = 1; t < ths->d; t++) \
00570 { \
00571 lprod *= NFCT_SUMMANDS; \
00572 eps *= nfct_phi(ths, K(0.0), t); \
00573 } \
00574 \
00575 for (ix = 0, ix_old = 0, j = 0; j < ths->M_total; j++) \
00576 { \
00577 MACRO_init_lb_lg_lc; \
00578 \
00579 for (l_L = 0; l_L < lprod; l_L++) \
00580 { \
00581 MACRO_update_phi_tilde_lg_plain(which_one, with_PRE_PSI); \
00582 \
00583 if (phi_tilde[ths->d] > eps) \
00584 { \
00585 index_g[ix] = lg_plain[ths->d]; \
00586 new_psi[ix] = phi_tilde[ths->d]; \
00587 \
00588 ix++; \
00589 if (ix == size_psi) \
00590 { \
00591 size_psi += ths->M_total; \
00592 index_g = (int*)realloc(index_g, size_psi * sizeof(int)); \
00593 new_psi = (R*)realloc(new_psi, size_psi * sizeof(R)); \
00594 } \
00595 } \
00596 \
00597 MACRO_count__lg_lc; \
00598 \
00599 } \
00600 \
00601 index_f[j] = ix - ix_old; \
00602 ix_old = ix; \
00603 \
00604 } \
00605 \
00606 Y(free)(ths->psi); \
00607 size_psi = ix; \
00608 ths->size_psi = size_psi; \
00609 \
00610 index_g = (int*)realloc(index_g, size_psi * sizeof(int)); \
00611 new_psi = (R*)realloc(new_psi, size_psi * sizeof(R)); \
00612 \
00613 ths->psi = new_psi; \
00614 ths->psi_index_g = index_g; \
00615 ths->psi_index_f = index_f; \
00616 \
00617 } \
00618 }
00619
00620 MACRO_nfct_full_psi(A)
00621 MACRO_nfct_full_psi(T)
00622
00623
00624
00625 void X(trafo)(X(plan) *ths)
00626 {
00627
00628 ths->g_hat = ths->g1;
00629 ths->g = ths->g2;
00630
00631
00632
00633 TIC(0)
00634 D_A(ths);
00635 TOC(0)
00636
00637
00638
00639
00640 TIC(1)
00641 Z(execute)(ths->my_fftw_r2r_plan);
00642 TOC(1)
00643
00644 if (ths->nfct_flags & PRE_FULL_PSI)
00645 full_psi__A(ths);
00646
00647
00648
00649 TIC(2)
00650 B_A(ths);
00651 TOC(2)
00652
00653 if (ths->nfct_flags & PRE_FULL_PSI)
00654 {
00655 Y(free)(ths->psi_index_g);
00656 Y(free)(ths->psi_index_f);
00657 }
00658 }
00659
00660 void X(adjoint)(X(plan) *ths)
00661 {
00662
00663 ths->g_hat = ths->g2;
00664 ths->g = ths->g1;
00665
00666 if (ths->nfct_flags & PRE_FULL_PSI)
00667 full_psi__T(ths);
00668
00669
00670
00671 TIC(2)
00672 B_T(ths);
00673 TOC(2)
00674
00675 if (ths->nfct_flags & PRE_FULL_PSI)
00676 {
00677 Y(free)(ths->psi_index_g);
00678 Y(free)(ths->psi_index_f);
00679 }
00680
00681
00682
00683
00684 TIC(1)
00685 Z(execute)(ths->my_fftw_r2r_plan);
00686 TOC(1)
00687
00688
00689
00690 TIC(0)
00691 D_T(ths);
00692 TOC(0)
00693
00694 }
00695
00696 static inline void precompute_phi_hut(X(plan) *ths)
00697 {
00698 int kg[ths->d];
00699 int t;
00700
00701 ths->c_phi_inv = (R**)Y(malloc)(ths->d * sizeof(R*));
00702
00703 for (t = 0; t < ths->d; t++)
00704 {
00705 ths->c_phi_inv[t] = (R*)Y(malloc)(ths->N[t] * sizeof(R));
00706
00707 for (kg[t] = 0; kg[t] < ths->N[t]; kg[t]++)
00708 {
00709 ths->c_phi_inv[t][kg[t]] = MACRO_compute_PHI_HUT_INV;
00710 }
00711 }
00712 }
00713
00714 void X(precompute_psi)(X(plan) *ths)
00715 {
00716 int t;
00717 int j;
00718 int lc[ths->d];
00719 int lb[ths->d];
00720
00721 for (t = 0; t < ths->d; t++)
00722 {
00723 for (j = 0; j < ths->M_total; j++)
00724 {
00725 MACRO_nfct__lower_boundary(j, t);
00726 for(lc[t] = 0; lc[t] < NFCT_SUMMANDS; lc[t]++)
00727 ths->psi[(j * ths->d + t) * NFCT_SUMMANDS + lc[t]] = MACRO_compute_PSI;
00728 }
00729 }
00730 }
00731
00732 static inline void init_help(X(plan) *ths)
00733 {
00734 int t;
00735
00736 ths->N_total = prod_int(ths->N, ths->d);
00737 ths->sigma = (R*)Y(malloc)(ths->d * sizeof(R));
00738
00739 for (t = 0; t < ths->d; t++)
00740 ths->sigma[t] = ((R)(ths->n[t] - 1)) / ths->N[t];
00741
00742
00743 ths->r2r_kind = (Z(r2r_kind)*)Y(malloc)(ths->d * sizeof (Z(r2r_kind)));
00744 for (t = 0; t < ths->d; t++)
00745 ths->r2r_kind[t] = FFTW_REDFT00;
00746
00747 NFCT_WINDOW_HELP_INIT;
00748
00749 if (ths->nfct_flags & MALLOC_X)
00750 ths->x = (R*)Y(malloc)(ths->d * ths->M_total * sizeof(R));
00751
00752 if (ths->nfct_flags & MALLOC_F_HAT)
00753 ths->f_hat = (R*)Y(malloc)(ths->N_total * sizeof(R));
00754
00755 if (ths->nfct_flags & MALLOC_F)
00756 ths->f = (R*)Y(malloc)(ths->M_total * sizeof(R));
00757
00758 if (ths->nfct_flags & PRE_PHI_HUT)
00759 precompute_phi_hut(ths);
00760
00761
00762 if (ths->nfct_flags & PRE_PSI)
00763 {
00764 ths->psi =
00765 (R*)Y(malloc)(ths->M_total * ths->d * NFCT_SUMMANDS * sizeof(R));
00766
00767
00768 ths->nfct_full_psi_eps = POW(K(10.0), K(-10.0));
00769 }
00770
00771 if (ths->nfct_flags & FFTW_INIT)
00772 {
00773 ths->g1 =
00774 (R*)Y(malloc)(prod_int(ths->n, ths->d) * sizeof(R));
00775
00776 if (ths->nfct_flags & FFT_OUT_OF_PLACE)
00777 ths->g2 =
00778 (R*) Y(malloc)(prod_int(ths->n, ths->d) * sizeof(R));
00779 else
00780 ths->g2 = ths->g1;
00781
00782 ths->my_fftw_r2r_plan =
00783 Z(plan_r2r)(ths->d, ths->n, ths->g1, ths->g2, ths->r2r_kind,
00784 ths->fftw_flags);
00785 }
00786
00787 ths->mv_trafo = (void (*) (void* ))X(trafo);
00788 ths->mv_adjoint = (void (*) (void* ))X(adjoint);
00789 }
00790
00791 void X(init)(X(plan) *ths, int d, int *N, int M_total)
00792 {
00793 int t;
00794
00795 ths->d = d;
00796 ths->M_total = M_total;
00797 ths->N = (int*) Y(malloc)(ths->d * sizeof(int));
00798
00799 for (t = 0;t < d; t++)
00800 ths->N[t] = N[t];
00801
00802 ths->n = (int*) Y(malloc)(ths->d * sizeof(int));
00803
00804 for (t = 0; t < d; t++)
00805 ths->n[t] = fftw_2N(Y(next_power_of_2)(ths->N[t]));
00806
00807
00808
00809
00810
00811 ths->nfct_flags = NFCT_DEFAULT_FLAGS;
00812 ths->fftw_flags = FFTW_DEFAULT_FLAGS;
00813
00814 init_help(ths);
00815 }
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830 void X(init_guru)(X(plan) *ths, int d, int *N, int M_total, int *n, int m,
00831 unsigned nfct_flags, unsigned fftw_flags)
00832 {
00833 int t;
00834
00835 ths->d = d;
00836 ths->M_total = M_total;
00837
00838 ths->N = (int*)Y(malloc)(ths->d * sizeof(int));
00839
00840 for (t = 0; t < d; t++)
00841 ths->N[t] = N[t];
00842
00843 ths->n = (int*)Y(malloc)(ths->d * sizeof(int));
00844
00845 for (t = 0; t < d; t++)
00846 ths->n[t] = n[t];
00847
00848 ths->m = m;
00849
00850 ths->nfct_flags = nfct_flags;
00851 ths->fftw_flags = fftw_flags;
00852
00853 init_help(ths);
00854 }
00855
00856 void X(init_1d)(X(plan) *ths, int N0, int M_total)
00857 {
00858 int N[1];
00859
00860 N[0] = N0;
00861 X(init)(ths, 1, N, M_total);
00862 }
00863
00864 void X(init_2d)(X(plan) *ths, int N0, int N1, int M_total)
00865 {
00866 int N[2];
00867
00868 N[0] = N0;
00869 N[1] = N1;
00870 X(init)(ths, 2, N, M_total);
00871 }
00872
00873 void X(init_3d)(X(plan) *ths, int N0, int N1, int N2, int M_total)
00874 {
00875 int N[3];
00876
00877 N[0] = N0;
00878 N[1] = N1;
00879 N[2] = N2;
00880 X(init)(ths, 3, N, M_total);
00881 }
00882
00883 void X(finalize)(X(plan) *ths)
00884 {
00885 int t;
00886
00887 if (ths->nfct_flags & FFTW_INIT)
00888 {
00889 Z(destroy_plan)(ths->my_fftw_r2r_plan);
00890
00891 if (ths->nfct_flags & FFT_OUT_OF_PLACE)
00892 Y(free)(ths->g2);
00893
00894 Y(free)(ths->g1);
00895 }
00896
00897
00898 if (ths->nfct_flags & PRE_PSI)
00899 {
00900 Y(free)(ths->psi);
00901 }
00902
00903 if (ths->nfct_flags & PRE_PHI_HUT)
00904 {
00905 for (t = 0; t < ths->d; t++)
00906 Y(free)(ths->c_phi_inv[t]);
00907 Y(free)(ths->c_phi_inv);
00908 }
00909
00910 if (ths->nfct_flags & MALLOC_F)
00911 Y(free)(ths->f);
00912
00913 if(ths->nfct_flags & MALLOC_F_HAT)
00914 Y(free)(ths->f_hat);
00915
00916 if (ths->nfct_flags & MALLOC_X)
00917 Y(free)(ths->x);
00918
00919 WINDOW_HELP_FINALIZE;
00920
00921 Y(free)(ths->N);
00922 Y(free)(ths->n);
00923 Y(free)(ths->sigma);
00924
00925 Y(free)(ths->r2r_kind);
00926 }