NFFT  3.3.0
nfsft.c
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2002, 2015 Jens Keiner, Stefan Kunis, Daniel Potts
3  *
4  * This program is free software; you can redistribute it and/or modify it under
5  * the terms of the GNU General Public License as published by the Free Software
6  * Foundation; either version 2 of the License, or (at your option) any later
7  * version.
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12  * details.
13  *
14  * You should have received a copy of the GNU General Public License along with
15  * this program; if not, write to the Free Software Foundation, Inc., 51
16  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  */
18 
19 /* $Id$ */
20 
27 #include "config.h"
28 
29 /* Include standard C headers. */
30 #include <math.h>
31 #include <stdlib.h>
32 #include <string.h>
33 #ifdef HAVE_COMPLEX_H
34 #include <complex.h>
35 #endif
36 
37 #ifdef _OPENMP
38 #include <omp.h>
39 #endif
40 
41 /* Include NFFT3 utilities header. */
42 
43 /* Include NFFT3 library header. */
44 #include "nfft3.h"
45 
46 #include "infft.h"
47 
48 /* Include private associated Legendre functions header. */
49 #include "legendre.h"
50 
51 /* Include private API header. */
52 #include "api.h"
53 
54 
64 #define NFSFT_DEFAULT_NFFT_CUTOFF 6
65 
71 #define NFSFT_DEFAULT_THRESHOLD 1000
72 
78 #define NFSFT_BREAK_EVEN 5
79 
86 static struct nfsft_wisdom wisdom = {false,0U,-1,-1,0,0,0,0,0};
87 
110 static inline void c2e(nfsft_plan *plan)
111 {
112  int k;
113  int n;
114  double _Complex last;
115  double _Complex act;
116  double _Complex *xp;
117  double _Complex *xm;
118  int low;
119  int up;
120  int lowe;
121  int upe;
123  /* Set the first row to order to zero since it is unused. */
124  memset(plan->f_hat_intern,0U,(2*plan->N+2)*sizeof(double _Complex));
125 
126  /* Determine lower and upper bounds for loop processing even terms. */
127  lowe = -plan->N + (plan->N%2);
128  upe = -lowe;
129 
130  /* Process even terms. */
131  for (n = lowe; n <= upe; n += 2)
132  {
133  /* Compute new coefficients \f$\left(c_k^n\right)_{k=-M,\ldots,M}\f$ from
134  * old coefficients $\left(b_k^n\right)_{k=0,\ldots,M}$. */
135  xm = &(plan->f_hat_intern[NFSFT_INDEX(-1,n,plan)]);
136  xp = &(plan->f_hat_intern[NFSFT_INDEX(+1,n,plan)]);
137  for(k = 1; k <= plan->N; k++)
138  {
139  *xp *= 0.5;
140  *xm-- = *xp++;
141  }
142  /* Set the first coefficient in the array corresponding to this order to
143  * zero since it is unused. */
144  *xm = 0.0;
145  }
146 
147  /* Determine lower and upper bounds for loop processing odd terms. */
148  low = -plan->N + (1-plan->N%2);
149  up = -low;
150 
151  /* Process odd terms incorporating the additional sine term
152  * \f$\sin \vartheta\f$. */
153  for (n = low; n <= up; n += 2)
154  {
155  /* Compute new coefficients \f$\left(c_k^n\right)_{k=-M,\ldots,M}\f$ from
156  * old coefficients $\left(b_k^n\right)_{k=0,\ldots,M-1}$ incorporating
157  * the additional term \f$\sin \vartheta\f$. */
158  plan->f_hat_intern[NFSFT_INDEX(0,n,plan)] *= 2.0;
159  xp = &(plan->f_hat_intern[NFSFT_INDEX(-plan->N-1,n,plan)]);
160  /* Set the first coefficient in the array corresponding to this order to zero
161  * since it is unused. */
162  *xp++ = 0.0;
163  xm = &(plan->f_hat_intern[NFSFT_INDEX(plan->N,n,plan)]);
164  last = *xm;
165  *xm = 0.5 * _Complex_I * (0.5*xm[-1]);
166  *xp++ = -(*xm--);
167  for (k = plan->N-1; k > 0; k--)
168  {
169  act = *xm;
170  *xm = 0.5 * _Complex_I * (0.5*(xm[-1] - last));
171  *xp++ = -(*xm--);
172  last = act;
173  }
174  *xm = 0.0;
175  }
176 }
177 
188 static inline void c2e_transposed(nfsft_plan *plan)
189 {
190  int k;
191  int n;
192  double _Complex last;
193  double _Complex act;
194  double _Complex *xp;
195  double _Complex *xm;
196  int low;
197  int up;
198  int lowe;
199  int upe;
201  /* Determine lower and upper bounds for loop processing even terms. */
202  lowe = -plan->N + (plan->N%2);
203  upe = -lowe;
204 
205  /* Process even terms. */
206  for (n = lowe; n <= upe; n += 2)
207  {
208  /* Compute new coefficients \f$\left(b_k^n\right)_{k=0,\ldots,M}\f$ from
209  * old coefficients $\left(c_k^n\right)_{k=-M,\ldots,M}$. */
210  xm = &(plan->f_hat[NFSFT_INDEX(-1,n,plan)]);
211  xp = &(plan->f_hat[NFSFT_INDEX(+1,n,plan)]);
212  for(k = 1; k <= plan->N; k++)
213  {
214  *xp += *xm--;
215  *xp++ *= 0.5;;
216  }
217  }
218 
219  /* Determine lower and upper bounds for loop processing odd terms. */
220  low = -plan->N + (1-plan->N%2);
221  up = -low;
222 
223  /* Process odd terms. */
224  for (n = low; n <= up; n += 2)
225  {
226  /* Compute new coefficients \f$\left(b_k^n\right)_{k=0,\ldots,M-1}\f$ from
227  * old coefficients $\left(c_k^n\right)_{k=0,\ldots,M-1}$. */
228  xm = &(plan->f_hat[NFSFT_INDEX(-1,n,plan)]);
229  xp = &(plan->f_hat[NFSFT_INDEX(+1,n,plan)]);
230  for(k = 1; k <= plan->N; k++)
231  {
232  *xp++ -= *xm--;
233  }
234 
235  plan->f_hat[NFSFT_INDEX(0,n,plan)] =
236  -0.25*_Complex_I*plan->f_hat[NFSFT_INDEX(1,n,plan)];
237  last = plan->f_hat[NFSFT_INDEX(1,n,plan)];
238  plan->f_hat[NFSFT_INDEX(1,n,plan)] =
239  -0.25*_Complex_I*plan->f_hat[NFSFT_INDEX(2,n,plan)];
240 
241  xp = &(plan->f_hat[NFSFT_INDEX(2,n,plan)]);
242  for (k = 2; k < plan->N; k++)
243  {
244  act = *xp;
245  *xp = -0.25 * _Complex_I * (xp[1] - last);
246  xp++;
247  last = act;
248  }
249  *xp = 0.25 * _Complex_I * last;
250 
251  plan->f_hat[NFSFT_INDEX(0,n,plan)] *= 2.0;
252  }
253 }
254 
255 void nfsft_init(nfsft_plan *plan, int N, int M)
256 {
257  /* Call nfsft_init_advanced with default flags. */
258  nfsft_init_advanced(plan, N, M, NFSFT_MALLOC_X | NFSFT_MALLOC_F |
259  NFSFT_MALLOC_F_HAT);
260 }
261 
262 void nfsft_init_advanced(nfsft_plan* plan, int N, int M,
263  unsigned int flags)
264 {
265  /* Call nfsft_init_guru with the flags and default NFFT cut-off. */
266  nfsft_init_guru(plan, N, M, flags, PRE_PHI_HUT | PRE_PSI | FFTW_INIT | NFFT_OMP_BLOCKWISE_ADJOINT |
267  FFT_OUT_OF_PLACE, NFSFT_DEFAULT_NFFT_CUTOFF);
268 }
269 
270 void nfsft_init_guru(nfsft_plan *plan, int N, int M, unsigned int flags,
271  unsigned int nfft_flags, int nfft_cutoff)
272 {
273  int *nfft_size; /*< NFFT size */
274  int *fftw_size; /*< FFTW size */
275 
276  /* Save the flags in the plan. */
277  plan->flags = flags;
278 
279  /* Save the bandwidth N and the number of samples M in the plan. */
280  plan->N = N;
281  plan->M_total = M;
282 
283  /* Calculate the next greater power of two with respect to the bandwidth N
284  * and the corresponding exponent. */
285  //next_power_of_2_exp(plan->N,&plan->NPT,&plan->t);
286 
287  /* Save length of array of Fourier coefficients. Owing to the data layout the
288  * length is (2N+2)(2N+2) */
289  plan->N_total = (2*plan->N+2)*(2*plan->N+2);
290 
291  /* Allocate memory for auxilliary array of spherical Fourier coefficients,
292  * if neccesary. */
293  if (plan->flags & NFSFT_PRESERVE_F_HAT)
294  {
295  plan->f_hat_intern = (double _Complex*) nfft_malloc(plan->N_total*
296  sizeof(double _Complex));
297  }
298 
299  /* Allocate memory for spherical Fourier coefficients, if neccesary. */
300  if (plan->flags & NFSFT_MALLOC_F_HAT)
301  {
302  plan->f_hat = (double _Complex*) nfft_malloc(plan->N_total*
303  sizeof(double _Complex));
304  }
305 
306  /* Allocate memory for samples, if neccesary. */
307  if (plan->flags & NFSFT_MALLOC_F)
308  {
309  plan->f = (double _Complex*) nfft_malloc(plan->M_total*sizeof(double _Complex));
310  }
311 
312  /* Allocate memory for nodes, if neccesary. */
313  if (plan->flags & NFSFT_MALLOC_X)
314  {
315  plan->x = (double*) nfft_malloc(plan->M_total*2*sizeof(double));
316  }
317 
318  /* Check if fast algorithm is activated. */
319  if (plan->flags & NFSFT_NO_FAST_ALGORITHM)
320  {
321  }
322  else
323  {
324  nfft_size = (int*)nfft_malloc(2*sizeof(int));
325  fftw_size = (int*)nfft_malloc(2*sizeof(int));
326 
328  nfft_size[0] = 2*plan->N+2;
329  nfft_size[1] = 2*plan->N+2;
330  fftw_size[0] = 4*plan->N;
331  fftw_size[1] = 4*plan->N;
332 
334  nfft_init_guru(&plan->plan_nfft, 2, nfft_size, plan->M_total, fftw_size,
335  nfft_cutoff, nfft_flags,
336  FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
337 
338  /* Assign angle array. */
339  plan->plan_nfft.x = plan->x;
340  /* Assign result array. */
341  plan->plan_nfft.f = plan->f;
342  /* Assign Fourier coefficients array. */
343  plan->plan_nfft.f_hat = plan->f_hat;
344 
347  /* Precompute. */
348  //nfft_precompute_one_psi(&plan->plan_nfft);
349 
350  /* Free auxilliary arrays. */
351  nfft_free(nfft_size);
352  nfft_free(fftw_size);
353  }
354 
355  plan->mv_trafo = (void (*) (void* ))nfsft_trafo;
356  plan->mv_adjoint = (void (*) (void* ))nfsft_adjoint;
357 }
358 
359 void nfsft_precompute(int N, double kappa, unsigned int nfsft_flags,
360  unsigned int fpt_flags)
361 {
362  int n; /*< The order n */
363 
364  /* Check if already initialized. */
365  if (wisdom.initialized == true)
366  {
367  return;
368  }
369 
370 #ifdef _OPENMP
371  #pragma omp parallel default(shared)
372  {
373  int nthreads = omp_get_num_threads();
374  int threadid = omp_get_thread_num();
375  #pragma omp single
376  {
377  wisdom.nthreads = nthreads;
378  }
379  }
380 #endif
381 
382  /* Save the precomputation flags. */
383  wisdom.flags = nfsft_flags;
384 
385  /* Compute and save N_max = 2^{\ceil{log_2 N}} as next greater
386  * power of two with respect to N. */
387  X(next_power_of_2_exp)(N,&wisdom.N_MAX,&wisdom.T_MAX);
388 
389  /* Check, if precomputation for direct algorithms needs to be performed. */
390  if (wisdom.flags & NFSFT_NO_DIRECT_ALGORITHM)
391  {
392  wisdom.alpha = NULL;
393  wisdom.beta = NULL;
394  wisdom.gamma = NULL;
395  }
396  else
397  {
398  /* Allocate memory for three-term recursion coefficients. */
399  wisdom.alpha = (double*) nfft_malloc((wisdom.N_MAX+1)*(wisdom.N_MAX+2)*
400  sizeof(double));
401  wisdom.beta = (double*) nfft_malloc((wisdom.N_MAX+1)*(wisdom.N_MAX+2)*
402  sizeof(double));
403  wisdom.gamma = (double*) nfft_malloc((wisdom.N_MAX+1)*(wisdom.N_MAX+2)*
404  sizeof(double));
406  /* Compute three-term recurrence coefficients alpha_k^n, beta_k^n, and
407  * gamma_k^n. */
408  alpha_al_all(wisdom.alpha,wisdom.N_MAX);
409  beta_al_all(wisdom.beta,wisdom.N_MAX);
410  gamma_al_all(wisdom.gamma,wisdom.N_MAX);
411  }
412 
413  /* Check, if precomputation for fast algorithms needs to be performed. */
414  if (wisdom.flags & NFSFT_NO_FAST_ALGORITHM)
415  {
416  }
417  else if (wisdom.N_MAX >= NFSFT_BREAK_EVEN)
418  {
419  /* Precompute data for DPT/FPT. */
420 
421  /* Check, if recursion coefficients have already been calculated. */
422  if (wisdom.alpha != NULL)
423  {
424 #ifdef _OPENMP
425  #pragma omp parallel default(shared) private(n)
426  {
427  int nthreads = omp_get_num_threads();
428  int threadid = omp_get_thread_num();
429  #pragma omp single
430  {
431  wisdom.nthreads = nthreads;
432  wisdom.set_threads = (fpt_set*) nfft_malloc(nthreads*sizeof(fpt_set));
433  }
434 
435  wisdom.set_threads[threadid] = fpt_init(wisdom.N_MAX+1, wisdom.T_MAX,
436  fpt_flags | FPT_AL_SYMMETRY | FPT_PERSISTENT_DATA);
437  for (n = 0; n <= wisdom.N_MAX; n++)
438  fpt_precompute(wisdom.set_threads[threadid],n,&wisdom.alpha[ROW(n)],
439  &wisdom.beta[ROW(n)], &wisdom.gamma[ROW(n)],n,kappa);
440  }
441 
442 #else
443  /* Use the recursion coefficients to precompute FPT data using persistent
444  * arrays. */
445  wisdom.set = fpt_init(wisdom.N_MAX+1, wisdom.T_MAX,
446  fpt_flags | FPT_AL_SYMMETRY | FPT_PERSISTENT_DATA);
447  for (n = 0; n <= wisdom.N_MAX; n++)
448  {
449  /*fprintf(stderr,"%d\n",n);
450  fflush(stderr);*/
451  /* Precompute data for FPT transformation for order n. */
452  fpt_precompute(wisdom.set,n,&wisdom.alpha[ROW(n)],&wisdom.beta[ROW(n)],
453  &wisdom.gamma[ROW(n)],n,kappa);
454  }
455 #endif
456  }
457  else
458  {
459 #ifdef _OPENMP
460  #pragma omp parallel default(shared) private(n)
461  {
462  double *alpha, *beta, *gamma;
463  int nthreads = omp_get_num_threads();
464  int threadid = omp_get_thread_num();
465  #pragma omp single
466  {
467  wisdom.nthreads = nthreads;
468  wisdom.set_threads = (fpt_set*) nfft_malloc(nthreads*sizeof(fpt_set));
469  }
470 
471  alpha = (double*) nfft_malloc((wisdom.N_MAX+2)*sizeof(double));
472  beta = (double*) nfft_malloc((wisdom.N_MAX+2)*sizeof(double));
473  gamma = (double*) nfft_malloc((wisdom.N_MAX+2)*sizeof(double));
474  wisdom.set_threads[threadid] = fpt_init(wisdom.N_MAX+1, wisdom.T_MAX,
475  fpt_flags | FPT_AL_SYMMETRY);
476 
477  for (n = 0; n <= wisdom.N_MAX; n++)
478  {
479  alpha_al_row(alpha,wisdom.N_MAX,n);
480  beta_al_row(beta,wisdom.N_MAX,n);
481  gamma_al_row(gamma,wisdom.N_MAX,n);
482 
483  /* Precompute data for FPT transformation for order n. */
484  fpt_precompute(wisdom.set_threads[threadid],n,alpha,beta,gamma,n,
485  kappa);
486  }
487  /* Free auxilliary arrays. */
488  nfft_free(alpha);
489  nfft_free(beta);
490  nfft_free(gamma);
491  }
492 #else
493  /* Allocate memory for three-term recursion coefficients. */
494  wisdom.alpha = (double*) nfft_malloc((wisdom.N_MAX+2)*sizeof(double));
495  wisdom.beta = (double*) nfft_malloc((wisdom.N_MAX+2)*sizeof(double));
496  wisdom.gamma = (double*) nfft_malloc((wisdom.N_MAX+2)*sizeof(double));
497  wisdom.set = fpt_init(wisdom.N_MAX+1, wisdom.T_MAX,
498  fpt_flags | FPT_AL_SYMMETRY);
499  for (n = 0; n <= wisdom.N_MAX; n++)
500  {
501  /*fprintf(stderr,"%d NO_DIRECT\n",n);
502  fflush(stderr);*/
503  /* Compute three-term recurrence coefficients alpha_k^n, beta_k^n, and
504  * gamma_k^n. */
505  alpha_al_row(wisdom.alpha,wisdom.N_MAX,n);
506  beta_al_row(wisdom.beta,wisdom.N_MAX,n);
507  gamma_al_row(wisdom.gamma,wisdom.N_MAX,n);
508 
509  /* Precompute data for FPT transformation for order n. */
510  fpt_precompute(wisdom.set,n,wisdom.alpha,wisdom.beta,wisdom.gamma,n,
511  kappa);
512  }
513  /* Free auxilliary arrays. */
514  nfft_free(wisdom.alpha);
515  nfft_free(wisdom.beta);
516  nfft_free(wisdom.gamma);
517 #endif
518  wisdom.alpha = NULL;
519  wisdom.beta = NULL;
520  wisdom.gamma = NULL;
521  }
522  }
523 
524  /* Wisdom has been initialised. */
525  wisdom.initialized = true;
526 }
527 
528 void nfsft_forget(void)
529 {
530  /* Check if wisdom has been initialised. */
531  if (wisdom.initialized == false)
532  {
533  /* Nothing to do. */
534  return;
535  }
536 
537  /* Check, if precomputation for direct algorithms has been performed. */
538  if (wisdom.flags & NFSFT_NO_DIRECT_ALGORITHM)
539  {
540  }
541  else
542  {
543  /* Free arrays holding three-term recurrence coefficients. */
544  nfft_free(wisdom.alpha);
545  nfft_free(wisdom.beta);
546  nfft_free(wisdom.gamma);
547  wisdom.alpha = NULL;
548  wisdom.beta = NULL;
549  wisdom.gamma = NULL;
550  }
551 
552  /* Check, if precomputation for fast algorithms has been performed. */
553  if (wisdom.flags & NFSFT_NO_FAST_ALGORITHM)
554  {
555  }
556  else if (wisdom.N_MAX >= NFSFT_BREAK_EVEN)
557  {
558 #ifdef _OPENMP
559  int k;
560  for (k = 0; k < wisdom.nthreads; k++)
561  fpt_finalize(wisdom.set_threads[k]);
562  nfft_free(wisdom.set_threads);
563 #else
564  /* Free precomputed data for FPT transformation. */
565  fpt_finalize(wisdom.set);
566 #endif
567  }
568 
569  /* Wisdom is now uninitialised. */
570  wisdom.initialized = false;
571 }
572 
573 
574 void nfsft_finalize(nfsft_plan *plan)
575 {
576  if (!plan)
577  return;
578 
579  /* Finalise the nfft plan. */
580  nfft_finalize(&plan->plan_nfft);
581 
582  /* De-allocate memory for auxilliary array of spherical Fourier coefficients,
583  * if neccesary. */
584  if (plan->flags & NFSFT_PRESERVE_F_HAT)
585  {
586  nfft_free(plan->f_hat_intern);
587  }
588 
589  /* De-allocate memory for spherical Fourier coefficients, if necessary. */
590  if (plan->flags & NFSFT_MALLOC_F_HAT)
591  {
592  //fprintf(stderr,"deallocating f_hat\n");
593  nfft_free(plan->f_hat);
594  }
595 
596  /* De-allocate memory for samples, if neccesary. */
597  if (plan->flags & NFSFT_MALLOC_F)
598  {
599  //fprintf(stderr,"deallocating f\n");
600  nfft_free(plan->f);
601  }
602 
603  /* De-allocate memory for nodes, if neccesary. */
604  if (plan->flags & NFSFT_MALLOC_X)
605  {
606  //fprintf(stderr,"deallocating x\n");
607  nfft_free(plan->x);
608  }
609 }
610 
611 void nfsft_trafo_direct(nfsft_plan *plan)
612 {
613  int m; /*< The node index */
614  int k; /*< The degree k */
615  int n; /*< The order n */
616  int n_abs; /*< The absolute value of the order n, ie n_abs = |n| */
617  double *alpha; /*< Pointer to current three-term recurrence
618  coefficient alpha_k^n for associated Legendre
619  functions P_k^n */
620  double *gamma; /*< Pointer to current three-term recurrence
621  coefficient beta_k^n for associated Legendre
622  functions P_k^n */
623  double _Complex *a; /*< Pointer to auxilliary array for Clenshaw algor. */
624  double _Complex it1; /*< Auxilliary variable for Clenshaw algorithm */
625  double _Complex it2; /*< Auxilliary variable for Clenshaw algorithm */
626  double _Complex temp; /*< Auxilliary variable for Clenshaw algorithm */
627  double _Complex f_m; /*< The final function value f_m = f(x_m) for a
628  single node. */
629  double stheta; /*< Current angle theta for Clenshaw algorithm */
630  double sphi; /*< Current angle phi for Clenshaw algorithm */
631 
632 #ifdef MEASURE_TIME
633  plan->MEASURE_TIME_t[0] = 0.0;
634  plan->MEASURE_TIME_t[1] = 0.0;
635  plan->MEASURE_TIME_t[2] = 0.0;
636 #endif
637 
638  if (wisdom.flags & NFSFT_NO_DIRECT_ALGORITHM)
639  {
640  return;
641  }
642 
643  /* Copy spherical Fourier coefficients, if necessary. */
644  if (plan->flags & NFSFT_PRESERVE_F_HAT)
645  {
646  memcpy(plan->f_hat_intern,plan->f_hat,plan->N_total*
647  sizeof(double _Complex));
648  }
649  else
650  {
651  plan->f_hat_intern = plan->f_hat;
652  }
653 
654  /* Check, if we compute with L^2-normalized spherical harmonics. If so,
655  * multiply spherical Fourier coefficients with corresponding normalization
656  * weight. */
657  if (plan->flags & NFSFT_NORMALIZED)
658  {
659  /* Traverse Fourier coefficients array. */
660 #ifdef _OPENMP
661  #pragma omp parallel for default(shared) private(k,n)
662 #endif
663  for (k = 0; k <= plan->N; k++)
664  {
665  for (n = -k; n <= k; n++)
666  {
667  /* Multiply with normalization weight. */
668  plan->f_hat_intern[NFSFT_INDEX(k,n,plan)] *=
669  sqrt((2*k+1)/(4.0*KPI));
670  }
671  }
672  }
673 
674  /* Distinguish by bandwidth M. */
675  if (plan->N == 0)
676  {
677  /* N = 0 */
678 
679  /* Constant function */
680  for (m = 0; m < plan->M_total; m++)
681  {
682  plan->f[m] = plan->f_hat_intern[NFSFT_INDEX(0,0,plan)];
683  }
684  }
685  else
686  {
687  /* N > 0 */
688 
689  /* Evaluate
690  * \sum_{k=0}^N \sum_{n=-k}^k a_k^n P_k^{|n|}(cos theta_m) e^{i n phi_m}
691  * = \sum_{n=-N}^N \sum_{k=|n|}^N a_k^n P_k^{|n|}(cos theta_m)
692  * e^{i n phi_m}.
693  */
694 #ifdef _OPENMP
695  #pragma omp parallel for default(shared) private(m,stheta,sphi,f_m,n,a,n_abs,alpha,gamma,it2,it1,k,temp)
696 #endif
697  for (m = 0; m < plan->M_total; m++)
698  {
699  /* Scale angle theta from [0,1/2] to [0,pi] and apply cosine. */
700  stheta = cos(2.0*KPI*plan->x[2*m+1]);
701  /* Scale angle phi from [-1/2,1/2] to [-pi,pi]. */
702  sphi = 2.0*KPI*plan->x[2*m];
703 
704  /* Initialize result for current node. */
705  f_m = 0.0;
706 
707  /* For n = -N,...,N, evaluate
708  * b_n := \sum_{k=|n|}^N a_k^n P_k^{|n|}(cos theta_m)
709  * using Clenshaw's algorithm.
710  */
711  for (n = -plan->N; n <= plan->N; n++)
712  {
713  /* Get pointer to Fourier coefficients vector for current order n. */
714  a = &(plan->f_hat_intern[NFSFT_INDEX(0,n,plan)]);
715 
716  /* Take absolute value of n. */
717  n_abs = abs(n);
718 
719  /* Get pointers to three-term recurrence coefficients arrays. */
720  alpha = &(wisdom.alpha[ROW(n_abs)]);
721  gamma = &(wisdom.gamma[ROW(n_abs)]);
722 
723  /* Clenshaw's algorithm */
724  it2 = a[plan->N];
725  it1 = a[plan->N-1];
726  for (k = plan->N; k > n_abs + 1; k--)
727  {
728  temp = a[k-2] + it2 * gamma[k];
729  it2 = it1 + it2 * alpha[k] * stheta;
730  it1 = temp;
731  }
732 
733  /* Compute final step if neccesary. */
734  if (n_abs < plan->N)
735  {
736  it2 = it1 + it2 * wisdom.alpha[ROWK(n_abs)+1] * stheta;
737  }
738 
739  /* Compute final result by multiplying the fixed part
740  * gamma_|n| (1-cos^2(theta))^{|n|/2}
741  * for order n and the exponential part
742  * e^{i n phi}
743  * and add the result to f_m.
744  */
745  f_m += it2 * wisdom.gamma[ROW(n_abs)] *
746  pow(1- stheta * stheta, 0.5*n_abs) * cexp(_Complex_I*n*sphi);
747  }
748 
749  /* Write result f_m for current node to array f. */
750  plan->f[m] = f_m;
751  }
752  }
753 }
754 
755 void nfsft_adjoint_direct(nfsft_plan *plan)
756 {
757  int m; /*< The node index */
758  int k; /*< The degree k */
759  int n; /*< The order n */
760  int n_abs; /*< The absolute value of the order n, ie n_abs = |n| */
761  double *alpha; /*< Pointer to current three-term recurrence
762  coefficient alpha_k^n for associated Legendre
763  functions P_k^n */
764  double *gamma; /*< Pointer to current three-term recurrence
765  coefficient beta_k^n for associated Legendre
766  functions P_k^n */
767  double _Complex it1; /*< Auxilliary variable for Clenshaw algorithm */
768  double _Complex it2; /*< Auxilliary variable for Clenshaw algorithm */
769  double _Complex temp; /*< Auxilliary variable for Clenshaw algorithm */
770  double stheta; /*< Current angle theta for Clenshaw algorithm */
771  double sphi; /*< Current angle phi for Clenshaw algorithm */
772 
773 #ifdef MEASURE_TIME
774  plan->MEASURE_TIME_t[0] = 0.0;
775  plan->MEASURE_TIME_t[1] = 0.0;
776  plan->MEASURE_TIME_t[2] = 0.0;
777 #endif
778 
779  if (wisdom.flags & NFSFT_NO_DIRECT_ALGORITHM)
780  {
781  return;
782  }
783 
784  /* Initialise spherical Fourier coefficients array with zeros. */
785  memset(plan->f_hat,0U,plan->N_total*sizeof(double _Complex));
786 
787  /* Distinguish by bandwidth N. */
788  if (plan->N == 0)
789  {
790  /* N == 0 */
791 
792  /* Constant function */
793  for (m = 0; m < plan->M_total; m++)
794  {
795  plan->f_hat[NFSFT_INDEX(0,0,plan)] += plan->f[m];
796  }
797  }
798  else
799  {
800  /* N > 0 */
801 
802 #ifdef _OPENMP
803  /* Traverse all orders n. */
804  #pragma omp parallel for default(shared) private(n,n_abs,alpha,gamma,m,stheta,sphi,it2,it1,k,temp)
805  for (n = -plan->N; n <= plan->N; n++)
806  {
807  /* Take absolute value of n. */
808  n_abs = abs(n);
809 
810  /* Get pointers to three-term recurrence coefficients arrays. */
811  alpha = &(wisdom.alpha[ROW(n_abs)]);
812  gamma = &(wisdom.gamma[ROW(n_abs)]);
813 
814  /* Traverse all nodes. */
815  for (m = 0; m < plan->M_total; m++)
816  {
817  /* Scale angle theta from [0,1/2] to [0,pi] and apply cosine. */
818  stheta = cos(2.0*KPI*plan->x[2*m+1]);
819  /* Scale angle phi from [-1/2,1/2] to [-pi,pi]. */
820  sphi = 2.0*KPI*plan->x[2*m];
821 
822  /* Transposed Clenshaw algorithm */
823 
824  /* Initial step */
825  it1 = plan->f[m] * wisdom.gamma[ROW(n_abs)] *
826  pow(1 - stheta * stheta, 0.5*n_abs) * cexp(-_Complex_I*n*sphi);
827  plan->f_hat[NFSFT_INDEX(n_abs,n,plan)] += it1;
828  it2 = 0.0;
829 
830  if (n_abs < plan->N)
831  {
832  it2 = it1 * wisdom.alpha[ROWK(n_abs)+1] * stheta;
833  plan->f_hat[NFSFT_INDEX(n_abs+1,n,plan)] += it2;
834  }
835 
836  /* Loop for transposed Clenshaw algorithm */
837  for (k = n_abs+2; k <= plan->N; k++)
838  {
839  temp = it2;
840  it2 = alpha[k] * stheta * it2 + gamma[k] * it1;
841  it1 = temp;
842  plan->f_hat[NFSFT_INDEX(k,n,plan)] += it2;
843  }
844  }
845  }
846 #else
847  /* Traverse all nodes. */
848  for (m = 0; m < plan->M_total; m++)
849  {
850  /* Scale angle theta from [0,1/2] to [0,pi] and apply cosine. */
851  stheta = cos(2.0*KPI*plan->x[2*m+1]);
852  /* Scale angle phi from [-1/2,1/2] to [-pi,pi]. */
853  sphi = 2.0*KPI*plan->x[2*m];
854 
855  /* Traverse all orders n. */
856  for (n = -plan->N; n <= plan->N; n++)
857  {
858  /* Take absolute value of n. */
859  n_abs = abs(n);
860 
861  /* Get pointers to three-term recurrence coefficients arrays. */
862  alpha = &(wisdom.alpha[ROW(n_abs)]);
863  gamma = &(wisdom.gamma[ROW(n_abs)]);
864 
865  /* Transposed Clenshaw algorithm */
866 
867  /* Initial step */
868  it1 = plan->f[m] * wisdom.gamma[ROW(n_abs)] *
869  pow(1 - stheta * stheta, 0.5*n_abs) * cexp(-_Complex_I*n*sphi);
870  plan->f_hat[NFSFT_INDEX(n_abs,n,plan)] += it1;
871  it2 = 0.0;
872 
873  if (n_abs < plan->N)
874  {
875  it2 = it1 * wisdom.alpha[ROWK(n_abs)+1] * stheta;
876  plan->f_hat[NFSFT_INDEX(n_abs+1,n,plan)] += it2;
877  }
878 
879  /* Loop for transposed Clenshaw algorithm */
880  for (k = n_abs+2; k <= plan->N; k++)
881  {
882  temp = it2;
883  it2 = alpha[k] * stheta * it2 + gamma[k] * it1;
884  it1 = temp;
885  plan->f_hat[NFSFT_INDEX(k,n,plan)] += it2;
886  }
887  }
888  }
889 #endif
890  }
891 
892  /* Check, if we compute with L^2-normalized spherical harmonics. If so,
893  * multiply spherical Fourier coefficients with corresponding normalization
894  * weight. */
895  if (plan->flags & NFSFT_NORMALIZED)
896  {
897  /* Traverse Fourier coefficients array. */
898 #ifdef _OPENMP
899  #pragma omp parallel for default(shared) private(k,n)
900 #endif
901  for (k = 0; k <= plan->N; k++)
902  {
903  for (n = -k; n <= k; n++)
904  {
905  /* Multiply with normalization weight. */
906  plan->f_hat[NFSFT_INDEX(k,n,plan)] *=
907  sqrt((2*k+1)/(4.0*KPI));
908  }
909  }
910  }
911 
912  /* Set unused coefficients to zero. */
913  if (plan->flags & NFSFT_ZERO_F_HAT)
914  {
915  for (n = -plan->N; n <= plan->N+1; n++)
916  {
917  memset(&plan->f_hat[NFSFT_INDEX(-plan->N-1,n,plan)],0U,
918  (plan->N+1+abs(n))*sizeof(double _Complex));
919  }
920  }
921 }
922 
923 void nfsft_trafo(nfsft_plan *plan)
924 {
925  int k; /*< The degree k */
926  int n; /*< The order n */
927 #ifdef MEASURE_TIME
928  ticks t0, t1;
929 #endif
930  #ifdef DEBUG
931  double t, t_pre, t_nfft, t_fpt, t_c2e, t_norm;
932  t_pre = 0.0;
933  t_norm = 0.0;
934  t_fpt = 0.0;
935  t_c2e = 0.0;
936  t_nfft = 0.0;
937  #endif
938 
939 #ifdef MEASURE_TIME
940  plan->MEASURE_TIME_t[0] = 0.0;
941  plan->MEASURE_TIME_t[1] = 0.0;
942  plan->MEASURE_TIME_t[2] = 0.0;
943 #endif
944 
945  if (wisdom.flags & NFSFT_NO_FAST_ALGORITHM)
946  {
947  return;
948  }
949 
950  /* Check, if precomputation was done and that the bandwidth N is not too
951  * big.
952  */
953  if (wisdom.initialized == 0 || plan->N > wisdom.N_MAX)
954  {
955  return;
956  }
957 
958  /* Check, if slow transformation should be used due to small bandwidth. */
959  if (plan->N < NFSFT_BREAK_EVEN)
960  {
961  /* Use NDSFT. */
962  nfsft_trafo_direct(plan);
963  }
964 
965  /* Check for correct value of the bandwidth N. */
966  else if (plan->N <= wisdom.N_MAX)
967  {
968  /* Copy spherical Fourier coefficients, if necessary. */
969  if (plan->flags & NFSFT_PRESERVE_F_HAT)
970  {
971  memcpy(plan->f_hat_intern,plan->f_hat,plan->N_total*
972  sizeof(double _Complex));
973  }
974  else
975  {
976  plan->f_hat_intern = plan->f_hat;
977  }
978 
979  /* Propagate pointer values to the internal NFFT plan to assure
980  * consistency. Pointers may have been modified externally.
981  */
982  plan->plan_nfft.x = plan->x;
983  plan->plan_nfft.f = plan->f;
984  plan->plan_nfft.f_hat = plan->f_hat_intern;
985 
986  /* Check, if we compute with L^2-normalized spherical harmonics. If so,
987  * multiply spherical Fourier coefficients with corresponding normalization
988  * weight. */
989  if (plan->flags & NFSFT_NORMALIZED)
990  {
991  /* Traverse Fourier coefficients array. */
992 #ifdef _OPENMP
993  #pragma omp parallel for default(shared) private(k,n)
994 #endif
995  for (k = 0; k <= plan->N; k++)
996  {
997  for (n = -k; n <= k; n++)
998  {
999  /* Multiply with normalization weight. */
1000  plan->f_hat_intern[NFSFT_INDEX(k,n,plan)] *=
1001  sqrt((2*k+1)/(4.0*KPI));
1002  }
1003  }
1004  }
1005 
1006 #ifdef MEASURE_TIME
1007  t0 = getticks();
1008 #endif
1009  /* Check, which polynomial transform algorithm should be used. */
1010  if (plan->flags & NFSFT_USE_DPT)
1011  {
1012 #ifdef _OPENMP
1013  #pragma omp parallel for default(shared) private(n) num_threads(wisdom.nthreads)
1014  for (n = -plan->N; n <= plan->N; n++)
1015  fpt_trafo_direct(wisdom.set_threads[omp_get_thread_num()],abs(n),
1016  &plan->f_hat_intern[NFSFT_INDEX(abs(n),n,plan)],
1017  &plan->f_hat_intern[NFSFT_INDEX(0,n,plan)],
1018  plan->N,0U);
1019 #else
1020  /* Use direct discrete polynomial transform DPT. */
1021  for (n = -plan->N; n <= plan->N; n++)
1022  {
1023  //fprintf(stderr,"nfsft_trafo: n = %d\n",n);
1024  //fflush(stderr);
1025  fpt_trafo_direct(wisdom.set,abs(n),
1026  &plan->f_hat_intern[NFSFT_INDEX(abs(n),n,plan)],
1027  &plan->f_hat_intern[NFSFT_INDEX(0,n,plan)],
1028  plan->N,0U);
1029  }
1030 #endif
1031  }
1032  else
1033  {
1034 #ifdef _OPENMP
1035  #pragma omp parallel for default(shared) private(n) num_threads(wisdom.nthreads)
1036  for (n = -plan->N; n <= plan->N; n++)
1037  fpt_trafo(wisdom.set_threads[omp_get_thread_num()],abs(n),
1038  &plan->f_hat_intern[NFSFT_INDEX(abs(n),n,plan)],
1039  &plan->f_hat_intern[NFSFT_INDEX(0,n,plan)],
1040  plan->N,0U);
1041 #else
1042  /* Use fast polynomial transform FPT. */
1043  for (n = -plan->N; n <= plan->N; n++)
1044  {
1045  //fprintf(stderr,"nfsft_trafo: n = %d\n",n);
1046  //fflush(stderr);
1047  fpt_trafo(wisdom.set,abs(n),
1048  &plan->f_hat_intern[NFSFT_INDEX(abs(n),n,plan)],
1049  &plan->f_hat_intern[NFSFT_INDEX(0,n,plan)],
1050  plan->N,0U);
1051  }
1052 #endif
1053  }
1054 #ifdef MEASURE_TIME
1055  t1 = getticks();
1056  plan->MEASURE_TIME_t[0] = nfft_elapsed_seconds(t1,t0);
1057 #endif
1058 
1059 #ifdef MEASURE_TIME
1060  t0 = getticks();
1061 #endif
1062  /* Convert Chebyshev coefficients to Fourier coefficients. */
1063  c2e(plan);
1064 #ifdef MEASURE_TIME
1065  t1 = getticks();
1066  plan->MEASURE_TIME_t[1] = nfft_elapsed_seconds(t1,t0);
1067 #endif
1068 
1069 #ifdef MEASURE_TIME
1070  t0 = getticks();
1071 #endif
1072  /* Check, which nonequispaced discrete Fourier transform algorithm should
1073  * be used.
1074  */
1075  if (plan->flags & NFSFT_USE_NDFT)
1076  {
1077  /* Use NDFT. */
1078  nfft_trafo_direct(&plan->plan_nfft);
1079  }
1080  else
1081  {
1082  /* Use NFFT. */
1083  //fprintf(stderr,"nfsft_adjoint: nfft_trafo\n");
1084  nfft_trafo_2d(&plan->plan_nfft);
1085  }
1086 #ifdef MEASURE_TIME
1087  t1 = getticks();
1088  plan->MEASURE_TIME_t[2] = nfft_elapsed_seconds(t1,t0);
1089 #endif
1090  }
1091 }
1092 
1093 void nfsft_adjoint(nfsft_plan *plan)
1094 {
1095  int k; /*< The degree k */
1096  int n; /*< The order n */
1097 #ifdef MEASURE_TIME
1098  ticks t0, t1;
1099 #endif
1100 
1101 #ifdef MEASURE_TIME
1102  plan->MEASURE_TIME_t[0] = 0.0;
1103  plan->MEASURE_TIME_t[1] = 0.0;
1104  plan->MEASURE_TIME_t[2] = 0.0;
1105 #endif
1106 
1107  if (wisdom.flags & NFSFT_NO_FAST_ALGORITHM)
1108  {
1109  return;
1110  }
1111 
1112  /* Check, if precomputation was done and that the bandwidth N is not too
1113  * big.
1114  */
1115  if (wisdom.initialized == 0 || plan->N > wisdom.N_MAX)
1116  {
1117  return;
1118  }
1119 
1120  /* Check, if slow transformation should be used due to small bandwidth. */
1121  if (plan->N < NFSFT_BREAK_EVEN)
1122  {
1123  /* Use adjoint NDSFT. */
1124  nfsft_adjoint_direct(plan);
1125  }
1126  /* Check for correct value of the bandwidth N. */
1127  else if (plan->N <= wisdom.N_MAX)
1128  {
1129  //fprintf(stderr,"nfsft_adjoint: Starting\n");
1130  //fflush(stderr);
1131  /* Propagate pointer values to the internal NFFT plan to assure
1132  * consistency. Pointers may have been modified externally.
1133  */
1134  plan->plan_nfft.x = plan->x;
1135  plan->plan_nfft.f = plan->f;
1136  plan->plan_nfft.f_hat = plan->f_hat;
1137 
1138 #ifdef MEASURE_TIME
1139  t0 = getticks();
1140 #endif
1141  /* Check, which adjoint nonequispaced discrete Fourier transform algorithm
1142  * should be used.
1143  */
1144  if (plan->flags & NFSFT_USE_NDFT)
1145  {
1146  //fprintf(stderr,"nfsft_adjoint: Executing nfft_adjoint_direct\n");
1147  //fflush(stderr);
1148  /* Use adjoint NDFT. */
1149  nfft_adjoint_direct(&plan->plan_nfft);
1150  }
1151  else
1152  {
1153  //fprintf(stderr,"nfsft_adjoint: Executing nfft_adjoint\n");
1154  //fflush(stderr);
1155  //fprintf(stderr,"nfsft_adjoint: nfft_adjoint\n");
1156  /* Use adjoint NFFT. */
1157  nfft_adjoint_2d(&plan->plan_nfft);
1158  }
1159 #ifdef MEASURE_TIME
1160  t1 = getticks();
1161  plan->MEASURE_TIME_t[2] = nfft_elapsed_seconds(t1,t0);
1162 #endif
1163 
1164  //fprintf(stderr,"nfsft_adjoint: Executing c2e_transposed\n");
1165  //fflush(stderr);
1166 #ifdef MEASURE_TIME
1167  t0 = getticks();
1168 #endif
1169  /* Convert Fourier coefficients to Chebyshev coefficients. */
1170  c2e_transposed(plan);
1171 #ifdef MEASURE_TIME
1172  t1 = getticks();
1173  plan->MEASURE_TIME_t[1] = nfft_elapsed_seconds(t1,t0);
1174 #endif
1175 
1176 #ifdef MEASURE_TIME
1177  t0 = getticks();
1178 #endif
1179  /* Check, which transposed polynomial transform algorithm should be used */
1180  if (plan->flags & NFSFT_USE_DPT)
1181  {
1182 #ifdef _OPENMP
1183  #pragma omp parallel for default(shared) private(n) num_threads(wisdom.nthreads)
1184  for (n = -plan->N; n <= plan->N; n++)
1185  fpt_transposed_direct(wisdom.set_threads[omp_get_thread_num()],abs(n),
1186  &plan->f_hat[NFSFT_INDEX(abs(n),n,plan)],
1187  &plan->f_hat[NFSFT_INDEX(0,n,plan)],
1188  plan->N,0U);
1189 #else
1190  /* Use transposed DPT. */
1191  for (n = -plan->N; n <= plan->N; n++)
1192  {
1193  //fprintf(stderr,"nfsft_adjoint: Executing dpt_transposed\n");
1194  //fflush(stderr);
1195  fpt_transposed_direct(wisdom.set,abs(n),
1196  &plan->f_hat[NFSFT_INDEX(abs(n),n,plan)],
1197  &plan->f_hat[NFSFT_INDEX(0,n,plan)],
1198  plan->N,0U);
1199  }
1200 #endif
1201  }
1202  else
1203  {
1204 #ifdef _OPENMP
1205  #pragma omp parallel for default(shared) private(n) num_threads(wisdom.nthreads)
1206  for (n = -plan->N; n <= plan->N; n++)
1207  fpt_transposed(wisdom.set_threads[omp_get_thread_num()],abs(n),
1208  &plan->f_hat[NFSFT_INDEX(abs(n),n,plan)],
1209  &plan->f_hat[NFSFT_INDEX(0,n,plan)],
1210  plan->N,0U);
1211 #else
1212  //fprintf(stderr,"nfsft_adjoint: fpt_transposed\n");
1213  /* Use transposed FPT. */
1214  for (n = -plan->N; n <= plan->N; n++)
1215  {
1216  //fprintf(stderr,"nfsft_adjoint: Executing fpt_transposed\n");
1217  //fflush(stderr);
1218  fpt_transposed(wisdom.set,abs(n),
1219  &plan->f_hat[NFSFT_INDEX(abs(n),n,plan)],
1220  &plan->f_hat[NFSFT_INDEX(0,n,plan)],
1221  plan->N,0U);
1222  }
1223 #endif
1224  }
1225 #ifdef MEASURE_TIME
1226  t1 = getticks();
1227  plan->MEASURE_TIME_t[0] = nfft_elapsed_seconds(t1,t0);
1228 #endif
1229 
1230  /* Check, if we compute with L^2-normalized spherical harmonics. If so,
1231  * multiply spherical Fourier coefficients with corresponding normalization
1232  * weight. */
1233  if (plan->flags & NFSFT_NORMALIZED)
1234  {
1235  //fprintf(stderr,"nfsft_adjoint: Normalizing\n");
1236  //fflush(stderr);
1237  /* Traverse Fourier coefficients array. */
1238 #ifdef _OPENMP
1239  #pragma omp parallel for default(shared) private(k,n)
1240 #endif
1241  for (k = 0; k <= plan->N; k++)
1242  {
1243  for (n = -k; n <= k; n++)
1244  {
1245  /* Multiply with normalization weight. */
1246  plan->f_hat[NFSFT_INDEX(k,n,plan)] *=
1247  sqrt((2*k+1)/(4.0*KPI));
1248  }
1249  }
1250  }
1251 
1252  /* Set unused coefficients to zero. */
1253  if (plan->flags & NFSFT_ZERO_F_HAT)
1254  {
1255  //fprintf(stderr,"nfsft_adjoint: Setting to zero\n");
1256  //fflush(stderr);
1257  for (n = -plan->N; n <= plan->N+1; n++)
1258  {
1259  memset(&plan->f_hat[NFSFT_INDEX(-plan->N-1,n,plan)],0U,
1260  (plan->N+1+abs(n))*sizeof(double _Complex));
1261  }
1262  }
1263  //fprintf(stderr,"nfsft_adjoint: Finished\n");
1264  //fflush(stderr);
1265  }
1266 }
1267 
1268 void nfsft_precompute_x(nfsft_plan *plan)
1269 {
1270  /* Pass angle array to NFFT plan. */
1271  plan->plan_nfft.x = plan->x;
1272 
1273  /* Precompute. */
1274  if (plan->plan_nfft.flags & PRE_ONE_PSI)
1275  nfft_precompute_one_psi(&plan->plan_nfft);
1276 }
1277 /* \} */
static void c2e(nfsft_plan *plan)
Converts coefficients with , from a linear combination of Chebyshev polynomials to coefficients m...
Definition: nfsft.c:110
double * gamma
Precomputed recursion coefficients /f$^n/f$ for /f$k = 0,/ldots, N_{{max}}; n=-k,/ldots,k/f$ of associated Legendre-functions /f$P_k^n/f$.
double * x
the nodes for ,
Definition: nfft3.h:576
void(* mv_trafo)(void *)
Transform.
Definition: nfft3.h:576
fftw_complex * f_hat
Fourier coefficients.
Definition: nfft3.h:194
Wisdom structure.
fftw_complex * f_hat
Fourier coefficients.
Definition: nfft3.h:576
R nfft_elapsed_seconds(ticks t1, ticks t0)
Return number of elapsed seconds between two time points.
int N
the bandwidth
Definition: nfft3.h:576
data structure for an NFSFT (nonequispaced fast spherical Fourier transform) plan with double precisi...
Definition: nfft3.h:576
fpt_set set
Structure for discrete polynomial transform (DPT)
#define NFSFT_BREAK_EVEN
The break-even bandwidth .
Definition: nfsft.c:78
int N_MAX
Stores precomputation flags.
bool initialized
Indicates wether the structure has been initialized.
void alpha_al_all(R *alpha, const int N)
Compute three-term-recurrence coefficients of associated Legendre functions for .
Definition: legendre.c:91
unsigned int flags
the planner flags
Definition: nfft3.h:576
Holds data for a set of cascade summations.
Definition: fpt.c:96
fftw_complex * f
Samples.
Definition: nfft3.h:194
#define NFSFT_DEFAULT_NFFT_CUTOFF
The default NFFT cutoff parameter.
Definition: nfsft.c:64
void nfft_free(void *p)
nfft_plan plan_nfft
the internal NFFT plan
Definition: nfft3.h:576
NFFT_INT M_total
Total number of samples.
Definition: nfft3.h:576
double MEASURE_TIME_t[3]
Measured time for each step if MEASURE_TIME is set.
Definition: nfft3.h:576
#define X(name)
Include header for C99 complex datatype.
Definition: fastsum.h:53
void * nfft_malloc(size_t n)
NFFT_INT N_total
Total number of Fourier coefficients.
Definition: nfft3.h:576
double * alpha
Precomputed recursion coefficients /f$^n/f$ for /f$k = 0,/ldots, N_{{max}}; n=-k,/ldots,k/f$ of associated Legendre-functions /f$P_k^n/f$.
fftw_complex * f
Samples.
Definition: nfft3.h:576
void(* mv_adjoint)(void *)
Adjoint transform.
Definition: nfft3.h:576
void beta_al_all(R *alpha, const int N)
Compute three-term-recurrence coefficients of associated Legendre functions for .
Definition: legendre.c:100
double * beta
Precomputed recursion coefficients /f$^n/f$ for /f$k = 0,/ldots, N_{{max}}; n=-k,/ldots,k/f$ of associated Legendre-functions /f$P_k^n/f$.
static void c2e_transposed(nfsft_plan *plan)
Transposed version of the function c2e.
Definition: nfsft.c:188
void gamma_al_all(R *alpha, const int N)
Compute three-term-recurrence coefficients of associated Legendre functions for .
Definition: legendre.c:109
double * x
Nodes in time/spatial domain, size is doubles.
Definition: nfft3.h:194
static struct nfsft_wisdom wisdom
The global wisdom structure for precomputed data.
Definition: nfsft.c:86
unsigned flags
Flags for precomputation, (de)allocation, and FFTW usage, default setting is PRE_PHI_HUT | PRE_PSI | ...
Definition: nfft3.h:194
fftw_complex * f_hat_intern
Internally used pointer to spherical Fourier coefficients.
Definition: nfft3.h:576
int T_MAX
The logarithm /f$t = N_{{max}}/f$ of the maximum bandwidth.