32 #define DEFAULT_NFFT_CUTOFF 6 33 #define FPT_THRESHOLD 1000 35 static fpt_set SO3_fpt_init(
int l,
unsigned int flags,
int kappa);
39 nfsoft_init_advanced(plan, N, M, NFSOFT_MALLOC_X | NFSOFT_MALLOC_F
40 | NFSOFT_MALLOC_F_HAT);
43 void nfsoft_init_advanced(
nfsoft_plan *plan,
int N,
int M,
44 unsigned int nfsoft_flags)
46 nfsoft_init_guru(plan, N, M, nfsoft_flags, PRE_PHI_HUT | PRE_PSI | MALLOC_X | NFFT_OMP_BLOCKWISE_ADJOINT
47 | MALLOC_F_HAT | MALLOC_F | FFTW_INIT | FFT_OUT_OF_PLACE,
48 DEFAULT_NFFT_CUTOFF, FPT_THRESHOLD);
51 void nfsoft_init_guru(
nfsoft_plan *plan,
int B,
int M,
52 unsigned int nfsoft_flags,
unsigned int nfft_flags,
int nfft_cutoff,
66 nfft_init_guru(&plan->
p_nfft, 3, N, M, n, nfft_cutoff, nfft_flags,
67 FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
69 if ((plan->
p_nfft).flags & PRE_LIN_PSI)
71 nfft_precompute_lin_psi(&(plan->
p_nfft));
76 plan->fpt_kappa = fpt_kappa;
77 plan->
flags = nfsoft_flags;
79 if (plan->
flags & NFSOFT_MALLOC_F_HAT)
82 if (plan->
f_hat == NULL ) printf(
"Allocation failed!\n");
85 if (plan->
flags & NFSOFT_MALLOC_X)
88 if (plan->
x == NULL ) printf(
"Allocation failed!\n");
90 if (plan->
flags & NFSOFT_MALLOC_F)
93 if (plan->
f == NULL ) printf(
"Allocation failed!\n");
100 if (plan->
wig_coeffs == NULL ) printf(
"Allocation failed!\n");
101 if (plan->
cheby == NULL ) printf(
"Allocation failed!\n");
102 if (plan->
aux == NULL ) printf(
"Allocation failed!\n");
104 plan->
mv_trafo = (void (*) (
void* ))nfsoft_trafo;
105 plan->
mv_adjoint = (void (*) (
void* ))nfsoft_adjoint;
120 my_plan->
cheby[0]=0.0;
122 for (j=1;j<my_plan->
N_total+1;j++)
131 aux[j]=my_plan->
cheby[j];
138 my_plan->
cheby[0]=(C) (-1.)/(2.*_Complex_I) * aux[1];
141 my_plan->
cheby[j]=(1./(2.*_Complex_I)*(aux[j+1]-aux[j-1]));
150 static fpt_set SO3_fpt_init(
int l,
unsigned int flags,
int kappa)
153 int N, t, k_start, k_end, k, m;
158 if (flags & NFSOFT_USE_DPT)
165 t = (int) log2(
X(next_power_of_2)(N));
174 N =
X(next_power_of_2)(l);
185 if (flags & NFSOFT_NO_STABILIZATION)
187 set = fpt_init((2* N + 1) * (2* N + 1), t, 0U | FPT_NO_STABILIZATION);
191 set = fpt_init((2* N + 1) * (2* N + 1), t, 0U);
194 for (k = -N; k <= N; k++)
195 for (m = -N; m <= N; m++)
198 k_start = (ABS(k) >= ABS(m)) ? ABS(k) : ABS(m);
205 fpt_precompute(
set, glo, alpha, beta, gamma, k_start, kappa);
219 static fpt_set SO3_single_fpt_init(
int l,
int k,
int m,
unsigned int flags,
int kappa)
221 int N, t, k_start, k_end;
226 if (flags & NFSOFT_USE_DPT)
233 t = (int) log2(
X(next_power_of_2)(N));
242 N =
X(next_power_of_2)(l);
254 unsigned int fptflags = 0U
255 | IF(flags & NFSOFT_USE_DPT,FPT_NO_FAST_ALGORITHM,IF(t > 1,FPT_NO_DIRECT_ALGORITHM,0U))
256 | IF(flags & NFSOFT_NO_STABILIZATION,FPT_NO_STABILIZATION,0U);
257 set = fpt_init(1, t, fptflags);
261 k_start = (ABS(k) >= ABS(m)) ? ABS(k) : ABS(m);
274 fpt_precompute(
set, 0, alpha, beta, gamma, k_start, kappa);
286 void SO3_fpt(C *coeffs,
fpt_set set,
int l,
int k,
int m,
unsigned int flags)
295 int k_start, k_end, j;
296 int function_values = 0;
299 if (flags & NFSOFT_USE_DPT)
310 N =
X(next_power_of_2)(l);
314 k_start = (ABS(k) >= ABS(m)) ? ABS(k) : ABS(m);
316 trafo_nr = (N + k) * (2* N + 1) + (m + N);
321 for (j = 0; j <= k_end; j++)
325 for (j = 0; j <= l - k_start; j++)
327 x[j + k_start] = coeffs[j];
329 for (j = l - k_start + 1; j <= k_end - k_start; j++)
331 x[j + k_start] = K(0.0);
337 if (flags & NFSOFT_USE_DPT)
339 fpt_trafo_direct(
set, trafo_nr, &x[k_start], y, k_end, 0U
340 | (function_values ? FPT_FUNCTION_VALUES : 0U));
344 fpt_trafo(
set, trafo_nr, &x[k_start], y, k_end, 0U
345 | (function_values ? FPT_FUNCTION_VALUES : 0U));
349 for (j = 0; j <= l; j++)
362 void SO3_fpt_transposed(C *coeffs,
fpt_set set,
int l,
int k,
int m,
365 int N, k_start, k_end, j;
367 int function_values = 0;
375 if (flags & NFSOFT_USE_DPT)
386 N =
X(next_power_of_2)(l);
390 k_start = (ABS(k) >= ABS(m)) ? ABS(k) : ABS(m);
392 trafo_nr = (N + k) * (2* N + 1) + (m + N);
399 for (j = 0; j <= l; j++)
403 for (j = l + 1; j <= k_end; j++)
408 if (flags & NFSOFT_USE_DPT)
410 fpt_transposed_direct(
set, trafo_nr, &x[k_start], y, k_end, 0U
411 | (function_values ? FPT_FUNCTION_VALUES : 0U));
415 fpt_transposed(
set, trafo_nr, &x[k_start], y, k_end, 0U
416 | (function_values ? FPT_FUNCTION_VALUES : 0U));
419 for (j = 0; j <= l; j++)
439 for (j = 0; j < M; j++)
441 plan3D->
p_nfft.
x[3* j ] = plan3D->
x[3* j + 2];
442 plan3D->
p_nfft.
x[3* j + 1] = plan3D->
x[3* j ];
443 plan3D->
p_nfft.
x[3* j + 2] = plan3D->
x[3* j + 1];
451 if ((plan3D->
p_nfft).flags & FG_PSI)
453 nfft_precompute_one_psi(&(plan3D->
p_nfft));
455 if ((plan3D->
p_nfft).flags & PRE_PSI)
457 nfft_precompute_one_psi(&(plan3D->
p_nfft));
464 int i, j, m, k, max, glo1, glo2;
476 for (j = 0; j < M; j++)
477 plan3D->
f[j] = plan3D->
f_hat[0];
484 for (k = -N; k <= N; k++)
486 for (m = -N; m <= N; m++)
489 max = (ABS(m) > ABS(k) ? ABS(m) : ABS(k));
491 for (j = 0; j <= N - max; j++)
495 if ((plan3D->
flags & NFSOFT_NORMALIZED))
498 * SQRT(0.5 * (2. * (max + j) + 1.));
501 if ((plan3D->
flags & NFSOFT_REPRESENT))
503 if ((k < 0) && (k % 2))
507 if ((m < 0) && (m % 2))
518 for (j = N - max + 1; j <
X(next_power_of_2)(N) + 1; j++)
523 c2e(plan3D, ABS((k + m) % 2));
525 for (i = 1; i <= 2* plan3D ->
N_total + 2; i++)
527 plan3D->
p_nfft.
f_hat[NFSOFT_INDEX(k, m, i - N - 1, N) - 1]
528 = plan3D->
cheby[i - 1];
536 if (plan3D->
flags & NFSOFT_USE_NDFT)
538 nfft_trafo_direct(&(plan3D->
p_nfft));
542 nfft_trafo(&(plan3D->
p_nfft));
545 for (j = 0; j < plan3D->
M_total; j++)
563 my_plan->
aux[0]= 1/(2*_Complex_I)*my_plan->
cheby[1];
567 my_plan->
aux[j]=1/(2*_Complex_I)*(my_plan->
cheby[j+1]-my_plan->
cheby[j-1]);
569 my_plan->
aux[N-1]=1/(2*_Complex_I)*(-my_plan->
cheby[j-1]);
580 for(j=1;j<=my_plan->
N_total;j++)
593 int i, j, m, k, max, glo1, glo2;
606 for (j = 0; j < M; j++)
607 plan3D->
f_hat[0] += plan3D->
f[j];
611 for (j = 0; j < M; j++)
616 if (plan3D->
flags & NFSOFT_USE_NDFT)
618 nfft_adjoint_direct(&(plan3D->
p_nfft));
622 nfft_adjoint(&(plan3D->
p_nfft));
629 for (k = -N; k <= N; k++)
631 for (m = -N; m <= N; m++)
634 max = (ABS(m) > ABS(k) ? ABS(m) : ABS(k));
636 for (i = 1; i < 2* plan3D ->
N_total + 3; i++)
644 e2c(plan3D, ABS((k + m) % 2));
653 for (j = max; j <= N; j++)
655 if ((plan3D->
flags & NFSOFT_REPRESENT))
657 if ((k < 0) && (k % 2))
661 if ((m < 0) && (m % 2))
671 if ((plan3D->
flags & NFSOFT_NORMALIZED))
673 plan3D->
f_hat[glo1] = plan3D->
f_hat[glo1] * (1 / (2. * KPI)) * SQRT(
674 0.5 * (2. * (j) + 1.));
687 nfft_finalize(&plan->
p_nfft);
695 if (plan->
flags & NFSOFT_MALLOC_F_HAT)
702 if (plan->
flags & NFSOFT_MALLOC_F)
709 if (plan->
flags & NFSOFT_MALLOC_X)
716 int posN(
int n,
int m,
int B)
721 pos = posN(n - 1, m, B) + B + 1 - MAX(ABS(m), ABS(n - 1));
static void c2e(nfsft_plan *plan)
Converts coefficients with , from a linear combination of Chebyshev polynomials to coefficients m...
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$.
fftw_complex * f_hat
Fourier coefficients.
fftw_complex * aux
used when converting Chebychev to Fourier coeffcients
fftw_complex * f_hat
Fourier coefficients.
void SO3_gamma_row(double *gamma, int N, int m, int n)
Compute three-term-recurrence coefficients of Wigner-d functions for all degrees ...
fftw_complex * wig_coeffs
contains a set of SO(3) Fourier coefficients for fixed orders m and n
void SO3_beta_row(double *beta, int N, int m, int n)
Compute three-term-recurrence coefficients of Wigner-d functions for all degrees ...
NFFT_INT M_total
Total number of samples.
Holds data for a set of cascade summations.
void(* mv_adjoint)(void *)
Adjoint transform.
void SO3_alpha_row(double *alpha, int N, int m, int n)
Compute three-term-recurrence coefficients of Wigner-d functions for all degrees ...
NFFT_INT N_total
Total number of Fourier coefficients.
fpt_set internal_fpt_set
the internal FPT plan
fftw_complex * cheby
contains a set of Chebychev coefficients for fixed orders m and n
NFFT_INT N_total
Total number of Fourier coefficients.
NFFT_INT M_total
Total number of samples.
#define X(name)
Include header for C99 complex datatype.
void * nfft_malloc(size_t n)
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$.
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$.
Header file for functions related to Wigner-d/D functions.
void(* mv_trafo)(void *)
Transform.
double * x
Nodes in time/spatial domain, size is doubles.
nfft_plan p_nfft
the internal NFFT plan
unsigned int flags
the planner flags