NFFT  3.3.1
nfsft.c
Go to the documentation of this file.
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 /* \} */