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