![]() |
NFFT
3.3.1
|
00001 /* 00002 * Copyright (c) 2002, 2016 Jens Keiner, Stefan Kunis, Daniel Potts 00003 * 00004 * This program is free software; you can redistribute it and/or modify it under 00005 * the terms of the GNU General Public License as published by the Free Software 00006 * Foundation; either version 2 of the License, or (at your option) any later 00007 * version. 00008 * 00009 * This program is distributed in the hope that it will be useful, but WITHOUT 00010 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 00011 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 00012 * details. 00013 * 00014 * You should have received a copy of the GNU General Public License along with 00015 * this program; if not, write to the Free Software Foundation, Inc., 51 00016 * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 00017 */ 00018 00025 #include "config.h" 00026 00027 /* Include standard C headers. */ 00028 #include <math.h> 00029 #include <stdlib.h> 00030 #include <string.h> 00031 #ifdef HAVE_COMPLEX_H 00032 #include <complex.h> 00033 #endif 00034 00035 #ifdef _OPENMP 00036 #include <omp.h> 00037 #endif 00038 00039 /* Include NFFT3 utilities header. */ 00040 00041 /* Include NFFT3 library header. */ 00042 #include "nfft3.h" 00043 00044 #include "infft.h" 00045 00046 /* Include private associated Legendre functions header. */ 00047 #include "legendre.h" 00048 00049 /* Include private API header. */ 00050 #include "api.h" 00051 00052 00062 #define NFSFT_DEFAULT_NFFT_CUTOFF 6 00063 00069 #define NFSFT_DEFAULT_THRESHOLD 1000 00070 00076 #define NFSFT_BREAK_EVEN 5 00077 00084 static struct nfsft_wisdom wisdom = {false,0U,-1,-1,0,0,0,0,0}; 00085 00108 static inline void c2e(nfsft_plan *plan) 00109 { 00110 int k; 00111 int n; 00112 double _Complex last; 00113 double _Complex act; 00114 double _Complex *xp; 00115 double _Complex *xm; 00116 int low; 00117 int up; 00118 int lowe; 00119 int upe; 00121 /* Set the first row to order to zero since it is unused. */ 00122 memset(plan->f_hat_intern,0U,(2*plan->N+2)*sizeof(double _Complex)); 00123 00124 /* Determine lower and upper bounds for loop processing even terms. */ 00125 lowe = -plan->N + (plan->N%2); 00126 upe = -lowe; 00127 00128 /* Process even terms. */ 00129 for (n = lowe; n <= upe; n += 2) 00130 { 00131 /* Compute new coefficients \f$\left(c_k^n\right)_{k=-M,\ldots,M}\f$ from 00132 * old coefficients $\left(b_k^n\right)_{k=0,\ldots,M}$. */ 00133 xm = &(plan->f_hat_intern[NFSFT_INDEX(-1,n,plan)]); 00134 xp = &(plan->f_hat_intern[NFSFT_INDEX(+1,n,plan)]); 00135 for(k = 1; k <= plan->N; k++) 00136 { 00137 *xp *= 0.5; 00138 *xm-- = *xp++; 00139 } 00140 /* Set the first coefficient in the array corresponding to this order to 00141 * zero since it is unused. */ 00142 *xm = 0.0; 00143 } 00144 00145 /* Determine lower and upper bounds for loop processing odd terms. */ 00146 low = -plan->N + (1-plan->N%2); 00147 up = -low; 00148 00149 /* Process odd terms incorporating the additional sine term 00150 * \f$\sin \vartheta\f$. */ 00151 for (n = low; n <= up; n += 2) 00152 { 00153 /* Compute new coefficients \f$\left(c_k^n\right)_{k=-M,\ldots,M}\f$ from 00154 * old coefficients $\left(b_k^n\right)_{k=0,\ldots,M-1}$ incorporating 00155 * the additional term \f$\sin \vartheta\f$. */ 00156 plan->f_hat_intern[NFSFT_INDEX(0,n,plan)] *= 2.0; 00157 xp = &(plan->f_hat_intern[NFSFT_INDEX(-plan->N-1,n,plan)]); 00158 /* Set the first coefficient in the array corresponding to this order to zero 00159 * since it is unused. */ 00160 *xp++ = 0.0; 00161 xm = &(plan->f_hat_intern[NFSFT_INDEX(plan->N,n,plan)]); 00162 last = *xm; 00163 *xm = 0.5 * _Complex_I * (0.5*xm[-1]); 00164 *xp++ = -(*xm--); 00165 for (k = plan->N-1; k > 0; k--) 00166 { 00167 act = *xm; 00168 *xm = 0.5 * _Complex_I * (0.5*(xm[-1] - last)); 00169 *xp++ = -(*xm--); 00170 last = act; 00171 } 00172 *xm = 0.0; 00173 } 00174 } 00175 00186 static inline void c2e_transposed(nfsft_plan *plan) 00187 { 00188 int k; 00189 int n; 00190 double _Complex last; 00191 double _Complex act; 00192 double _Complex *xp; 00193 double _Complex *xm; 00194 int low; 00195 int up; 00196 int lowe; 00197 int upe; 00199 /* Determine lower and upper bounds for loop processing even terms. */ 00200 lowe = -plan->N + (plan->N%2); 00201 upe = -lowe; 00202 00203 /* Process even terms. */ 00204 for (n = lowe; n <= upe; n += 2) 00205 { 00206 /* Compute new coefficients \f$\left(b_k^n\right)_{k=0,\ldots,M}\f$ from 00207 * old coefficients $\left(c_k^n\right)_{k=-M,\ldots,M}$. */ 00208 xm = &(plan->f_hat[NFSFT_INDEX(-1,n,plan)]); 00209 xp = &(plan->f_hat[NFSFT_INDEX(+1,n,plan)]); 00210 for(k = 1; k <= plan->N; k++) 00211 { 00212 *xp += *xm--; 00213 *xp++ *= 0.5;; 00214 } 00215 } 00216 00217 /* Determine lower and upper bounds for loop processing odd terms. */ 00218 low = -plan->N + (1-plan->N%2); 00219 up = -low; 00220 00221 /* Process odd terms. */ 00222 for (n = low; n <= up; n += 2) 00223 { 00224 /* Compute new coefficients \f$\left(b_k^n\right)_{k=0,\ldots,M-1}\f$ from 00225 * old coefficients $\left(c_k^n\right)_{k=0,\ldots,M-1}$. */ 00226 xm = &(plan->f_hat[NFSFT_INDEX(-1,n,plan)]); 00227 xp = &(plan->f_hat[NFSFT_INDEX(+1,n,plan)]); 00228 for(k = 1; k <= plan->N; k++) 00229 { 00230 *xp++ -= *xm--; 00231 } 00232 00233 plan->f_hat[NFSFT_INDEX(0,n,plan)] = 00234 -0.25*_Complex_I*plan->f_hat[NFSFT_INDEX(1,n,plan)]; 00235 last = plan->f_hat[NFSFT_INDEX(1,n,plan)]; 00236 plan->f_hat[NFSFT_INDEX(1,n,plan)] = 00237 -0.25*_Complex_I*plan->f_hat[NFSFT_INDEX(2,n,plan)]; 00238 00239 xp = &(plan->f_hat[NFSFT_INDEX(2,n,plan)]); 00240 for (k = 2; k < plan->N; k++) 00241 { 00242 act = *xp; 00243 *xp = -0.25 * _Complex_I * (xp[1] - last); 00244 xp++; 00245 last = act; 00246 } 00247 *xp = 0.25 * _Complex_I * last; 00248 00249 plan->f_hat[NFSFT_INDEX(0,n,plan)] *= 2.0; 00250 } 00251 } 00252 00253 void nfsft_init(nfsft_plan *plan, int N, int M) 00254 { 00255 /* Call nfsft_init_advanced with default flags. */ 00256 nfsft_init_advanced(plan, N, M, NFSFT_MALLOC_X | NFSFT_MALLOC_F | 00257 NFSFT_MALLOC_F_HAT); 00258 } 00259 00260 void nfsft_init_advanced(nfsft_plan* plan, int N, int M, 00261 unsigned int flags) 00262 { 00263 /* Call nfsft_init_guru with the flags and default NFFT cut-off. */ 00264 nfsft_init_guru(plan, N, M, flags, PRE_PHI_HUT | PRE_PSI | FFTW_INIT | NFFT_OMP_BLOCKWISE_ADJOINT | 00265 FFT_OUT_OF_PLACE, NFSFT_DEFAULT_NFFT_CUTOFF); 00266 } 00267 00268 void nfsft_init_guru(nfsft_plan *plan, int N, int M, unsigned int flags, 00269 unsigned int nfft_flags, int nfft_cutoff) 00270 { 00271 int *nfft_size; /*< NFFT size */ 00272 int *fftw_size; /*< FFTW size */ 00273 00274 /* Save the flags in the plan. */ 00275 plan->flags = flags; 00276 00277 /* Save the bandwidth N and the number of samples M in the plan. */ 00278 plan->N = N; 00279 plan->M_total = M; 00280 00281 /* Calculate the next greater power of two with respect to the bandwidth N 00282 * and the corresponding exponent. */ 00283 //next_power_of_2_exp_int(plan->N,&plan->NPT,&plan->t); 00284 00285 /* Save length of array of Fourier coefficients. Owing to the data layout the 00286 * length is (2N+2)(2N+2) */ 00287 plan->N_total = (2*plan->N+2)*(2*plan->N+2); 00288 00289 /* Allocate memory for auxilliary array of spherical Fourier coefficients, 00290 * if neccesary. */ 00291 if (plan->flags & NFSFT_PRESERVE_F_HAT) 00292 { 00293 plan->f_hat_intern = (double _Complex*) nfft_malloc(plan->N_total* 00294 sizeof(double _Complex)); 00295 } 00296 00297 /* Allocate memory for spherical Fourier coefficients, if neccesary. */ 00298 if (plan->flags & NFSFT_MALLOC_F_HAT) 00299 { 00300 plan->f_hat = (double _Complex*) nfft_malloc(plan->N_total* 00301 sizeof(double _Complex)); 00302 } 00303 00304 /* Allocate memory for samples, if neccesary. */ 00305 if (plan->flags & NFSFT_MALLOC_F) 00306 { 00307 plan->f = (double _Complex*) nfft_malloc(plan->M_total*sizeof(double _Complex)); 00308 } 00309 00310 /* Allocate memory for nodes, if neccesary. */ 00311 if (plan->flags & NFSFT_MALLOC_X) 00312 { 00313 plan->x = (double*) nfft_malloc(plan->M_total*2*sizeof(double)); 00314 } 00315 00316 /* Check if fast algorithm is activated. */ 00317 if (plan->flags & NFSFT_NO_FAST_ALGORITHM) 00318 { 00319 } 00320 else 00321 { 00322 nfft_size = (int*)nfft_malloc(2*sizeof(int)); 00323 fftw_size = (int*)nfft_malloc(2*sizeof(int)); 00324 00326 nfft_size[0] = 2*plan->N+2; 00327 nfft_size[1] = 2*plan->N+2; 00328 fftw_size[0] = 4*plan->N; 00329 fftw_size[1] = 4*plan->N; 00330 00332 nfft_init_guru(&plan->plan_nfft, 2, nfft_size, plan->M_total, fftw_size, 00333 nfft_cutoff, nfft_flags, 00334 FFTW_ESTIMATE | FFTW_DESTROY_INPUT); 00335 00336 /* Assign angle array. */ 00337 plan->plan_nfft.x = plan->x; 00338 /* Assign result array. */ 00339 plan->plan_nfft.f = plan->f; 00340 /* Assign Fourier coefficients array. */ 00341 plan->plan_nfft.f_hat = plan->f_hat; 00342 00345 /* Precompute. */ 00346 //nfft_precompute_one_psi(&plan->plan_nfft); 00347 00348 /* Free auxilliary arrays. */ 00349 nfft_free(nfft_size); 00350 nfft_free(fftw_size); 00351 } 00352 00353 plan->mv_trafo = (void (*) (void* ))nfsft_trafo; 00354 plan->mv_adjoint = (void (*) (void* ))nfsft_adjoint; 00355 } 00356 00357 void nfsft_precompute(int N, double kappa, unsigned int nfsft_flags, 00358 unsigned int fpt_flags) 00359 { 00360 int n; /*< The order n */ 00361 00362 /* Check if already initialized. */ 00363 if (wisdom.initialized == true) 00364 { 00365 return; 00366 } 00367 00368 #ifdef _OPENMP 00369 #pragma omp parallel default(shared) 00370 { 00371 int nthreads = omp_get_num_threads(); 00372 int threadid = omp_get_thread_num(); 00373 #pragma omp single 00374 { 00375 wisdom.nthreads = nthreads; 00376 } 00377 } 00378 #endif 00379 00380 /* Save the precomputation flags. */ 00381 wisdom.flags = nfsft_flags; 00382 00383 /* Compute and save N_max = 2^{\ceil{log_2 N}} as next greater 00384 * power of two with respect to N. */ 00385 X(next_power_of_2_exp_int)(N,&wisdom.N_MAX,&wisdom.T_MAX); 00386 00387 /* Check, if precomputation for direct algorithms needs to be performed. */ 00388 if (wisdom.flags & NFSFT_NO_DIRECT_ALGORITHM) 00389 { 00390 wisdom.alpha = NULL; 00391 wisdom.beta = NULL; 00392 wisdom.gamma = NULL; 00393 } 00394 else 00395 { 00396 /* Allocate memory for three-term recursion coefficients. */ 00397 wisdom.alpha = (double*) nfft_malloc((wisdom.N_MAX+1)*(wisdom.N_MAX+2)* 00398 sizeof(double)); 00399 wisdom.beta = (double*) nfft_malloc((wisdom.N_MAX+1)*(wisdom.N_MAX+2)* 00400 sizeof(double)); 00401 wisdom.gamma = (double*) nfft_malloc((wisdom.N_MAX+1)*(wisdom.N_MAX+2)* 00402 sizeof(double)); 00404 /* Compute three-term recurrence coefficients alpha_k^n, beta_k^n, and 00405 * gamma_k^n. */ 00406 alpha_al_all(wisdom.alpha,wisdom.N_MAX); 00407 beta_al_all(wisdom.beta,wisdom.N_MAX); 00408 gamma_al_all(wisdom.gamma,wisdom.N_MAX); 00409 } 00410 00411 /* Check, if precomputation for fast algorithms needs to be performed. */ 00412 if (wisdom.flags & NFSFT_NO_FAST_ALGORITHM) 00413 { 00414 } 00415 else if (wisdom.N_MAX >= NFSFT_BREAK_EVEN) 00416 { 00417 /* Precompute data for DPT/FPT. */ 00418 00419 /* Check, if recursion coefficients have already been calculated. */ 00420 if (wisdom.alpha != NULL) 00421 { 00422 #ifdef _OPENMP 00423 #pragma omp parallel default(shared) private(n) 00424 { 00425 int nthreads = omp_get_num_threads(); 00426 int threadid = omp_get_thread_num(); 00427 #pragma omp single 00428 { 00429 wisdom.nthreads = nthreads; 00430 wisdom.set_threads = (fpt_set*) nfft_malloc(nthreads*sizeof(fpt_set)); 00431 } 00432 00433 wisdom.set_threads[threadid] = fpt_init(wisdom.N_MAX+1, wisdom.T_MAX, 00434 fpt_flags | FPT_AL_SYMMETRY | FPT_PERSISTENT_DATA); 00435 for (n = 0; n <= wisdom.N_MAX; n++) 00436 fpt_precompute(wisdom.set_threads[threadid],n,&wisdom.alpha[ROW(n)], 00437 &wisdom.beta[ROW(n)], &wisdom.gamma[ROW(n)],n,kappa); 00438 } 00439 00440 #else 00441 /* Use the recursion coefficients to precompute FPT data using persistent 00442 * arrays. */ 00443 wisdom.set = fpt_init(wisdom.N_MAX+1, wisdom.T_MAX, 00444 fpt_flags | FPT_AL_SYMMETRY | FPT_PERSISTENT_DATA); 00445 for (n = 0; n <= wisdom.N_MAX; n++) 00446 { 00447 /*fprintf(stderr,"%d\n",n); 00448 fflush(stderr);*/ 00449 /* Precompute data for FPT transformation for order n. */ 00450 fpt_precompute(wisdom.set,n,&wisdom.alpha[ROW(n)],&wisdom.beta[ROW(n)], 00451 &wisdom.gamma[ROW(n)],n,kappa); 00452 } 00453 #endif 00454 } 00455 else 00456 { 00457 #ifdef _OPENMP 00458 #pragma omp parallel default(shared) private(n) 00459 { 00460 double *alpha, *beta, *gamma; 00461 int nthreads = omp_get_num_threads(); 00462 int threadid = omp_get_thread_num(); 00463 #pragma omp single 00464 { 00465 wisdom.nthreads = nthreads; 00466 wisdom.set_threads = (fpt_set*) nfft_malloc(nthreads*sizeof(fpt_set)); 00467 } 00468 00469 alpha = (double*) nfft_malloc((wisdom.N_MAX+2)*sizeof(double)); 00470 beta = (double*) nfft_malloc((wisdom.N_MAX+2)*sizeof(double)); 00471 gamma = (double*) nfft_malloc((wisdom.N_MAX+2)*sizeof(double)); 00472 wisdom.set_threads[threadid] = fpt_init(wisdom.N_MAX+1, wisdom.T_MAX, 00473 fpt_flags | FPT_AL_SYMMETRY); 00474 00475 for (n = 0; n <= wisdom.N_MAX; n++) 00476 { 00477 alpha_al_row(alpha,wisdom.N_MAX,n); 00478 beta_al_row(beta,wisdom.N_MAX,n); 00479 gamma_al_row(gamma,wisdom.N_MAX,n); 00480 00481 /* Precompute data for FPT transformation for order n. */ 00482 fpt_precompute(wisdom.set_threads[threadid],n,alpha,beta,gamma,n, 00483 kappa); 00484 } 00485 /* Free auxilliary arrays. */ 00486 nfft_free(alpha); 00487 nfft_free(beta); 00488 nfft_free(gamma); 00489 } 00490 #else 00491 /* Allocate memory for three-term recursion coefficients. */ 00492 wisdom.alpha = (double*) nfft_malloc((wisdom.N_MAX+2)*sizeof(double)); 00493 wisdom.beta = (double*) nfft_malloc((wisdom.N_MAX+2)*sizeof(double)); 00494 wisdom.gamma = (double*) nfft_malloc((wisdom.N_MAX+2)*sizeof(double)); 00495 wisdom.set = fpt_init(wisdom.N_MAX+1, wisdom.T_MAX, 00496 fpt_flags | FPT_AL_SYMMETRY); 00497 for (n = 0; n <= wisdom.N_MAX; n++) 00498 { 00499 /*fprintf(stderr,"%d NO_DIRECT\n",n); 00500 fflush(stderr);*/ 00501 /* Compute three-term recurrence coefficients alpha_k^n, beta_k^n, and 00502 * gamma_k^n. */ 00503 alpha_al_row(wisdom.alpha,wisdom.N_MAX,n); 00504 beta_al_row(wisdom.beta,wisdom.N_MAX,n); 00505 gamma_al_row(wisdom.gamma,wisdom.N_MAX,n); 00506 00507 /* Precompute data for FPT transformation for order n. */ 00508 fpt_precompute(wisdom.set,n,wisdom.alpha,wisdom.beta,wisdom.gamma,n, 00509 kappa); 00510 } 00511 /* Free auxilliary arrays. */ 00512 nfft_free(wisdom.alpha); 00513 nfft_free(wisdom.beta); 00514 nfft_free(wisdom.gamma); 00515 #endif 00516 wisdom.alpha = NULL; 00517 wisdom.beta = NULL; 00518 wisdom.gamma = NULL; 00519 } 00520 } 00521 00522 /* Wisdom has been initialised. */ 00523 wisdom.initialized = true; 00524 } 00525 00526 void nfsft_forget(void) 00527 { 00528 /* Check if wisdom has been initialised. */ 00529 if (wisdom.initialized == false) 00530 { 00531 /* Nothing to do. */ 00532 return; 00533 } 00534 00535 /* Check, if precomputation for direct algorithms has been performed. */ 00536 if (wisdom.flags & NFSFT_NO_DIRECT_ALGORITHM) 00537 { 00538 } 00539 else 00540 { 00541 /* Free arrays holding three-term recurrence coefficients. */ 00542 nfft_free(wisdom.alpha); 00543 nfft_free(wisdom.beta); 00544 nfft_free(wisdom.gamma); 00545 wisdom.alpha = NULL; 00546 wisdom.beta = NULL; 00547 wisdom.gamma = NULL; 00548 } 00549 00550 /* Check, if precomputation for fast algorithms has been performed. */ 00551 if (wisdom.flags & NFSFT_NO_FAST_ALGORITHM) 00552 { 00553 } 00554 else if (wisdom.N_MAX >= NFSFT_BREAK_EVEN) 00555 { 00556 #ifdef _OPENMP 00557 int k; 00558 for (k = 0; k < wisdom.nthreads; k++) 00559 fpt_finalize(wisdom.set_threads[k]); 00560 nfft_free(wisdom.set_threads); 00561 #else 00562 /* Free precomputed data for FPT transformation. */ 00563 fpt_finalize(wisdom.set); 00564 #endif 00565 } 00566 00567 /* Wisdom is now uninitialised. */ 00568 wisdom.initialized = false; 00569 } 00570 00571 00572 void nfsft_finalize(nfsft_plan *plan) 00573 { 00574 if (!plan) 00575 return; 00576 00577 /* Finalise the nfft plan. */ 00578 nfft_finalize(&plan->plan_nfft); 00579 00580 /* De-allocate memory for auxilliary array of spherical Fourier coefficients, 00581 * if neccesary. */ 00582 if (plan->flags & NFSFT_PRESERVE_F_HAT) 00583 { 00584 nfft_free(plan->f_hat_intern); 00585 } 00586 00587 /* De-allocate memory for spherical Fourier coefficients, if necessary. */ 00588 if (plan->flags & NFSFT_MALLOC_F_HAT) 00589 { 00590 //fprintf(stderr,"deallocating f_hat\n"); 00591 nfft_free(plan->f_hat); 00592 } 00593 00594 /* De-allocate memory for samples, if neccesary. */ 00595 if (plan->flags & NFSFT_MALLOC_F) 00596 { 00597 //fprintf(stderr,"deallocating f\n"); 00598 nfft_free(plan->f); 00599 } 00600 00601 /* De-allocate memory for nodes, if neccesary. */ 00602 if (plan->flags & NFSFT_MALLOC_X) 00603 { 00604 //fprintf(stderr,"deallocating x\n"); 00605 nfft_free(plan->x); 00606 } 00607 } 00608 00609 void nfsft_trafo_direct(nfsft_plan *plan) 00610 { 00611 int m; /*< The node index */ 00612 int k; /*< The degree k */ 00613 int n; /*< The order n */ 00614 int n_abs; /*< The absolute value of the order n, ie n_abs = |n| */ 00615 double *alpha; /*< Pointer to current three-term recurrence 00616 coefficient alpha_k^n for associated Legendre 00617 functions P_k^n */ 00618 double *gamma; /*< Pointer to current three-term recurrence 00619 coefficient beta_k^n for associated Legendre 00620 functions P_k^n */ 00621 double _Complex *a; /*< Pointer to auxilliary array for Clenshaw algor. */ 00622 double _Complex it1; /*< Auxilliary variable for Clenshaw algorithm */ 00623 double _Complex it2; /*< Auxilliary variable for Clenshaw algorithm */ 00624 double _Complex temp; /*< Auxilliary variable for Clenshaw algorithm */ 00625 double _Complex f_m; /*< The final function value f_m = f(x_m) for a 00626 single node. */ 00627 double stheta; /*< Current angle theta for Clenshaw algorithm */ 00628 double sphi; /*< Current angle phi for Clenshaw algorithm */ 00629 00630 #ifdef MEASURE_TIME 00631 plan->MEASURE_TIME_t[0] = 0.0; 00632 plan->MEASURE_TIME_t[1] = 0.0; 00633 plan->MEASURE_TIME_t[2] = 0.0; 00634 #endif 00635 00636 if (wisdom.flags & NFSFT_NO_DIRECT_ALGORITHM) 00637 { 00638 return; 00639 } 00640 00641 /* Copy spherical Fourier coefficients, if necessary. */ 00642 if (plan->flags & NFSFT_PRESERVE_F_HAT) 00643 { 00644 memcpy(plan->f_hat_intern,plan->f_hat,plan->N_total* 00645 sizeof(double _Complex)); 00646 } 00647 else 00648 { 00649 plan->f_hat_intern = plan->f_hat; 00650 } 00651 00652 /* Check, if we compute with L^2-normalized spherical harmonics. If so, 00653 * multiply spherical Fourier coefficients with corresponding normalization 00654 * weight. */ 00655 if (plan->flags & NFSFT_NORMALIZED) 00656 { 00657 /* Traverse Fourier coefficients array. */ 00658 #ifdef _OPENMP 00659 #pragma omp parallel for default(shared) private(k,n) 00660 #endif 00661 for (k = 0; k <= plan->N; k++) 00662 { 00663 for (n = -k; n <= k; n++) 00664 { 00665 /* Multiply with normalization weight. */ 00666 plan->f_hat_intern[NFSFT_INDEX(k,n,plan)] *= 00667 sqrt((2*k+1)/(4.0*KPI)); 00668 } 00669 } 00670 } 00671 00672 /* Distinguish by bandwidth M. */ 00673 if (plan->N == 0) 00674 { 00675 /* N = 0 */ 00676 00677 /* Constant function */ 00678 for (m = 0; m < plan->M_total; m++) 00679 { 00680 plan->f[m] = plan->f_hat_intern[NFSFT_INDEX(0,0,plan)]; 00681 } 00682 } 00683 else 00684 { 00685 /* N > 0 */ 00686 00687 /* Evaluate 00688 * \sum_{k=0}^N \sum_{n=-k}^k a_k^n P_k^{|n|}(cos theta_m) e^{i n phi_m} 00689 * = \sum_{n=-N}^N \sum_{k=|n|}^N a_k^n P_k^{|n|}(cos theta_m) 00690 * e^{i n phi_m}. 00691 */ 00692 #ifdef _OPENMP 00693 #pragma omp parallel for default(shared) private(m,stheta,sphi,f_m,n,a,n_abs,alpha,gamma,it2,it1,k,temp) 00694 #endif 00695 for (m = 0; m < plan->M_total; m++) 00696 { 00697 /* Scale angle theta from [0,1/2] to [0,pi] and apply cosine. */ 00698 stheta = cos(2.0*KPI*plan->x[2*m+1]); 00699 /* Scale angle phi from [-1/2,1/2] to [-pi,pi]. */ 00700 sphi = 2.0*KPI*plan->x[2*m]; 00701 00702 /* Initialize result for current node. */ 00703 f_m = 0.0; 00704 00705 /* For n = -N,...,N, evaluate 00706 * b_n := \sum_{k=|n|}^N a_k^n P_k^{|n|}(cos theta_m) 00707 * using Clenshaw's algorithm. 00708 */ 00709 for (n = -plan->N; n <= plan->N; n++) 00710 { 00711 /* Get pointer to Fourier coefficients vector for current order n. */ 00712 a = &(plan->f_hat_intern[NFSFT_INDEX(0,n,plan)]); 00713 00714 /* Take absolute value of n. */ 00715 n_abs = abs(n); 00716 00717 /* Get pointers to three-term recurrence coefficients arrays. */ 00718 alpha = &(wisdom.alpha[ROW(n_abs)]); 00719 gamma = &(wisdom.gamma[ROW(n_abs)]); 00720 00721 /* Clenshaw's algorithm */ 00722 it2 = a[plan->N]; 00723 it1 = a[plan->N-1]; 00724 for (k = plan->N; k > n_abs + 1; k--) 00725 { 00726 temp = a[k-2] + it2 * gamma[k]; 00727 it2 = it1 + it2 * alpha[k] * stheta; 00728 it1 = temp; 00729 } 00730 00731 /* Compute final step if neccesary. */ 00732 if (n_abs < plan->N) 00733 { 00734 it2 = it1 + it2 * wisdom.alpha[ROWK(n_abs)+1] * stheta; 00735 } 00736 00737 /* Compute final result by multiplying the fixed part 00738 * gamma_|n| (1-cos^2(theta))^{|n|/2} 00739 * for order n and the exponential part 00740 * e^{i n phi} 00741 * and add the result to f_m. 00742 */ 00743 f_m += it2 * wisdom.gamma[ROW(n_abs)] * 00744 pow(1- stheta * stheta, 0.5*n_abs) * cexp(_Complex_I*n*sphi); 00745 } 00746 00747 /* Write result f_m for current node to array f. */ 00748 plan->f[m] = f_m; 00749 } 00750 } 00751 } 00752 00753 void nfsft_adjoint_direct(nfsft_plan *plan) 00754 { 00755 int m; /*< The node index */ 00756 int k; /*< The degree k */ 00757 int n; /*< The order n */ 00758 int n_abs; /*< The absolute value of the order n, ie n_abs = |n| */ 00759 double *alpha; /*< Pointer to current three-term recurrence 00760 coefficient alpha_k^n for associated Legendre 00761 functions P_k^n */ 00762 double *gamma; /*< Pointer to current three-term recurrence 00763 coefficient beta_k^n for associated Legendre 00764 functions P_k^n */ 00765 double _Complex it1; /*< Auxilliary variable for Clenshaw algorithm */ 00766 double _Complex it2; /*< Auxilliary variable for Clenshaw algorithm */ 00767 double _Complex temp; /*< Auxilliary variable for Clenshaw algorithm */ 00768 double stheta; /*< Current angle theta for Clenshaw algorithm */ 00769 double sphi; /*< Current angle phi for Clenshaw algorithm */ 00770 00771 #ifdef MEASURE_TIME 00772 plan->MEASURE_TIME_t[0] = 0.0; 00773 plan->MEASURE_TIME_t[1] = 0.0; 00774 plan->MEASURE_TIME_t[2] = 0.0; 00775 #endif 00776 00777 if (wisdom.flags & NFSFT_NO_DIRECT_ALGORITHM) 00778 { 00779 return; 00780 } 00781 00782 /* Initialise spherical Fourier coefficients array with zeros. */ 00783 memset(plan->f_hat,0U,plan->N_total*sizeof(double _Complex)); 00784 00785 /* Distinguish by bandwidth N. */ 00786 if (plan->N == 0) 00787 { 00788 /* N == 0 */ 00789 00790 /* Constant function */ 00791 for (m = 0; m < plan->M_total; m++) 00792 { 00793 plan->f_hat[NFSFT_INDEX(0,0,plan)] += plan->f[m]; 00794 } 00795 } 00796 else 00797 { 00798 /* N > 0 */ 00799 00800 #ifdef _OPENMP 00801 /* Traverse all orders n. */ 00802 #pragma omp parallel for default(shared) private(n,n_abs,alpha,gamma,m,stheta,sphi,it2,it1,k,temp) 00803 for (n = -plan->N; n <= plan->N; n++) 00804 { 00805 /* Take absolute value of n. */ 00806 n_abs = abs(n); 00807 00808 /* Get pointers to three-term recurrence coefficients arrays. */ 00809 alpha = &(wisdom.alpha[ROW(n_abs)]); 00810 gamma = &(wisdom.gamma[ROW(n_abs)]); 00811 00812 /* Traverse all nodes. */ 00813 for (m = 0; m < plan->M_total; m++) 00814 { 00815 /* Scale angle theta from [0,1/2] to [0,pi] and apply cosine. */ 00816 stheta = cos(2.0*KPI*plan->x[2*m+1]); 00817 /* Scale angle phi from [-1/2,1/2] to [-pi,pi]. */ 00818 sphi = 2.0*KPI*plan->x[2*m]; 00819 00820 /* Transposed Clenshaw algorithm */ 00821 00822 /* Initial step */ 00823 it1 = plan->f[m] * wisdom.gamma[ROW(n_abs)] * 00824 pow(1 - stheta * stheta, 0.5*n_abs) * cexp(-_Complex_I*n*sphi); 00825 plan->f_hat[NFSFT_INDEX(n_abs,n,plan)] += it1; 00826 it2 = 0.0; 00827 00828 if (n_abs < plan->N) 00829 { 00830 it2 = it1 * wisdom.alpha[ROWK(n_abs)+1] * stheta; 00831 plan->f_hat[NFSFT_INDEX(n_abs+1,n,plan)] += it2; 00832 } 00833 00834 /* Loop for transposed Clenshaw algorithm */ 00835 for (k = n_abs+2; k <= plan->N; k++) 00836 { 00837 temp = it2; 00838 it2 = alpha[k] * stheta * it2 + gamma[k] * it1; 00839 it1 = temp; 00840 plan->f_hat[NFSFT_INDEX(k,n,plan)] += it2; 00841 } 00842 } 00843 } 00844 #else 00845 /* Traverse all nodes. */ 00846 for (m = 0; m < plan->M_total; m++) 00847 { 00848 /* Scale angle theta from [0,1/2] to [0,pi] and apply cosine. */ 00849 stheta = cos(2.0*KPI*plan->x[2*m+1]); 00850 /* Scale angle phi from [-1/2,1/2] to [-pi,pi]. */ 00851 sphi = 2.0*KPI*plan->x[2*m]; 00852 00853 /* Traverse all orders n. */ 00854 for (n = -plan->N; n <= plan->N; n++) 00855 { 00856 /* Take absolute value of n. */ 00857 n_abs = abs(n); 00858 00859 /* Get pointers to three-term recurrence coefficients arrays. */ 00860 alpha = &(wisdom.alpha[ROW(n_abs)]); 00861 gamma = &(wisdom.gamma[ROW(n_abs)]); 00862 00863 /* Transposed Clenshaw algorithm */ 00864 00865 /* Initial step */ 00866 it1 = plan->f[m] * wisdom.gamma[ROW(n_abs)] * 00867 pow(1 - stheta * stheta, 0.5*n_abs) * cexp(-_Complex_I*n*sphi); 00868 plan->f_hat[NFSFT_INDEX(n_abs,n,plan)] += it1; 00869 it2 = 0.0; 00870 00871 if (n_abs < plan->N) 00872 { 00873 it2 = it1 * wisdom.alpha[ROWK(n_abs)+1] * stheta; 00874 plan->f_hat[NFSFT_INDEX(n_abs+1,n,plan)] += it2; 00875 } 00876 00877 /* Loop for transposed Clenshaw algorithm */ 00878 for (k = n_abs+2; k <= plan->N; k++) 00879 { 00880 temp = it2; 00881 it2 = alpha[k] * stheta * it2 + gamma[k] * it1; 00882 it1 = temp; 00883 plan->f_hat[NFSFT_INDEX(k,n,plan)] += it2; 00884 } 00885 } 00886 } 00887 #endif 00888 } 00889 00890 /* Check, if we compute with L^2-normalized spherical harmonics. If so, 00891 * multiply spherical Fourier coefficients with corresponding normalization 00892 * weight. */ 00893 if (plan->flags & NFSFT_NORMALIZED) 00894 { 00895 /* Traverse Fourier coefficients array. */ 00896 #ifdef _OPENMP 00897 #pragma omp parallel for default(shared) private(k,n) 00898 #endif 00899 for (k = 0; k <= plan->N; k++) 00900 { 00901 for (n = -k; n <= k; n++) 00902 { 00903 /* Multiply with normalization weight. */ 00904 plan->f_hat[NFSFT_INDEX(k,n,plan)] *= 00905 sqrt((2*k+1)/(4.0*KPI)); 00906 } 00907 } 00908 } 00909 00910 /* Set unused coefficients to zero. */ 00911 if (plan->flags & NFSFT_ZERO_F_HAT) 00912 { 00913 for (n = -plan->N; n <= plan->N+1; n++) 00914 { 00915 memset(&plan->f_hat[NFSFT_INDEX(-plan->N-1,n,plan)],0U, 00916 (plan->N+1+abs(n))*sizeof(double _Complex)); 00917 } 00918 } 00919 } 00920 00921 void nfsft_trafo(nfsft_plan *plan) 00922 { 00923 int k; /*< The degree k */ 00924 int n; /*< The order n */ 00925 #ifdef MEASURE_TIME 00926 ticks t0, t1; 00927 #endif 00928 #ifdef DEBUG 00929 double t, t_pre, t_nfft, t_fpt, t_c2e, t_norm; 00930 t_pre = 0.0; 00931 t_norm = 0.0; 00932 t_fpt = 0.0; 00933 t_c2e = 0.0; 00934 t_nfft = 0.0; 00935 #endif 00936 00937 #ifdef MEASURE_TIME 00938 plan->MEASURE_TIME_t[0] = 0.0; 00939 plan->MEASURE_TIME_t[1] = 0.0; 00940 plan->MEASURE_TIME_t[2] = 0.0; 00941 #endif 00942 00943 if (wisdom.flags & NFSFT_NO_FAST_ALGORITHM) 00944 { 00945 return; 00946 } 00947 00948 /* Check, if precomputation was done and that the bandwidth N is not too 00949 * big. 00950 */ 00951 if (wisdom.initialized == 0 || plan->N > wisdom.N_MAX) 00952 { 00953 return; 00954 } 00955 00956 /* Check, if slow transformation should be used due to small bandwidth. */ 00957 if (plan->N < NFSFT_BREAK_EVEN) 00958 { 00959 /* Use NDSFT. */ 00960 nfsft_trafo_direct(plan); 00961 } 00962 00963 /* Check for correct value of the bandwidth N. */ 00964 else if (plan->N <= wisdom.N_MAX) 00965 { 00966 /* Copy spherical Fourier coefficients, if necessary. */ 00967 if (plan->flags & NFSFT_PRESERVE_F_HAT) 00968 { 00969 memcpy(plan->f_hat_intern,plan->f_hat,plan->N_total* 00970 sizeof(double _Complex)); 00971 } 00972 else 00973 { 00974 plan->f_hat_intern = plan->f_hat; 00975 } 00976 00977 /* Propagate pointer values to the internal NFFT plan to assure 00978 * consistency. Pointers may have been modified externally. 00979 */ 00980 plan->plan_nfft.x = plan->x; 00981 plan->plan_nfft.f = plan->f; 00982 plan->plan_nfft.f_hat = plan->f_hat_intern; 00983 00984 /* Check, if we compute with L^2-normalized spherical harmonics. If so, 00985 * multiply spherical Fourier coefficients with corresponding normalization 00986 * weight. */ 00987 if (plan->flags & NFSFT_NORMALIZED) 00988 { 00989 /* Traverse Fourier coefficients array. */ 00990 #ifdef _OPENMP 00991 #pragma omp parallel for default(shared) private(k,n) 00992 #endif 00993 for (k = 0; k <= plan->N; k++) 00994 { 00995 for (n = -k; n <= k; n++) 00996 { 00997 /* Multiply with normalization weight. */ 00998 plan->f_hat_intern[NFSFT_INDEX(k,n,plan)] *= 00999 sqrt((2*k+1)/(4.0*KPI)); 01000 } 01001 } 01002 } 01003 01004 #ifdef MEASURE_TIME 01005 t0 = getticks(); 01006 #endif 01007 /* Check, which polynomial transform algorithm should be used. */ 01008 if (plan->flags & NFSFT_USE_DPT) 01009 { 01010 #ifdef _OPENMP 01011 #pragma omp parallel for default(shared) private(n) num_threads(wisdom.nthreads) 01012 for (n = -plan->N; n <= plan->N; n++) 01013 fpt_trafo_direct(wisdom.set_threads[omp_get_thread_num()],abs(n), 01014 &plan->f_hat_intern[NFSFT_INDEX(abs(n),n,plan)], 01015 &plan->f_hat_intern[NFSFT_INDEX(0,n,plan)], 01016 plan->N,0U); 01017 #else 01018 /* Use direct discrete polynomial transform DPT. */ 01019 for (n = -plan->N; n <= plan->N; n++) 01020 { 01021 //fprintf(stderr,"nfsft_trafo: n = %d\n",n); 01022 //fflush(stderr); 01023 fpt_trafo_direct(wisdom.set,abs(n), 01024 &plan->f_hat_intern[NFSFT_INDEX(abs(n),n,plan)], 01025 &plan->f_hat_intern[NFSFT_INDEX(0,n,plan)], 01026 plan->N,0U); 01027 } 01028 #endif 01029 } 01030 else 01031 { 01032 #ifdef _OPENMP 01033 #pragma omp parallel for default(shared) private(n) num_threads(wisdom.nthreads) 01034 for (n = -plan->N; n <= plan->N; n++) 01035 fpt_trafo(wisdom.set_threads[omp_get_thread_num()],abs(n), 01036 &plan->f_hat_intern[NFSFT_INDEX(abs(n),n,plan)], 01037 &plan->f_hat_intern[NFSFT_INDEX(0,n,plan)], 01038 plan->N,0U); 01039 #else 01040 /* Use fast polynomial transform FPT. */ 01041 for (n = -plan->N; n <= plan->N; n++) 01042 { 01043 //fprintf(stderr,"nfsft_trafo: n = %d\n",n); 01044 //fflush(stderr); 01045 fpt_trafo(wisdom.set,abs(n), 01046 &plan->f_hat_intern[NFSFT_INDEX(abs(n),n,plan)], 01047 &plan->f_hat_intern[NFSFT_INDEX(0,n,plan)], 01048 plan->N,0U); 01049 } 01050 #endif 01051 } 01052 #ifdef MEASURE_TIME 01053 t1 = getticks(); 01054 plan->MEASURE_TIME_t[0] = nfft_elapsed_seconds(t1,t0); 01055 #endif 01056 01057 #ifdef MEASURE_TIME 01058 t0 = getticks(); 01059 #endif 01060 /* Convert Chebyshev coefficients to Fourier coefficients. */ 01061 c2e(plan); 01062 #ifdef MEASURE_TIME 01063 t1 = getticks(); 01064 plan->MEASURE_TIME_t[1] = nfft_elapsed_seconds(t1,t0); 01065 #endif 01066 01067 #ifdef MEASURE_TIME 01068 t0 = getticks(); 01069 #endif 01070 /* Check, which nonequispaced discrete Fourier transform algorithm should 01071 * be used. 01072 */ 01073 if (plan->flags & NFSFT_USE_NDFT) 01074 { 01075 /* Use NDFT. */ 01076 nfft_trafo_direct(&plan->plan_nfft); 01077 } 01078 else 01079 { 01080 /* Use NFFT. */ 01081 //fprintf(stderr,"nfsft_adjoint: nfft_trafo\n"); 01082 nfft_trafo_2d(&plan->plan_nfft); 01083 } 01084 #ifdef MEASURE_TIME 01085 t1 = getticks(); 01086 plan->MEASURE_TIME_t[2] = nfft_elapsed_seconds(t1,t0); 01087 #endif 01088 } 01089 } 01090 01091 void nfsft_adjoint(nfsft_plan *plan) 01092 { 01093 int k; /*< The degree k */ 01094 int n; /*< The order n */ 01095 #ifdef MEASURE_TIME 01096 ticks t0, t1; 01097 #endif 01098 01099 #ifdef MEASURE_TIME 01100 plan->MEASURE_TIME_t[0] = 0.0; 01101 plan->MEASURE_TIME_t[1] = 0.0; 01102 plan->MEASURE_TIME_t[2] = 0.0; 01103 #endif 01104 01105 if (wisdom.flags & NFSFT_NO_FAST_ALGORITHM) 01106 { 01107 return; 01108 } 01109 01110 /* Check, if precomputation was done and that the bandwidth N is not too 01111 * big. 01112 */ 01113 if (wisdom.initialized == 0 || plan->N > wisdom.N_MAX) 01114 { 01115 return; 01116 } 01117 01118 /* Check, if slow transformation should be used due to small bandwidth. */ 01119 if (plan->N < NFSFT_BREAK_EVEN) 01120 { 01121 /* Use adjoint NDSFT. */ 01122 nfsft_adjoint_direct(plan); 01123 } 01124 /* Check for correct value of the bandwidth N. */ 01125 else if (plan->N <= wisdom.N_MAX) 01126 { 01127 //fprintf(stderr,"nfsft_adjoint: Starting\n"); 01128 //fflush(stderr); 01129 /* Propagate pointer values to the internal NFFT plan to assure 01130 * consistency. Pointers may have been modified externally. 01131 */ 01132 plan->plan_nfft.x = plan->x; 01133 plan->plan_nfft.f = plan->f; 01134 plan->plan_nfft.f_hat = plan->f_hat; 01135 01136 #ifdef MEASURE_TIME 01137 t0 = getticks(); 01138 #endif 01139 /* Check, which adjoint nonequispaced discrete Fourier transform algorithm 01140 * should be used. 01141 */ 01142 if (plan->flags & NFSFT_USE_NDFT) 01143 { 01144 //fprintf(stderr,"nfsft_adjoint: Executing nfft_adjoint_direct\n"); 01145 //fflush(stderr); 01146 /* Use adjoint NDFT. */ 01147 nfft_adjoint_direct(&plan->plan_nfft); 01148 } 01149 else 01150 { 01151 //fprintf(stderr,"nfsft_adjoint: Executing nfft_adjoint\n"); 01152 //fflush(stderr); 01153 //fprintf(stderr,"nfsft_adjoint: nfft_adjoint\n"); 01154 /* Use adjoint NFFT. */ 01155 nfft_adjoint_2d(&plan->plan_nfft); 01156 } 01157 #ifdef MEASURE_TIME 01158 t1 = getticks(); 01159 plan->MEASURE_TIME_t[2] = nfft_elapsed_seconds(t1,t0); 01160 #endif 01161 01162 //fprintf(stderr,"nfsft_adjoint: Executing c2e_transposed\n"); 01163 //fflush(stderr); 01164 #ifdef MEASURE_TIME 01165 t0 = getticks(); 01166 #endif 01167 /* Convert Fourier coefficients to Chebyshev coefficients. */ 01168 c2e_transposed(plan); 01169 #ifdef MEASURE_TIME 01170 t1 = getticks(); 01171 plan->MEASURE_TIME_t[1] = nfft_elapsed_seconds(t1,t0); 01172 #endif 01173 01174 #ifdef MEASURE_TIME 01175 t0 = getticks(); 01176 #endif 01177 /* Check, which transposed polynomial transform algorithm should be used */ 01178 if (plan->flags & NFSFT_USE_DPT) 01179 { 01180 #ifdef _OPENMP 01181 #pragma omp parallel for default(shared) private(n) num_threads(wisdom.nthreads) 01182 for (n = -plan->N; n <= plan->N; n++) 01183 fpt_transposed_direct(wisdom.set_threads[omp_get_thread_num()],abs(n), 01184 &plan->f_hat[NFSFT_INDEX(abs(n),n,plan)], 01185 &plan->f_hat[NFSFT_INDEX(0,n,plan)], 01186 plan->N,0U); 01187 #else 01188 /* Use transposed DPT. */ 01189 for (n = -plan->N; n <= plan->N; n++) 01190 { 01191 //fprintf(stderr,"nfsft_adjoint: Executing dpt_transposed\n"); 01192 //fflush(stderr); 01193 fpt_transposed_direct(wisdom.set,abs(n), 01194 &plan->f_hat[NFSFT_INDEX(abs(n),n,plan)], 01195 &plan->f_hat[NFSFT_INDEX(0,n,plan)], 01196 plan->N,0U); 01197 } 01198 #endif 01199 } 01200 else 01201 { 01202 #ifdef _OPENMP 01203 #pragma omp parallel for default(shared) private(n) num_threads(wisdom.nthreads) 01204 for (n = -plan->N; n <= plan->N; n++) 01205 fpt_transposed(wisdom.set_threads[omp_get_thread_num()],abs(n), 01206 &plan->f_hat[NFSFT_INDEX(abs(n),n,plan)], 01207 &plan->f_hat[NFSFT_INDEX(0,n,plan)], 01208 plan->N,0U); 01209 #else 01210 //fprintf(stderr,"nfsft_adjoint: fpt_transposed\n"); 01211 /* Use transposed FPT. */ 01212 for (n = -plan->N; n <= plan->N; n++) 01213 { 01214 //fprintf(stderr,"nfsft_adjoint: Executing fpt_transposed\n"); 01215 //fflush(stderr); 01216 fpt_transposed(wisdom.set,abs(n), 01217 &plan->f_hat[NFSFT_INDEX(abs(n),n,plan)], 01218 &plan->f_hat[NFSFT_INDEX(0,n,plan)], 01219 plan->N,0U); 01220 } 01221 #endif 01222 } 01223 #ifdef MEASURE_TIME 01224 t1 = getticks(); 01225 plan->MEASURE_TIME_t[0] = nfft_elapsed_seconds(t1,t0); 01226 #endif 01227 01228 /* Check, if we compute with L^2-normalized spherical harmonics. If so, 01229 * multiply spherical Fourier coefficients with corresponding normalization 01230 * weight. */ 01231 if (plan->flags & NFSFT_NORMALIZED) 01232 { 01233 //fprintf(stderr,"nfsft_adjoint: Normalizing\n"); 01234 //fflush(stderr); 01235 /* Traverse Fourier coefficients array. */ 01236 #ifdef _OPENMP 01237 #pragma omp parallel for default(shared) private(k,n) 01238 #endif 01239 for (k = 0; k <= plan->N; k++) 01240 { 01241 for (n = -k; n <= k; n++) 01242 { 01243 /* Multiply with normalization weight. */ 01244 plan->f_hat[NFSFT_INDEX(k,n,plan)] *= 01245 sqrt((2*k+1)/(4.0*KPI)); 01246 } 01247 } 01248 } 01249 01250 /* Set unused coefficients to zero. */ 01251 if (plan->flags & NFSFT_ZERO_F_HAT) 01252 { 01253 //fprintf(stderr,"nfsft_adjoint: Setting to zero\n"); 01254 //fflush(stderr); 01255 for (n = -plan->N; n <= plan->N+1; n++) 01256 { 01257 memset(&plan->f_hat[NFSFT_INDEX(-plan->N-1,n,plan)],0U, 01258 (plan->N+1+abs(n))*sizeof(double _Complex)); 01259 } 01260 } 01261 //fprintf(stderr,"nfsft_adjoint: Finished\n"); 01262 //fflush(stderr); 01263 } 01264 } 01265 01266 void nfsft_precompute_x(nfsft_plan *plan) 01267 { 01268 /* Pass angle array to NFFT plan. */ 01269 plan->plan_nfft.x = plan->x; 01270 01271 /* Precompute. */ 01272 if (plan->plan_nfft.flags & PRE_ONE_PSI) 01273 nfft_precompute_one_psi(&plan->plan_nfft); 01274 } 01275 /* \} */