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