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 NFST_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 NFST_SUMMANDS ( 2 * ths->m + 2)
00034 #define NODE(p,r) ( ths->x[(p) * ths->d + (r)])
00035
00036 #define MACRO_ndst_init_result_trafo \
00037 memset( f, 0, ths->M_total * sizeof( double));
00038 #define MACRO_ndst_init_result_adjoint \
00039 memset( f_hat, 0, ths->N_total * sizeof( double));
00040
00041
00042 #define MACRO_nfst_D_init_result_A \
00043 memset(g_hat, 0, nfst_prod_minus_a_int( ths->n, 0, ths->d) * sizeof( double));
00044 #define MACRO_nfst_D_init_result_T \
00045 memset(f_hat, 0, ths->N_total * sizeof( double));
00046
00047 #define MACRO_nfst_B_init_result_A \
00048 memset(f, 0, ths->M_total * sizeof( double));
00049 #define MACRO_nfst_B_init_result_T \
00050 memset(g, 0, nfst_prod_minus_a_int( ths->n, 0, ths->d) * sizeof( double));
00051
00052
00053 #define NFST_PRE_WINFUN( d) ths->N[d] = 2 * ths->N[d]; \
00054 ths->n[d] = nfst_fftw_2N( ths->n[d]);
00055
00056 #define NFST_POST_WINFUN( d) ths->N[d] = ((int)(0.5 * ths->N[d])); \
00057 ths->n[d] = nfst_fftw_2N_rev( ths->n[d]);
00058
00059
00060 #define NFST_WINDOW_HELP_INIT WINDOW_HELP_INIT
00061
00062
00063 double nfst_phi_hut( nfst_plan *ths, int k, int d)
00064 {
00065 NFST_PRE_WINFUN( d);
00066 double phi_hut_tmp = PHI_HUT( k, d);
00067 NFST_POST_WINFUN( d);
00068
00069 return phi_hut_tmp;
00070 }
00071
00072 double nfst_phi( nfst_plan *ths, double x, int d)
00073 {
00074 NFST_PRE_WINFUN( d);
00075 double phi_tmp = PHI( x, d);
00076 NFST_POST_WINFUN( d);
00077
00078 return phi_tmp;
00079 }
00080
00081 int nfst_fftw_2N( int n)
00082 {
00083 return 2 * ( n + 1);
00084 }
00085
00086 int nfst_fftw_2N_rev( int n)
00087 {
00088 div_t n_div;
00089
00090 n_div = div(n, 2);
00091 return n_div.quot - 1;
00092 }
00093
00094 #define MACRO_with_sin_vec sin_vec[t][ka[t]]
00095 #define MACRO_without_sin_vec sin( 2.0 * PI * (ka[t]+1) * NODE(j,t))
00096
00097
00098 #define MACRO_with_PRE_PHI_HUT ths->c_phi_inv[t][kg[t]];
00099 #define MACRO_compute_PHI_HUT_INV (1.0 / (nfst_phi_hut( ths, kg[t]+1, t)))
00100
00101 #define MACRO_with_PRE_PSI ths->psi[(j * ths->d + t) * NFST_SUMMANDS + lc[t]];
00102 #define MACRO_compute_PSI \
00103 nfst_phi( ths, NODE(j,t) - (( double)(lc[t] + lb[t])) / nfst_fftw_2N( ths->n[t]), t)
00104
00105
00106
00122 #define MACRO_ndst_malloc__sin_vec \
00123 \
00124 double **sin_vec; \
00125 sin_vec = (double**)malloc( ths->d * sizeof( double*)); \
00126 for( t = 0; t < ths->d; t++) \
00127 sin_vec[t] = (double*)malloc( ( ths->N[t] - 1) * sizeof( double)); \
00128
00129
00130
00131
00132 #define MACRO_ndst_free__sin_vec \
00133 { \
00134 \
00135 for( t = 0; t < ths->d; t++) \
00136 free( sin_vec[t]); \
00137 free( sin_vec); \
00138 }
00139
00140
00141
00142 #define MACRO_ndst_init__sin_vec \
00143 { \
00144 for( t = 0; t < ths->d; t++) \
00145 { \
00146 cos_x[t] = cos( 2.0 * PI * NODE(j,t)); \
00147 sin_vec[t][0] = sin( 2.0 * PI * NODE(j,t)); \
00148 sin_vec[t][1] = sin( 4.0 * PI * NODE(j,t)); \
00149 for( k = 2; k < ths->N[t] - 1; k++) \
00150 sin_vec[t][k] = 2.0 * cos_x[t] * sin_vec[t][k-1] \
00151 - sin_vec[t][k-2]; \
00152 } \
00153 }
00154
00155
00156 #define MACRO_ndst_init__k__sin_k( which_one) \
00157 { \
00158 sin_k[0] = 1.0; \
00159 for( t = 0; t < ths->d; t++) \
00160 ka[t] = 0; \
00161 \
00162 for( t = 0; t < ths->d; t++) \
00163 { \
00164 sin_k[t+1] = sin_k[t] * MACRO_ ##which_one; \
00165 } \
00166 }
00167
00168
00169 #define MACRO_ndst_count__k__sin_k( which_one) \
00170 { \
00171 ka[ths->d-1]++; \
00172 i = ths->d - 1; \
00173 while( ( ka[i] == ths->N[i] - 1) && ( i > 0)) \
00174 { \
00175 ka[i - 1]++; \
00176 ka[i] = 0; \
00177 \
00178 i--; \
00179 } \
00180 for( t = i; t < ths->d; t++) \
00181 sin_k[t+1] = sin_k[t] * MACRO_ ##which_one; \
00182 }
00183
00184
00185 #define MACRO_ndst_compute__trafo \
00186 { \
00187 f[j] += f_hat[k] * sin_k[ths->d]; \
00188 }
00189
00190 #define MACRO_ndst_compute__adjoint \
00191 { \
00192 f_hat[k] += f[j] * sin_k[ths->d]; \
00193 }
00194
00195
00196
00197 #define MACRO_ndst( which_one) \
00198 void ndst_ ## which_one ( nfst_plan *ths) \
00199 { \
00200 int j, k, t, i; \
00201 int ka[ths->d]; \
00202 double sin_k[ths->d+1]; \
00203 double cos_x[ths->d]; \
00204 \
00205 double *f = ths->f; \
00206 double *f_hat = ths->f_hat; \
00207 \
00208 MACRO_ndst_init_result_ ## which_one; \
00209 \
00210 if( ths->d == 1) \
00211 for( j = 0; j < ths->M_total; j++) \
00212 for( k = 0; k < ths->N_total; k++) \
00213 { \
00214 sin_k[ths->d] = sin( 2.0 * PI * (k+1) * NODE(j,0)); \
00215 MACRO_ndst_compute__ ## which_one; \
00216 } \
00217 else \
00218 if( 1 == 0) \
00219 \
00220 for( j = 0; j < ths->M_total; j++) \
00221 { \
00222 MACRO_ndst_init__k__sin_k(without_sin_vec); \
00223 \
00224 for( k = 0; k < ths->N_total; k++) \
00225 { \
00226 MACRO_ndst_compute__ ## which_one; \
00227 \
00228 MACRO_ndst_count__k__sin_k(without_sin_vec); \
00229 } \
00230 } \
00231 else \
00232 { \
00233 \
00234 MACRO_ndst_malloc__sin_vec; \
00235 \
00236 for( j = 0; j < ths->M_total; j++) \
00237 { \
00238 MACRO_ndst_init__sin_vec; \
00239 \
00240 MACRO_ndst_init__k__sin_k(with_sin_vec); \
00241 \
00242 for( k = 0; k < ths->N_total; k++) \
00243 { \
00244 MACRO_ndst_compute__ ## which_one; \
00245 \
00246 MACRO_ndst_count__k__sin_k(with_sin_vec); \
00247 } \
00248 } \
00249 MACRO_ndst_free__sin_vec; \
00250 } \
00251 }
00252
00253
00254 MACRO_ndst(trafo)
00255 MACRO_ndst(adjoint)
00256
00257
00258
00259
00274 #define MACRO_nfst__lower_boundary( j,act_dim) \
00275 { \
00276 lb[(act_dim)] = \
00277 ((int)(NODE((j),(act_dim)) * nfst_fftw_2N( ths->n[(act_dim)]))) - ths->m; \
00278 }
00279
00280 #define MACRO_nfst_D_compute_A \
00281 { \
00282 g_hat[kg_plain[ths->d]] = f_hat[k_L] * c_phi_inv_k[ths->d]; \
00283 }
00284
00285 #define MACRO_nfst_D_compute_T \
00286 { \
00287 f_hat[k_L] = g_hat[kg_plain[ths->d]] * c_phi_inv_k[ths->d]; \
00288 }
00289
00290
00291 #define MACRO_init__kg \
00292 { \
00293 for( t = 0; t < ths->d; t++) \
00294 kg[t] = 0; \
00295 \
00296 i = 0; \
00297 }
00298
00299
00300 #define MACRO_count__kg \
00301 { \
00302 kg[ths->d - 1]++; \
00303 i = ths->d - 1; \
00304 while( ( kg[i] == ths->N[i] - 1) && ( i > 0)) \
00305 { \
00306 kg[i - 1]++; \
00307 kg[i] = 0; \
00308 \
00309 i--; \
00310 } \
00311 }
00312
00313
00314 #define MACRO_update__c_phi_inv_k__lg_plain( which_one, which_phi) \
00315 { \
00316 for( t = i; t < ths->d; t++) { \
00317 MACRO__c_phi_inv_k( which_phi); \
00318 kg_plain[t+1] = kg_plain[t] * ths->n[t] + kg[t]; \
00319 } \
00320 }
00321
00322
00323 #define MACRO__c_phi_inv_k( which_phi) \
00324 { \
00325 c_phi_inv_k[t+1] = 0.5 * c_phi_inv_k[t] * MACRO_ ## which_phi; \
00326 }
00327
00328
00329 #define MACRO_nfst_D(which_one) \
00330 inline void nfst_D_ ## which_one (nfst_plan *ths) \
00331 { \
00332 int k_L; \
00333 \
00334 int i, t; \
00335 int kg[ths->d]; \
00336 double c_phi_inv_k[ths->d+1]; \
00337 int kg_plain[ths->d+1]; \
00338 \
00339 double *g_hat, *f_hat; \
00340 \
00341 g_hat = ths->g_hat; \
00342 f_hat = ths->f_hat; \
00343 \
00344 MACRO_nfst_D_init_result_ ## which_one \
00345 \
00346 c_phi_inv_k[0] = 1; \
00347 kg_plain[0] = 0; \
00348 \
00349 MACRO_init__kg; \
00350 \
00351 if( ths->nfst_flags & PRE_PHI_HUT) \
00352 \
00353 for( k_L = 0; k_L < ths->N_total; k_L++) \
00354 { \
00355 MACRO_update__c_phi_inv_k__lg_plain( which_one, with_PRE_PHI_HUT); \
00356 \
00357 MACRO_nfst_D_compute_ ## which_one; \
00358 \
00359 MACRO_count__kg; \
00360 \
00361 } \
00362 \
00363 else \
00364 \
00365 for( k_L = 0; k_L < ths->N_total; k_L++) \
00366 { \
00367 MACRO_update__c_phi_inv_k__lg_plain( which_one, compute_PHI_HUT_INV); \
00368 \
00369 MACRO_nfst_D_compute_ ## which_one; \
00370 \
00371 MACRO_count__kg \
00372 \
00373 } \
00374 }
00375
00376 MACRO_nfst_D(A)
00377 MACRO_nfst_D(T)
00378
00379
00380
00381
00382
00383
00384
00388 #define MACRO_nfst_B_PRE_FULL_PSI_compute_A \
00389 { \
00390 (*fj) += ths->psi[ix] * g[ths->psi_index_g[ix]]; \
00391 }
00392
00393 #define MACRO_nfst_B_PRE_FULL_PSI_compute_T \
00394 { \
00395 g[ths->psi_index_g[ix]] += ths->psi[ix] * (*fj); \
00396 }
00397
00398
00399
00400 #define MACRO_nfst_B_compute_A \
00401 { \
00402 (*fj) += phi_tilde[ths->d] * g[lg_plain[ths->d]]; \
00403 }
00404
00405 #define MACRO_nfst_B_compute_T \
00406 { \
00407 g[lg_plain[ths->d]] += phi_tilde[ths->d] * (*fj); \
00408 }
00409
00410
00411
00412 #define MACRO_compute_lg_offset__count_lg( i0) \
00413 { \
00414 \
00415 if( lb[(i0)] < 0) \
00416 { \
00417 lg_offset[(i0)] = \
00418 (lb[(i0)] % nfst_fftw_2N( ths->n[(i0)])) + nfst_fftw_2N( ths->n[(i0)]); \
00419 } \
00420 else \
00421 { \
00422 lg_offset[(i0)] = lb[(i0)] % nfst_fftw_2N( ths->n[(i0)]); \
00423 } \
00424 \
00425 if( lg_offset[(i0)] > ths->n[(i0)]+1) \
00426 lg_offset[(i0)] = -( nfst_fftw_2N( ths->n[(i0)]) - lg_offset[(i0)]); \
00427 }
00428
00429
00430
00431 #define MACRO_set__lg__to__lg_offset \
00432 { \
00433 if( lg_offset[i] <= 0) \
00434 { \
00435 lg[i] = -lg_offset[i]; \
00436 count_lg[i] = -1; \
00437 } \
00438 else \
00439 { \
00440 lg[i] = +lg_offset[i]; \
00441 count_lg[i] = +1; \
00442 } \
00443 }
00444
00445
00446
00447 #define MACRO_count__lg(dim) \
00448 { \
00449 \
00450 if( ((lg[(dim)] == 0) || (lg[(dim)] == (ths->n[(dim)] + 1))) ) \
00451 count_lg[(dim)] *= -1; \
00452 \
00453 lg[(dim)] += count_lg[(dim)]; \
00454 }
00455
00456
00457
00458 #define MACRO_init_lb_lg_lc_phi_tilde_lg_plain( which_psi) \
00459 { \
00460 for( i = 0; i < ths->d; i++) \
00461 { \
00462 MACRO_nfst__lower_boundary( j, i); \
00463 \
00464 MACRO_compute_lg_offset__count_lg( i); \
00465 MACRO_set__lg__to__lg_offset; \
00466 \
00467 \
00468 lc[i] = 0; \
00469 } \
00470 \
00471 for( t = 0; t < ths->d; t++) \
00472 { \
00473 if( lg[t] == 0) \
00474 { \
00475 lg_plain[t+1] = lg_plain[t] * ths->n[t]; \
00476 phi_tilde[t+1] = 0.0; \
00477 } \
00478 else \
00479 if( lg[t] == ths->n[t]+1) \
00480 { \
00481 lg_plain[t+1] = lg_plain[t] * ths->n[t] + ths->n[t]-1; \
00482 phi_tilde[t+1] = 0.0; \
00483 } \
00484 else \
00485 { \
00486 MACRO__phi_tilde( which_psi); \
00487 lg_plain[t+1] = lg_plain[t] * ths->n[t] + lg[t]-1; \
00488 } \
00489 } \
00490 \
00491 i = 0; \
00492 }
00493
00494
00495
00496 #define MACRO_count__lg_lc \
00497 { \
00498 MACRO_count__lg( ths->d-1); \
00499 \
00500 lc[ths->d - 1]++; \
00501 i = ths->d - 1; \
00502 \
00503 while( (lc[i] == NFST_SUMMANDS) && (i > 0)) \
00504 { \
00505 lc[i - 1]++; \
00506 lc[i] = 0; \
00507 \
00508 \
00509 MACRO_count__lg( i - 1); \
00510 \
00511 MACRO_set__lg__to__lg_offset; \
00512 \
00513 i--; \
00514 } \
00515 }
00516
00517
00518 #define MACRO_update__phi_tilde__lg_plain( which_psi) \
00519 { \
00520 for( t = i; t < ths->d; t++) \
00521 { \
00522 if( (lg[t] != 0) && (lg[t] != ths->n[t]+1)) \
00523 { \
00524 MACRO__phi_tilde( which_psi); \
00525 lg_plain[t+1] = lg_plain[t] * ths->n[t] + lg[t]-1; \
00526 } \
00527 else \
00528 phi_tilde[t+1] = 0.0; \
00529 } \
00530 }
00531
00532
00533
00534 #define MACRO__phi_tilde( which_psi) \
00535 { \
00536 phi_tilde[t+1] = (double)count_lg[t] * phi_tilde[t] * MACRO_ ## which_psi; \
00537 }
00538
00539
00540
00541
00542 #define MACRO_nfst_B( which_one) \
00543 inline void nfst_B_ ## which_one ( nfst_plan *ths) \
00544 { \
00545 int lb[ths->d]; \
00546 int j, t, i; \
00547 int lprod, l_L, ix; \
00548 int lc[ths->d]; \
00549 int lg[ths->d]; \
00550 int lg_offset[ths->d]; \
00551 int count_lg[ths->d]; \
00552 int lg_plain[ths->d+1]; \
00553 double *f, *g; \
00554 double phi_tilde[ths->d+1]; \
00555 double *fj; \
00556 \
00557 f = ths->f; g = ths->g; \
00558 \
00559 MACRO_nfst_B_init_result_ ## which_one \
00560 \
00561 \
00562 if( (ths->nfst_flags & PRE_PSI) && (ths->nfst_flags & PRE_FULL_PSI)) \
00563 { \
00564 for( ix = 0, j = 0, fj = &f[0]; j < ths->M_total; j++, fj += 1) \
00565 for( l_L = 0; l_L < ths->psi_index_f[j]; l_L++, ix++) \
00566 { \
00567 MACRO_nfst_B_PRE_FULL_PSI_compute_ ## which_one; \
00568 } \
00569 } \
00570 else \
00571 { \
00572 phi_tilde[0] = 1; \
00573 lg_plain[0] = 0; \
00574 \
00575 for( t = 0, lprod = 1; t < ths->d; t++) \
00576 lprod *= NFST_SUMMANDS; \
00577 \
00578 \
00579 if( ths->nfst_flags & PRE_PSI) \
00580 { \
00581 for( j = 0, fj = &f[0]; j < ths->M_total; j++, fj += 1) \
00582 { \
00583 MACRO_init_lb_lg_lc_phi_tilde_lg_plain( with_PRE_PSI); \
00584 \
00585 for( l_L = 0; l_L < lprod; l_L++) \
00586 { \
00587 MACRO_update__phi_tilde__lg_plain( with_PRE_PSI); \
00588 \
00589 MACRO_nfst_B_compute_ ## which_one; \
00590 \
00591 MACRO_count__lg_lc; \
00592 \
00593 } \
00594 } \
00595 } \
00596 \
00597 \
00598 else \
00599 { \
00600 for( j = 0, fj = &f[0]; j < ths->M_total; j++, fj += 1) \
00601 { \
00602 MACRO_init_lb_lg_lc_phi_tilde_lg_plain( compute_PSI); \
00603 \
00604 for( l_L = 0; l_L < lprod; l_L++) \
00605 { \
00606 MACRO_update__phi_tilde__lg_plain( compute_PSI); \
00607 \
00608 MACRO_nfst_B_compute_ ## which_one; \
00609 \
00610 MACRO_count__lg_lc; \
00611 \
00612 } \
00613 } \
00614 } \
00615 } \
00616 }
00617
00618 MACRO_nfst_B(A)
00619 MACRO_nfst_B(T)
00620
00621
00622
00623
00624
00625
00630 void nfst_trafo( nfst_plan *ths)
00631 {
00636 ths->g_hat = ths->g1;
00637 ths->g = ths->g2;
00638
00639
00645 TIC(0)
00646 nfst_D_A( ths);
00647 TOC(0)
00648
00649
00650
00656 TIC(1)
00657 fftw_execute( ths->my_fftw_r2r_plan);
00658 TOC(1)
00659
00660
00661
00666 TIC(2)
00667 nfst_B_A( ths);
00668 TOC(2)
00669
00670 }
00671
00672
00673
00674
00675 void nfst_adjoint( nfst_plan *ths)
00676 {
00681 ths->g_hat = ths->g2;
00682 ths->g = ths->g1;
00683
00684
00690 TIC(2)
00691 nfst_B_T( ths);
00692 TOC(2)
00693
00694
00695
00701 TIC(1)
00702 fftw_execute( ths->my_fftw_r2r_plan);
00703 TOC(1)
00704
00705
00706
00711 TIC(0)
00712 nfst_D_T( ths);
00713 TOC(0)
00714
00715 }
00716
00717
00718
00723 void nfst_precompute_phi_hut( nfst_plan *ths)
00724 {
00725 int kg[ths->d];
00726 int t;
00728 ths->c_phi_inv = (double**)fftw_malloc( ths->d * sizeof( double*));
00729
00730 for( t = 0; t < ths->d; t++)
00731 {
00732 ths->c_phi_inv[t] = (double*)fftw_malloc( ( ths->N[t] - 1) * sizeof( double));
00733
00734 for( kg[t] = 0; kg[t] < ths->N[t] - 1; kg[t]++)
00735 {
00736 ths->c_phi_inv[t][kg[t]] = MACRO_compute_PHI_HUT_INV;
00737 }
00738 }
00739 }
00740
00741
00742
00743 void nfst_precompute_psi( nfst_plan *ths)
00744 {
00745 int t;
00746 int j;
00747 int lc[ths->d];
00748 int lb[ths->d];
00750 for (t = 0; t < ths->d; t++)
00751 {
00752 for(j = 0; j < ths->M_total; j++)
00753 {
00754 MACRO_nfst__lower_boundary( j, t);
00755
00756 for( lc[t] = 0; lc[t] < NFST_SUMMANDS; lc[t]++)
00757 ths->psi[(j * ths->d + t) * NFST_SUMMANDS + lc[t]] = MACRO_compute_PSI;
00758
00759 }
00760 }
00761
00762
00763 if ( ths->nfst_flags & PRE_FULL_PSI)
00764 nfst_full_psi( ths, ths->nfst_full_psi_eps);
00765
00766 }
00767
00768
00769
00771 void nfst_full_psi(nfst_plan *ths, double eps)
00772 {
00773 int t, i;
00774 int j;
00775 int l_L;
00776 int lc[ths->d];
00777 int lg_plain[ths->d+1];
00778 int count_lg[ths->d];
00779 int lg_offset[ths->d];
00780 int lg[ths->d];
00781 int lprod;
00782 int lb[ths->d];
00784 double phi_tilde[ths->d+1];
00785
00786 int *index_g, *index_f;
00787 double *new_psi;
00788 int ix, ix_old, size_psi;
00789
00790 phi_tilde[0] = 1.0;
00791 lg_plain[0] = 0;
00792
00793 if(ths->nfst_flags & PRE_PSI)
00794 {
00795 size_psi = ths->M_total;
00796 index_f = (int*)malloc( ths->M_total * sizeof( int));
00797 index_g = (int*)malloc( size_psi * sizeof( int));
00798 new_psi = (double*)malloc( size_psi * sizeof( double));
00799
00800 for( t = 0,lprod = 1; t < ths->d; t++)
00801 {
00802 lprod *= NFST_SUMMANDS;
00803 eps *= PHI( 0, t);
00804 }
00805
00806 for( ix = 0, ix_old = 0, j = 0; j < ths->M_total; j++)
00807 {
00808 MACRO_init_lb_lg_lc_phi_tilde_lg_plain( with_PRE_PSI);
00809
00810 for( l_L = 0; l_L < lprod; l_L++)
00811 {
00812 MACRO_update__phi_tilde__lg_plain( with_PRE_PSI);
00813
00814 if( fabs(phi_tilde[ths->d]) > eps)
00815 {
00816 index_g[ix] = lg_plain[ths->d];
00817 new_psi[ix] = phi_tilde[ths->d];
00818
00819 ix++;
00820 if( ix == size_psi)
00821 {
00822 size_psi += ths->M_total;
00823 index_g = (int*)realloc( index_g, size_psi * sizeof( int));
00824 new_psi = (double*)realloc( new_psi, size_psi * sizeof( double));
00825 }
00826 }
00827 MACRO_count__lg_lc;
00828
00829 }
00830
00831 index_f[j] = ix - ix_old;
00832 ix_old = ix;
00833
00834 }
00835
00836 free( ths->psi);
00837
00838 size_psi = ix;
00839 ths->size_psi = size_psi;
00840 index_g = (int*)realloc( index_g, size_psi * sizeof( int));
00841 new_psi = (double*)realloc( new_psi, size_psi * sizeof( double));
00842
00843 ths->psi = new_psi;
00844 ths->psi_index_g = index_g;
00845 ths->psi_index_f = index_f;
00846
00847 }
00848 }
00849
00850
00851
00852
00853 void nfst_init_help( nfst_plan *ths)
00854 {
00855 int t;
00857 ths->N_total = nfst_prod_minus_a_int( ths->N, 1, ths->d);
00858
00859 ths->sigma = (double*)fftw_malloc( ths->d * sizeof( double));
00860
00861 for( t = 0; t < ths->d; t++)
00862
00863 ths->sigma[t] = ((double)ths->n[t] + 1) / ths->N[t];
00864
00865
00866 ths->r2r_kind = (fftw_r2r_kind*) fftw_malloc ( ths->d * sizeof( fftw_r2r_kind));
00867 for (t = 0; t < ths->d; t++)
00868 ths->r2r_kind[t] = FFTW_RODFT00;
00869
00870
00871 WINDOW_HELP_INIT;
00872
00873 if(ths->nfst_flags & MALLOC_X)
00874 ths->x = (double*)fftw_malloc( ths->d * ths->M_total * sizeof( double));
00875
00876 if(ths->nfst_flags & MALLOC_F_HAT)
00877 ths->f_hat = (double*)fftw_malloc( ths->N_total * sizeof( double));
00878
00879 if(ths->nfst_flags & MALLOC_F)
00880 ths->f = (double*)fftw_malloc( ths->M_total * sizeof( double));
00881
00882 if(ths->nfst_flags & PRE_PHI_HUT)
00883 nfst_precompute_phi_hut( ths);
00884
00885
00886 if(ths->nfst_flags & PRE_PSI)
00887 {
00888 ths->psi =
00889 (double*)malloc( ths->M_total * ths->d * NFST_SUMMANDS * sizeof( double));
00890
00894 ths->nfst_full_psi_eps = pow(10, -10);
00895 }
00896
00897 if(ths->nfst_flags & FFTW_INIT)
00898 {
00899 ths->g1 =
00900 (double*)fftw_malloc( nfst_prod_minus_a_int( ths->n, 0, ths->d) * sizeof( double));
00901
00902 if(ths->nfst_flags & FFT_OUT_OF_PLACE)
00903 ths->g2 =
00904 (double*)fftw_malloc( nfst_prod_minus_a_int( ths->n, 0, ths->d) * sizeof( double));
00905 else
00906 ths->g2 = ths->g1;
00907
00908 ths->my_fftw_r2r_plan =
00909 fftw_plan_r2r( ths->d, ths->n, ths->g1, ths->g2, ths->r2r_kind, ths->fftw_flags);
00910 }
00911 }
00912
00913 void nfst_init( nfst_plan *ths, int d, int *N, int M_total)
00914 {
00915 int t;
00917 ths->d = d;
00918
00919 ths->N = (int*)fftw_malloc( ths->d * sizeof( int));
00920
00921 for(t = 0;t < d; t++)
00922 ths->N[t] = N[t];
00923
00924 ths->n = (int*)fftw_malloc( ths->d * sizeof( int));
00925
00926 for( t = 0; t < d; t++)
00927 ths->n[t] = 2 * nfft_next_power_of_2( ths->N[t]) - 1;
00928
00929 ths->M_total = M_total;
00930
00931 WINDOW_HELP_ESTIMATE_m;
00932
00933 ths->nfst_flags = NFST_DEFAULT_FLAGS;
00934 ths->fftw_flags = FFTW_DEFAULT_FLAGS;
00935
00936 nfst_init_help( ths);
00937 }
00938
00939
00940 void nfst_init_m( nfst_plan *ths, int d, int *N, int M_total, int m)
00941 {
00942 int t, n[d];
00943
00944 for( t = 0; t < d; t++)
00945 n[t] = nfst_fftw_2N( nfft_next_power_of_2( N[t]));
00946
00947 nfst_init_guru( ths, d, N, M_total, n, m, NFST_DEFAULT_FLAGS, FFTW_DEFAULT_FLAGS);
00948 }
00949
00950
00951 void nfst_init_guru( nfst_plan *ths, int d, int *N,
00952 int M_total, int *n, int m,
00953 unsigned nfst_flags, unsigned fftw_flags)
00954 {
00955 int t;
00957 ths->d = d;
00958 ths->M_total = M_total;
00959
00960 ths->N = (int*)fftw_malloc( ths->d * sizeof( int));
00961
00962 for( t = 0; t < d; t++)
00963 ths->N[t] = N[t];
00964
00965 ths->n = (int*)fftw_malloc( ths->d * sizeof( int));
00966
00967 for( t = 0; t < d; t++)
00968 ths->n[t] = n[t];
00969
00970 ths->m = m;
00971
00972 ths->nfst_flags = nfst_flags;
00973 ths->fftw_flags = fftw_flags;
00974
00975 nfst_init_help( ths);
00976 }
00977
00978
00979 void nfst_init_1d( nfst_plan *ths, int N0, int M_total)
00980 {
00981 int N[1];
00982
00983 N[0] = N0;
00984 nfst_init( ths, 1, N, M_total);
00985 }
00986
00987 void nfst_init_2d( nfst_plan *ths, int N0, int N1, int M_total)
00988 {
00989 int N[2];
00990
00991 N[0] = N0;
00992 N[1] = N1;
00993 nfst_init( ths, 2, N, M_total);
00994 }
00995
00996 void nfst_init_3d( nfst_plan *ths, int N0, int N1, int N2, int M_total)
00997 {
00998 int N[3];
00999
01000 N[0] = N0;
01001 N[1] = N1;
01002 N[2] = N2;
01003 nfst_init( ths, 3, N, M_total);
01004 }
01005
01006 void nfst_finalize( nfst_plan *ths)
01007 {
01008 int t;
01009
01010 if( ths->nfst_flags & FFTW_INIT)
01011 {
01012 fftw_destroy_plan( ths->my_fftw_r2r_plan);
01013
01014 if( ths->nfst_flags & FFT_OUT_OF_PLACE)
01015 fftw_free( ths->g2);
01016
01017 fftw_free( ths->g1);
01018 }
01019
01020
01021 if( ths->nfst_flags & PRE_PSI)
01022 {
01023 if( ths->nfst_flags & PRE_FULL_PSI)
01024 {
01025 free( ths->psi_index_g);
01026 free( ths->psi_index_f);
01027 }
01028
01029 free( ths->psi);
01030 }
01031
01032 if( ths->nfst_flags & PRE_PHI_HUT) {
01033 for( t = 0; t < ths->d; t++)
01034 fftw_free( ths->c_phi_inv[t]);
01035 fftw_free( ths->c_phi_inv);
01036 }
01037
01038 if( ths->nfst_flags & MALLOC_F)
01039 fftw_free( ths->f);
01040
01041 if( ths->nfst_flags & MALLOC_F_HAT)
01042 fftw_free( ths->f_hat);
01043
01044 if( ths->nfst_flags & MALLOC_X)
01045 fftw_free( ths->x);
01046
01047 WINDOW_HELP_FINALIZE;
01048
01049 fftw_free( ths->N);
01050 fftw_free( ths->n);
01051 fftw_free( ths->sigma);
01052
01053 }
01054
