00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00027 #include <stdio.h>
00028 #include <math.h>
00029 #include <string.h>
00030 #include <stdlib.h>
00031
00032 #include "nfft3util.h"
00033 #include "nfft3.h"
00034 #include "infft.h"
00035
00039 #define NFCT_DEFAULT_FLAGS PRE_PHI_HUT|\
00040 PRE_PSI|\
00041 MALLOC_X|\
00042 MALLOC_F_HAT|\
00043 MALLOC_F|\
00044 FFTW_INIT|\
00045 FFT_OUT_OF_PLACE
00046
00047 #define FFTW_DEFAULT_FLAGS FFTW_ESTIMATE|\
00048 FFTW_DESTROY_INPUT
00049
00050 #define NFCT_SUMMANDS (2 * ths->m + 2)
00051 #define NODE(p,r) (ths->x[(p) * ths->d + (r)])
00052
00053 #define MACRO_ndct_init_result_trafo \
00054 memset(f, 0, ths->M_total * sizeof(double));
00055 #define MACRO_ndct_init_result_adjoint \
00056 memset(f_hat, 0, ths->N_total * sizeof(double));
00057
00058
00059 #define MACRO_nfct_D_init_result_A \
00060 memset(g_hat, 0, nfct_prod_int(ths->n, ths->d) * sizeof(double));
00061 #define MACRO_nfct_D_init_result_T \
00062 memset(f_hat, 0, ths->N_total * sizeof(double));
00063
00064 #define MACRO_nfct_B_init_result_A \
00065 memset(f, 0, ths->M_total * sizeof(double));
00066 #define MACRO_nfct_B_init_result_T \
00067 memset(g, 0, nfct_prod_int(ths->n, ths->d) * sizeof(double));
00068
00069
00070 #define NFCT_PRE_WINFUN(d) ths->N[d] = 2 * ths->N[d]; \
00071 ths->n[d] = nfct_fftw_2N(ths->n[d]);
00072
00073 #define NFCT_POST_WINFUN(d) ths->N[d] = LRINT(0.5 * ths->N[d]); \
00074 ths->n[d] = nfct_fftw_2N_rev(ths->n[d]);
00075
00076
00077 #define NFCT_WINDOW_HELP_INIT WINDOW_HELP_INIT
00078
00079
00080 double nfct_phi_hut(nfct_plan *ths, int k, int d)
00081 {
00082 NFCT_PRE_WINFUN(d);
00083 double phi_hut_tmp = PHI_HUT(k, d);
00084 NFCT_POST_WINFUN(d);
00085
00086 return phi_hut_tmp;
00087 }
00088
00089 double nfct_phi(nfct_plan *ths, double x, int d)
00090 {
00091 NFCT_PRE_WINFUN(d);
00092 double phi_tmp = PHI(x, d);
00093 NFCT_POST_WINFUN(d);
00094
00095 return phi_tmp;
00096 }
00097
00098 int nfct_fftw_2N(int n)
00099 {
00100 return 2 * (n - 1);
00101 }
00102
00103 int nfct_fftw_2N_rev(int n)
00104 {
00105 return (LRINT(0.5 * n) + 1);
00106 }
00107
00108
00109 #define MACRO_with_cos_vec cos_vec[t][ka[t]]
00110 #define MACRO_without_cos_vec cos(2.0 * PI * ka[t] * NODE(j,t))
00111
00112
00113 #define MACRO_with_PRE_PHI_HUT ths->c_phi_inv[t][kg[t]];
00114 #define MACRO_compute_PHI_HUT_INV (1.0 / (nfct_phi_hut(ths, kg[t], t)))
00115
00116 #define MACRO_with_PRE_PSI ths->psi[(j * ths->d + t) * NFCT_SUMMANDS + lc[t]]
00117 #define MACRO_compute_PSI nfct_phi(ths, NODE(j,t) - ((double)(lc[t] + lb[t])) / (double)nfct_fftw_2N(ths->n[t]), t)
00118
00119
00120
00121
00138 #define MACRO_ndct_malloc__cos_vec \
00139 \
00140 double **cos_vec; \
00141 cos_vec = (double**)nfft_malloc(ths->d * sizeof(double*)); \
00142 for(t = 0; t < ths->d; t++) \
00143 cos_vec[t] = (double*)nfft_malloc(ths->N[t] * sizeof(double)); \
00144
00145
00146
00147
00148 #define MACRO_ndct_free__cos_vec \
00149 { \
00150 \
00151 for(t = 0; t < ths->d; t++) \
00152 nfft_free(cos_vec[t]); \
00153 nfft_free(cos_vec); \
00154 }
00155
00156
00157
00158 #define MACRO_ndct_init__cos_vec \
00159 { \
00160 for(t = 0; t < ths->d; t++) \
00161 { \
00162 cos_vec[t][0] = 1.0; \
00163 cos_vec[t][1] = cos(2.0 * PI * NODE(j,t)); \
00164 for(k = 2; k < ths->N[t]; k++) \
00165 cos_vec[t][k] = 2.0 * cos_vec[t][1] * cos_vec[t][k-1] \
00166 - cos_vec[t][k-2]; \
00167 } \
00168 }
00169
00170
00171
00172 #define MACRO_ndct_init__k__cos_k(which_one) \
00173 { \
00174 cos_k[0] = 1.0; \
00175 for(t = 0; t < ths->d; t++) \
00176 ka[t] = 0; \
00177 \
00178 for(t = 0; t < ths->d; t++) \
00179 { \
00180 cos_k[t+1] = cos_k[t] * MACRO_ ##which_one; \
00181 } \
00182 }
00183
00184
00185
00186 #define MACRO_ndct_count__k__cos_k(which_one) \
00187 { \
00188 ka[ths->d-1]++; \
00189 i = ths->d - 1; \
00190 while((ka[i] == ths->N[i]) && (i>0)) \
00191 { \
00192 ka[i - 1]++; \
00193 ka[i] = 0; \
00194 \
00195 i--; \
00196 } \
00197 for(t = i; t < ths->d; t++) \
00198 cos_k[t+1] = cos_k[t] * MACRO_ ##which_one; \
00199 }
00200
00201
00202
00203 #define MACRO_ndct_compute__trafo \
00204 { \
00205 f[j] += f_hat[k] * cos_k[ths->d]; \
00206 }
00207
00208 #define MACRO_ndct_compute__adjoint \
00209 { \
00210 f_hat[k] += f[j] * cos_k[ths->d]; \
00211 }
00212
00213
00214 #define MACRO_ndct(which_one) \
00215 void ndct_ ## which_one (nfct_plan *ths) \
00216 { \
00217 int j, k, t, i; \
00218 int ka[ths->d]; \
00219 double cos_k[ths->d+1]; \
00220 \
00221 double *f = ths->f; \
00222 double *f_hat = ths->f_hat; \
00223 \
00224 MACRO_ndct_init_result_ ## which_one; \
00225 if(ths->d == 1) \
00226 for(j = 0; j < ths->M_total; j++) \
00227 { \
00228 for(k = 0; k < ths->N[0]; k++) \
00229 { \
00230 cos_k[ths->d] = cos(2.0 * PI * k * NODE(j,0)); \
00231 MACRO_ndct_compute__ ## which_one; \
00232 } \
00233 } \
00234 else \
00235 if(1 == 0) \
00236 \
00237 for(j = 0; j < ths->M_total; j++) \
00238 { \
00239 MACRO_ndct_init__k__cos_k(without_cos_vec); \
00240 \
00241 for(k = 0; k < ths->N_total; k++) \
00242 { \
00243 MACRO_ndct_compute__ ## which_one; \
00244 \
00245 MACRO_ndct_count__k__cos_k(without_cos_vec); \
00246 } \
00247 } \
00248 else \
00249 { \
00250 \
00251 MACRO_ndct_malloc__cos_vec; \
00252 \
00253 for(j = 0; j < ths->M_total; j++) \
00254 { \
00255 \
00256 MACRO_ndct_init__cos_vec; \
00257 \
00258 MACRO_ndct_init__k__cos_k(with_cos_vec); \
00259 \
00260 for(k = 0; k < ths->N_total; k++) \
00261 { \
00262 \
00263 MACRO_ndct_compute__ ## which_one; \
00264 \
00265 MACRO_ndct_count__k__cos_k(with_cos_vec); \
00266 } \
00267 } \
00268 MACRO_ndct_free__cos_vec; \
00269 } \
00270 }
00271
00272
00273 MACRO_ndct(trafo)
00274 MACRO_ndct(adjoint)
00275
00276
00277
00278
00295 #define MACRO_nfct__lower_boundary(j,act_dim) \
00296 { \
00297 lb[(act_dim)] = \
00298 LRINT(NODE((j),(act_dim)) * nfct_fftw_2N(ths->n[(act_dim)])) - ths->m; \
00299 }
00300
00301 #define MACRO_nfct_D_compute_A \
00302 { \
00303 g_hat[kg_plain[ths->d]] = f_hat[k_L] * c_phi_inv_k[ths->d]; \
00304 }
00305
00306 #define MACRO_nfct_D_compute_T \
00307 { \
00308 f_hat[k_L] = g_hat[kg_plain[ths->d]] * c_phi_inv_k[ths->d]; \
00309 }
00310
00311
00312
00313 #define MACRO_init__kg \
00314 { \
00315 for(t = 0; t < ths->d; t++) \
00316 kg[t] = 0; \
00317 \
00318 i = 0; \
00319 }
00320
00321
00322 #define MACRO_count__kg \
00323 { \
00324 \
00325 kg[ths->d - 1]++; \
00326 i = ths->d - 1; \
00327 while((kg[i] == ths->N[i]) && (i > 0)) \
00328 { \
00329 kg[i - 1]++; \
00330 kg[i] = 0; \
00331 \
00332 i--; \
00333 } \
00334 }
00335
00336
00337 #define MACRO_update__phi_inv_k__kg_plain(what_kind, which_phi) \
00338 { \
00339 for(t = i; t < ths->d; t++) \
00340 { \
00341 MACRO__c_phi_inv_k__ ## what_kind(which_phi); \
00342 kg_plain[t+1] = kg_plain[t] * ths->n[t] + kg[t]; \
00343 } \
00344 }
00345
00346
00347 #define MACRO__c_phi_inv_k__A(which_phi) \
00348 { \
00349 if(kg[t] == 0) \
00350 { \
00351 c_phi_inv_k[t+1] = c_phi_inv_k[t] * MACRO_ ## which_phi; \
00352 } \
00353 else \
00354 { \
00355 c_phi_inv_k[t+1] = 0.5 * c_phi_inv_k[t] * MACRO_ ## which_phi; \
00356 } \
00357 }
00358
00359
00360 #define MACRO__c_phi_inv_k__T(which_phi) \
00361 { \
00362 c_phi_inv_k[t+1] = c_phi_inv_k[t] * MACRO_ ## which_phi; \
00363 }
00364
00365
00366
00367 #define MACRO_nfct_D(which_one) \
00368 static inline void nfct_D_ ## which_one (nfct_plan *ths) \
00369 { \
00370 \
00371 int k_L; \
00372 \
00373 int i, t; \
00374 int kg[ths->d]; \
00375 double c_phi_inv_k[ths->d+1]; \
00376 int kg_plain[ths->d+1]; \
00377 \
00378 double *g_hat, *f_hat; \
00379 \
00380 g_hat = ths->g_hat; \
00381 f_hat = ths->f_hat; \
00382 \
00383 MACRO_nfct_D_init_result_ ## which_one \
00384 \
00385 c_phi_inv_k[0] = 1.0; \
00386 kg_plain[0] = 0; \
00387 \
00388 MACRO_init__kg; \
00389 \
00390 if(ths->nfct_flags & PRE_PHI_HUT) \
00391 for(k_L = 0; k_L < ths->N_total; k_L++) \
00392 { \
00393 MACRO_update__phi_inv_k__kg_plain(which_one, with_PRE_PHI_HUT); \
00394 \
00395 MACRO_nfct_D_compute_ ## which_one; \
00396 \
00397 MACRO_count__kg; \
00398 \
00399 } \
00400 else \
00401 for(k_L = 0; k_L < ths->N_total; k_L++) \
00402 { \
00403 MACRO_update__phi_inv_k__kg_plain(which_one,compute_PHI_HUT_INV); \
00404 \
00405 MACRO_nfct_D_compute_ ## which_one; \
00406 \
00407 MACRO_count__kg; \
00408 \
00409 } \
00410 }
00411
00412 MACRO_nfct_D(A)
00413 MACRO_nfct_D(T)
00414
00415
00416
00417
00418
00419
00420
00424 #define MACRO_nfct_B_PRE_FULL_PSI_compute_A \
00425 { \
00426 (*fj) += ths->psi[ix] * g[ths->psi_index_g[ix]]; \
00427 }
00428
00429 #define MACRO_nfct_B_PRE_FULL_PSI_compute_T \
00430 { \
00431 g[ths->psi_index_g[ix]] += ths->psi[ix] * (*fj); \
00432 }
00433
00434
00435
00436 #define MACRO_nfct_B_compute_A \
00437 { \
00438 (*fj) += phi_tilde[ths->d] * g[lg_plain[ths->d]]; \
00439 }
00440
00441 #define MACRO_nfct_B_compute_T \
00442 { \
00443 g[lg_plain[ths->d]] += phi_tilde[ths->d] * (*fj); \
00444 }
00445
00446
00447 #define MACRO_compute_lg_offset__count_lg(i0) \
00448 { \
00449 \
00450 if(lb[(i0)] < 0) \
00451 lg_offset[(i0)] = \
00452 (lb[(i0)] % nfct_fftw_2N(ths->n[(i0)])) + nfct_fftw_2N(ths->n[(i0)]); \
00453 else \
00454 lg_offset[(i0)] = lb[(i0)] % (nfct_fftw_2N(ths->n[(i0)])); \
00455 if(lg_offset[(i0)] >= ths->n[(i0)]) \
00456 lg_offset[(i0)] = -(nfct_fftw_2N(ths->n[(i0)]) - lg_offset[(i0)]); \
00457 }
00458
00459 #define MACRO_set__lg__to__lg_offset \
00460 { \
00461 if(lg_offset[i] <= 0) \
00462 { \
00463 lg[i] = -lg_offset[i]; \
00464 count_lg[i] = -1; \
00465 } \
00466 else \
00467 { \
00468 lg[i] = +lg_offset[i]; \
00469 count_lg[i] = +1; \
00470 } \
00471 }
00472
00473 #define MACRO_count__lg(dim) \
00474 { \
00475 \
00476 \
00477 if((lg[(dim)] == 0) || (lg[(dim)] == ths->n[(dim)]-1)) \
00478 count_lg[(dim)] *= -1; \
00479 \
00480 lg[(dim)] += count_lg[(dim)]; \
00481 }
00482
00483
00484
00485
00486 #define MACRO_init_lb_lg_lc \
00487 { \
00488 for(i = 0; i < ths->d; i++) \
00489 { \
00490 MACRO_nfct__lower_boundary(j, i); \
00491 \
00492 MACRO_compute_lg_offset__count_lg(i); \
00493 MACRO_set__lg__to__lg_offset; \
00494 \
00495 \
00496 lc[i] = 0; \
00497 } \
00498 \
00499 i = 0; \
00500 }
00501
00502
00503
00504 #define MACRO_count__lg_lc \
00505 { \
00506 MACRO_count__lg(ths->d-1); \
00507 \
00508 lc[ths->d - 1]++; \
00509 i = ths->d - 1; \
00510 while((lc[i] == NFCT_SUMMANDS) && (i > 0)) \
00511 { \
00512 lc[i - 1]++; \
00513 lc[i] = 0; \
00514 \
00515 \
00516 MACRO_count__lg(i - 1); \
00517 \
00518 MACRO_set__lg__to__lg_offset; \
00519 \
00520 i--; \
00521 } \
00522 }
00523
00524 #define MACRO_update_phi_tilde_lg_plain(which_one, which_psi) \
00525 { \
00526 for(t = i; t < ths->d; t++) \
00527 { \
00528 MACRO__phi_tilde__ ## which_one(which_psi); \
00529 lg_plain[t+1] = lg_plain[t] * ths->n[t] + lg[t]; \
00530 } \
00531 }
00532
00533 #define MACRO__phi_tilde__A(which_psi) \
00534 { \
00535 phi_tilde[t+1] = phi_tilde[t] * MACRO_ ## which_psi; \
00536 }
00537
00538 #define MACRO__phi_tilde__T(which_psi) \
00539 { \
00540 if(lg[t] == 0 || lg[t] == ths->n[t] - 1) \
00541 { \
00542 phi_tilde[t+1] = phi_tilde[t] * MACRO_ ## which_psi; \
00543 } \
00544 else \
00545 { \
00546 phi_tilde[t+1] = 0.5 * phi_tilde[t] * MACRO_ ## which_psi; \
00547 } \
00548 }
00549
00550
00551
00552 #define MACRO_nfct_B(which_one) \
00553 static inline void nfct_B_ ## which_one (nfct_plan *ths) \
00554 { \
00555 int lb[ths->d]; \
00556 int j, t, i; \
00557 int lprod, l_L, ix; \
00558 int lc[ths->d]; \
00559 int lg[ths->d]; \
00560 int lg_offset[ths->d]; \
00561 int count_lg[ths->d]; \
00562 int lg_plain[ths->d+1]; \
00563 double *f, *g; \
00564 double phi_tilde[ths->d+1]; \
00565 double *fj; \
00566 \
00567 f = ths->f; g = ths->g; \
00568 \
00569 MACRO_nfct_B_init_result_ ## which_one \
00570 \
00571 \
00572 if((ths->nfct_flags & PRE_PSI)&&(ths->nfct_flags & PRE_FULL_PSI)) \
00573 { \
00574 for(ix = 0, j = 0, fj = &f[0]; j < ths->M_total; j++, fj += 1) \
00575 for(l_L = 0; l_L < ths->psi_index_f[j]; l_L++, ix++) \
00576 { \
00577 MACRO_nfct_B_PRE_FULL_PSI_compute_ ## which_one; \
00578 } \
00579 } \
00580 else \
00581 { \
00582 phi_tilde[0] = 1; \
00583 lg_plain[0] = 0; \
00584 \
00585 for(t = 0, lprod = 1; t < ths->d; t++) \
00586 lprod *= NFCT_SUMMANDS; \
00587 \
00588 \
00589 \
00590 \
00591 if(ths->nfct_flags & PRE_PSI) \
00592 for(j = 0, fj = &f[0]; j < ths->M_total; j++, fj += 1) \
00593 { \
00594 MACRO_init_lb_lg_lc; \
00595 \
00596 for(l_L = 0; l_L < lprod; l_L++) \
00597 { \
00598 MACRO_update_phi_tilde_lg_plain(which_one, with_PRE_PSI); \
00599 \
00600 MACRO_nfct_B_compute_ ## which_one; \
00601 \
00602 MACRO_count__lg_lc; \
00603 } \
00604 } \
00605 \
00606 \
00607 \
00608 \
00609 else \
00610 for(j = 0, fj = &f[0]; j < ths->M_total; j++, fj += 1) \
00611 { \
00612 MACRO_init_lb_lg_lc; \
00613 \
00614 for(l_L = 0; l_L < lprod; l_L++) \
00615 { \
00616 MACRO_update_phi_tilde_lg_plain(which_one,compute_PSI); \
00617 \
00618 MACRO_nfct_B_compute_ ## which_one; \
00619 \
00620 MACRO_count__lg_lc; \
00621 \
00622 } \
00623 } \
00624 } \
00625 }
00626
00627 MACRO_nfct_B(A)
00628 MACRO_nfct_B(T)
00629
00630
00631
00632
00634 #define MACRO_nfct_full_psi(which_one) \
00635 static void nfct_full_psi__ ## which_one(nfct_plan *ths) \
00636 { \
00637 int t, i; \
00638 int j; \
00639 int l_L; \
00640 int lc[ths->d]; \
00641 int lg_plain[ths->d+1]; \
00642 int count_lg[ths->d]; \
00643 int lg_offset[ths->d]; \
00644 int lg[ths->d]; \
00645 int lprod; \
00646 int lb[ths->d]; \
00647 \
00648 double phi_tilde[ths->d+1]; \
00649 double eps = ths->nfct_full_psi_eps; \
00650 \
00651 int *index_g, *index_f; \
00652 double *new_psi; \
00653 int ix, ix_old, size_psi; \
00654 \
00655 phi_tilde[0] = 1.0; \
00656 lg_plain[0] = 0; \
00657 \
00658 if(ths->nfct_flags & PRE_PSI) \
00659 { \
00660 size_psi = ths->M_total; \
00661 index_f = (int*)nfft_malloc(ths->M_total * sizeof(int)); \
00662 index_g = (int*)nfft_malloc(size_psi * sizeof(int)); \
00663 new_psi = (double*)nfft_malloc(size_psi * sizeof(double)); \
00664 \
00665 for(t = 0,lprod = 1; t < ths->d; t++) \
00666 { \
00667 lprod *= NFCT_SUMMANDS; \
00668 eps *= nfct_phi(ths, 0, t); \
00669 } \
00670 \
00671 for(ix = 0, ix_old = 0, j = 0; j < ths->M_total; j++) \
00672 { \
00673 MACRO_init_lb_lg_lc; \
00674 \
00675 for(l_L = 0; l_L < lprod; l_L++) \
00676 { \
00677 MACRO_update_phi_tilde_lg_plain(which_one, with_PRE_PSI); \
00678 \
00679 if(phi_tilde[ths->d] > eps) \
00680 { \
00681 index_g[ix] = lg_plain[ths->d]; \
00682 new_psi[ix] = phi_tilde[ths->d]; \
00683 \
00684 ix++; \
00685 if(ix == size_psi) \
00686 { \
00687 size_psi += ths->M_total; \
00688 index_g = (int*)realloc(index_g, \
00689 size_psi * sizeof(int)); \
00690 new_psi = (double*)realloc(new_psi, \
00691 size_psi * sizeof(double)); \
00692 } \
00693 } \
00694 \
00695 MACRO_count__lg_lc; \
00696 \
00697 } \
00698 \
00699 index_f[j] = ix - ix_old; \
00700 ix_old = ix; \
00701 \
00702 } \
00703 \
00704 nfft_free(ths->psi); \
00705 size_psi = ix; \
00706 ths->size_psi = size_psi; \
00707 \
00708 index_g = (int*)realloc(index_g, size_psi * sizeof(int)); \
00709 new_psi = (double*)realloc(new_psi, size_psi * sizeof(double)); \
00710 \
00711 ths->psi = new_psi; \
00712 ths->psi_index_g = index_g; \
00713 ths->psi_index_f = index_f; \
00714 \
00715 } \
00716 }
00717
00718 MACRO_nfct_full_psi(A)
00719 MACRO_nfct_full_psi(T)
00720
00721
00722
00723
00724
00728 void nfct_trafo(nfct_plan *ths)
00729 {
00734 ths->g_hat = ths->g1;
00735 ths->g = ths->g2;
00736
00737
00743 TIC(0)
00744 nfct_D_A(ths);
00745 TOC(0)
00746
00747
00748
00754 TIC(1)
00755 fftw_execute(ths->my_fftw_r2r_plan);
00756 TOC(1)
00757
00758
00759 if(ths->nfct_flags & PRE_FULL_PSI)
00760 nfct_full_psi__A(ths);
00761
00762
00768 TIC(2)
00769 nfct_B_A(ths);
00770 TOC(2)
00771
00772 if(ths->nfct_flags & PRE_FULL_PSI) {
00773 nfft_free(ths->psi_index_g);
00774 nfft_free(ths->psi_index_f);
00775 }
00776
00777 }
00778
00779
00780
00781
00782 void nfct_adjoint(nfct_plan *ths)
00783 {
00788 ths->g_hat = ths->g2;
00789 ths->g = ths->g1;
00790
00791 if(ths->nfct_flags & PRE_FULL_PSI)
00792 nfct_full_psi__T(ths);
00793
00799 TIC(2)
00800 nfct_B_T(ths);
00801 TOC(2)
00802
00803 if(ths->nfct_flags & PRE_FULL_PSI) {
00804 nfft_free(ths->psi_index_g);
00805 nfft_free(ths->psi_index_f);
00806 }
00807
00814 TIC(1)
00815 fftw_execute(ths->my_fftw_r2r_plan);
00816 TOC(1)
00817
00818
00823 TIC(0)
00824 nfct_D_T(ths);
00825 TOC(0)
00826
00827 }
00828
00829
00830
00835 static void nfct_precompute_phi_hut(nfct_plan *ths)
00836 {
00837 int kg[ths->d];
00838 int t;
00840 ths->c_phi_inv = (double**)nfft_malloc(ths->d * sizeof(double*));
00841
00842 for(t = 0; t < ths->d; t++)
00843 {
00844 ths->c_phi_inv[t] = (double*)nfft_malloc(ths->N[t] * sizeof(double));
00845
00846 for(kg[t] = 0; kg[t] < ths->N[t]; kg[t]++)
00847 {
00848 ths->c_phi_inv[t][kg[t]] = MACRO_compute_PHI_HUT_INV;
00849 }
00850 }
00851 }
00852
00853
00854
00855 void nfct_precompute_psi(nfct_plan *ths)
00856 {
00857 int t;
00858 int j;
00859 int lc[ths->d];
00860 int lb[ths->d];
00862 for (t = 0; t < ths->d; t++)
00863 {
00864 for(j = 0; j < ths->M_total; j++)
00865 {
00866
00867 MACRO_nfct__lower_boundary(j, t);
00868
00869 for(lc[t] = 0; lc[t] < NFCT_SUMMANDS; lc[t]++)
00870 ths->psi[(j * ths->d + t) * NFCT_SUMMANDS + lc[t]] = MACRO_compute_PSI;
00871
00872 }
00873 }
00874 }
00875
00876
00877
00878
00879
00880
00881
00882 static void nfct_init_help(nfct_plan *ths)
00883 {
00884 int t;
00886 ths->N_total = nfct_prod_int(ths->N, ths->d);
00887
00888 ths->sigma = (double*)nfft_malloc(ths->d * sizeof(double));
00889
00890 for(t = 0; t < ths->d; t++)
00891 ths->sigma[t] = ((double)(ths->n[t] - 1)) / ths->N[t];
00892
00896 ths->r2r_kind = (fftw_r2r_kind*)nfft_malloc (ths->d * sizeof (fftw_r2r_kind));
00897 for (t = 0; t < ths->d; t++)
00898 ths->r2r_kind[t] = FFTW_REDFT00;
00899
00900
00901 NFCT_WINDOW_HELP_INIT;
00902
00903 if(ths->nfct_flags & MALLOC_X)
00904 ths->x = (double*)nfft_malloc(ths->d * ths->M_total * sizeof(double));
00905
00906 if(ths->nfct_flags & MALLOC_F_HAT)
00907 ths->f_hat = (double*)nfft_malloc(ths->N_total * sizeof(double));
00908
00909 if(ths->nfct_flags & MALLOC_F)
00910 ths->f = (double*)nfft_malloc(ths->M_total * sizeof(double));
00911
00912 if(ths->nfct_flags & PRE_PHI_HUT)
00913 nfct_precompute_phi_hut(ths);
00914
00915
00916 if(ths->nfct_flags & PRE_PSI)
00917 {
00918 ths->psi =
00919 (double*)nfft_malloc(ths->M_total * ths->d * NFCT_SUMMANDS * sizeof(double));
00920
00924 ths->nfct_full_psi_eps = pow((R)10, (R)(-10));
00925 }
00926
00927 if(ths->nfct_flags & FFTW_INIT)
00928 {
00929 ths->g1 =
00930 (double*)nfft_malloc(nfct_prod_int(ths->n, ths->d) * sizeof(double));
00931
00932 if(ths->nfct_flags & FFT_OUT_OF_PLACE)
00933 ths->g2 =
00934 (double*) nfft_malloc(nfct_prod_int(ths->n, ths->d) * sizeof(double));
00935 else
00936 ths->g2 = ths->g1;
00937
00938 ths->my_fftw_r2r_plan =
00939 fftw_plan_r2r(ths->d, ths->n, ths->g1, ths->g2,
00940 ths->r2r_kind, ths->fftw_flags);
00941 }
00942
00943 ths->mv_trafo = (void (*) (void* ))nfct_trafo;
00944 ths->mv_adjoint = (void (*) (void* ))nfct_adjoint;
00945 }
00946
00947
00948
00949
00950
00951 void nfct_init(nfct_plan *ths, int d, int *N, int M_total)
00952 {
00953 int t;
00954
00955 ths->d = d;
00956 ths->M_total = M_total;
00957
00958 ths->N = (int*) nfft_malloc(ths->d * sizeof(int));
00959
00960 for(t = 0;t < d; t++)
00961 ths->N[t] = N[t];
00962
00963 ths->n = (int*) nfft_malloc(ths->d * sizeof(int));
00964
00965 for(t = 0; t < d; t++)
00966 ths->n[t] = nfct_fftw_2N(nfft_next_power_of_2(ths->N[t]));
00967
00968 WINDOW_HELP_ESTIMATE_m;
00969
00970 ths->nfct_flags = NFCT_DEFAULT_FLAGS;
00971 ths->fftw_flags = FFTW_DEFAULT_FLAGS;
00972
00973 nfct_init_help(ths);
00974 }
00975
00976
00977 void nfct_init_m(nfct_plan *ths, int d, int *N, int M_total, int m)
00978 {
00979 int t, n[d];
00980
00981 for(t = 0; t < d; t++)
00982 n[t] = nfct_fftw_2N(nfft_next_power_of_2(N[t]));
00983
00984 nfct_init_guru(ths, d, N, M_total, n, m, NFCT_DEFAULT_FLAGS, FFTW_DEFAULT_FLAGS);
00985 }
00986
00987
00988 void nfct_init_guru(nfct_plan *ths, int d, int *N,
00989 int M_total, int *n, int m,
00990 unsigned nfct_flags, unsigned fftw_flags)
00991 {
00992 int t;
00994 ths->d = d;
00995 ths->M_total = M_total;
00996
00997 ths->N = (int*)nfft_malloc(ths->d * sizeof(int));
00998
00999 for(t = 0; t < d; t++)
01000 ths->N[t] = N[t];
01001
01002 ths->n = (int*)nfft_malloc(ths->d * sizeof(int));
01003
01004 for(t = 0; t < d; t++)
01005 ths->n[t] = n[t];
01006
01007 ths->m = m;
01008
01009 ths->nfct_flags = nfct_flags;
01010 ths->fftw_flags = fftw_flags;
01011
01012 nfct_init_help(ths);
01013 }
01014
01015
01016 void nfct_init_1d(nfct_plan *ths, int N0, int M_total)
01017 {
01018 int N[1];
01019
01020 N[0] = N0;
01021 nfct_init(ths, 1, N, M_total);
01022 }
01023
01024 void nfct_init_2d(nfct_plan *ths, int N0, int N1, int M_total)
01025 {
01026 int N[2];
01027
01028 N[0] = N0;
01029 N[1] = N1;
01030 nfct_init(ths, 2, N, M_total);
01031 }
01032
01033 void nfct_init_3d(nfct_plan *ths, int N0, int N1, int N2, int M_total)
01034 {
01035 int N[3];
01036
01037 N[0] = N0;
01038 N[1] = N1;
01039 N[2] = N2;
01040 nfct_init(ths, 3, N, M_total);
01041 }
01042
01043
01047 void nfct_finalize(nfct_plan *ths)
01048 {
01052 int t;
01053
01054 if(ths->nfct_flags & FFTW_INIT)
01055 {
01056 fftw_destroy_plan(ths->my_fftw_r2r_plan);
01057
01058 if(ths->nfct_flags & FFT_OUT_OF_PLACE)
01059 nfft_free(ths->g2);
01060
01061 nfft_free(ths->g1);
01062 }
01063
01064
01065 if(ths->nfct_flags & PRE_PSI)
01066 {
01067 nfft_free(ths->psi);
01068 }
01069
01070 if(ths->nfct_flags & PRE_PHI_HUT)
01071 {
01072 for(t = 0; t < ths->d; t++)
01073 nfft_free(ths->c_phi_inv[t]);
01074 nfft_free(ths->c_phi_inv);
01075 }
01076
01077 if(ths->nfct_flags & MALLOC_F)
01078 nfft_free(ths->f);
01079
01080 if(ths->nfct_flags & MALLOC_F_HAT)
01081 nfft_free(ths->f_hat);
01082
01083 if(ths->nfct_flags & MALLOC_X)
01084 nfft_free(ths->x);
01085
01086 WINDOW_HELP_FINALIZE;
01087
01088 nfft_free(ths->N);
01089 nfft_free(ths->n);
01090 nfft_free(ths->sigma);
01091
01092 nfft_free(ths->r2r_kind);
01093 }
01094