34 #define DEFAULT_NFFT_CUTOFF 6 35 #define FPT_THRESHOLD 1000 37 static fpt_set SO3_fpt_init(
int l,
unsigned int flags,
int kappa);
41 nfsoft_init_advanced(plan, N, M, NFSOFT_MALLOC_X | NFSOFT_MALLOC_F
42 | NFSOFT_MALLOC_F_HAT);
45 void nfsoft_init_advanced(
nfsoft_plan *plan,
int N,
int M,
46 unsigned int nfsoft_flags)
48 nfsoft_init_guru(plan, N, M, nfsoft_flags, PRE_PHI_HUT | PRE_PSI | MALLOC_X | NFFT_OMP_BLOCKWISE_ADJOINT
49 | MALLOC_F_HAT | MALLOC_F | FFTW_INIT | FFT_OUT_OF_PLACE,
50 DEFAULT_NFFT_CUTOFF, FPT_THRESHOLD);
53 void nfsoft_init_guru(
nfsoft_plan *plan,
int B,
int M,
54 unsigned int nfsoft_flags,
unsigned int nfft_flags,
int nfft_cutoff,
68 nfft_init_guru(&plan->
p_nfft, 3, N, M, n, nfft_cutoff, nfft_flags,
69 FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
71 if ((plan->
p_nfft).flags & PRE_LIN_PSI)
73 nfft_precompute_lin_psi(&(plan->
p_nfft));
78 plan->fpt_kappa = fpt_kappa;
79 plan->
flags = nfsoft_flags;
81 if (plan->
flags & NFSOFT_MALLOC_F_HAT)
84 if (plan->
f_hat == NULL ) printf(
"Allocation failed!\n");
87 if (plan->
flags & NFSOFT_MALLOC_X)
90 if (plan->
x == NULL ) printf(
"Allocation failed!\n");
92 if (plan->
flags & NFSOFT_MALLOC_F)
95 if (plan->
f == NULL ) printf(
"Allocation failed!\n");
102 if (plan->
wig_coeffs == NULL ) printf(
"Allocation failed!\n");
103 if (plan->
cheby == NULL ) printf(
"Allocation failed!\n");
104 if (plan->
aux == NULL ) printf(
"Allocation failed!\n");
106 plan->
mv_trafo = (void (*) (
void* ))nfsoft_trafo;
107 plan->
mv_adjoint = (void (*) (
void* ))nfsoft_adjoint;
122 my_plan->
cheby[0]=0.0;
124 for (j=1;j<my_plan->
N_total+1;j++)
133 aux[j]=my_plan->
cheby[j];
140 my_plan->
cheby[0]=(C) (-1.)/(2.*_Complex_I) * aux[1];
143 my_plan->
cheby[j]=(1./(2.*_Complex_I)*(aux[j+1]-aux[j-1]));
152 static fpt_set SO3_fpt_init(
int l,
unsigned int flags,
int kappa)
155 int N, t, k_start, k_end, k, m;
160 if (flags & NFSOFT_USE_DPT)
167 t = (int) log2(
X(next_power_of_2)(N));
176 N =
X(next_power_of_2)(l);
187 if (flags & NFSOFT_NO_STABILIZATION)
189 set = fpt_init((2* N + 1) * (2* N + 1), t, 0U | FPT_NO_STABILIZATION);
193 set = fpt_init((2* N + 1) * (2* N + 1), t, 0U);
196 for (k = -N; k <= N; k++)
197 for (m = -N; m <= N; m++)
200 k_start = (ABS(k) >= ABS(m)) ? ABS(k) : ABS(m);
207 fpt_precompute(
set, glo, alpha, beta, gamma, k_start, kappa);
221 static fpt_set SO3_single_fpt_init(
int l,
int k,
int m,
unsigned int flags,
int kappa)
223 int N, t, k_start, k_end;
228 if (flags & NFSOFT_USE_DPT)
235 t = (int) log2(
X(next_power_of_2)(N));
244 N =
X(next_power_of_2)(l);
256 unsigned int fptflags = 0U
257 | IF(flags & NFSOFT_USE_DPT,FPT_NO_FAST_ALGORITHM,IF(t > 1,FPT_NO_DIRECT_ALGORITHM,0U))
258 | IF(flags & NFSOFT_NO_STABILIZATION,FPT_NO_STABILIZATION,0U);
259 set = fpt_init(1, t, fptflags);
263 k_start = (ABS(k) >= ABS(m)) ? ABS(k) : ABS(m);
276 fpt_precompute(
set, 0, alpha, beta, gamma, k_start, kappa);
288 void SO3_fpt(C *coeffs,
fpt_set set,
int l,
int k,
int m,
unsigned int flags)
297 int k_start, k_end, j;
298 int function_values = 0;
301 if (flags & NFSOFT_USE_DPT)
312 N =
X(next_power_of_2)(l);
316 k_start = (ABS(k) >= ABS(m)) ? ABS(k) : ABS(m);
318 trafo_nr = (N + k) * (2* N + 1) + (m + N);
323 for (j = 0; j <= k_end; j++)
327 for (j = 0; j <= l - k_start; j++)
329 x[j + k_start] = coeffs[j];
331 for (j = l - k_start + 1; j <= k_end - k_start; j++)
333 x[j + k_start] = K(0.0);
339 if (flags & NFSOFT_USE_DPT)
341 fpt_trafo_direct(
set, trafo_nr, &x[k_start], y, k_end, 0U
342 | (function_values ? FPT_FUNCTION_VALUES : 0U));
346 fpt_trafo(
set, trafo_nr, &x[k_start], y, k_end, 0U
347 | (function_values ? FPT_FUNCTION_VALUES : 0U));
351 for (j = 0; j <= l; j++)
364 void SO3_fpt_transposed(C *coeffs,
fpt_set set,
int l,
int k,
int m,
367 int N, k_start, k_end, j;
369 int function_values = 0;
377 if (flags & NFSOFT_USE_DPT)
388 N =
X(next_power_of_2)(l);
392 k_start = (ABS(k) >= ABS(m)) ? ABS(k) : ABS(m);
394 trafo_nr = (N + k) * (2* N + 1) + (m + N);
401 for (j = 0; j <= l; j++)
405 for (j = l + 1; j <= k_end; j++)
410 if (flags & NFSOFT_USE_DPT)
412 fpt_transposed_direct(
set, trafo_nr, &x[k_start], y, k_end, 0U
413 | (function_values ? FPT_FUNCTION_VALUES : 0U));
417 fpt_transposed(
set, trafo_nr, &x[k_start], y, k_end, 0U
418 | (function_values ? FPT_FUNCTION_VALUES : 0U));
421 for (j = 0; j <= l; j++)
441 for (j = 0; j < M; j++)
443 plan3D->
p_nfft.
x[3* j ] = plan3D->
x[3* j + 2];
444 plan3D->
p_nfft.
x[3* j + 1] = plan3D->
x[3* j ];
445 plan3D->
p_nfft.
x[3* j + 2] = plan3D->
x[3* j + 1];
453 if ((plan3D->
p_nfft).flags & FG_PSI)
455 nfft_precompute_one_psi(&(plan3D->
p_nfft));
457 if ((plan3D->
p_nfft).flags & PRE_PSI)
459 nfft_precompute_one_psi(&(plan3D->
p_nfft));
466 int i, j, m, k, max, glo1, glo2;
478 for (j = 0; j < M; j++)
479 plan3D->
f[j] = plan3D->
f_hat[0];
486 for (k = -N; k <= N; k++)
488 for (m = -N; m <= N; m++)
491 max = (ABS(m) > ABS(k) ? ABS(m) : ABS(k));
493 for (j = 0; j <= N - max; j++)
497 if ((plan3D->
flags & NFSOFT_NORMALIZED))
500 * SQRT(0.5 * (2. * (max + j) + 1.));
503 if ((plan3D->
flags & NFSOFT_REPRESENT))
505 if ((k < 0) && (k % 2))
509 if ((m < 0) && (m % 2))
520 for (j = N - max + 1; j <
X(next_power_of_2)(N) + 1; j++)
525 c2e(plan3D, ABS((k + m) % 2));
527 for (i = 1; i <= 2* plan3D ->
N_total + 2; i++)
529 plan3D->
p_nfft.
f_hat[NFSOFT_INDEX(k, m, i - N - 1, N) - 1]
530 = plan3D->
cheby[i - 1];
538 if (plan3D->
flags & NFSOFT_USE_NDFT)
540 nfft_trafo_direct(&(plan3D->
p_nfft));
544 nfft_trafo(&(plan3D->
p_nfft));
547 for (j = 0; j < plan3D->
M_total; j++)
565 my_plan->
aux[0]= 1/(2*_Complex_I)*my_plan->
cheby[1];
569 my_plan->
aux[j]=1/(2*_Complex_I)*(my_plan->
cheby[j+1]-my_plan->
cheby[j-1]);
571 my_plan->
aux[N-1]=1/(2*_Complex_I)*(-my_plan->
cheby[j-1]);
582 for(j=1;j<=my_plan->
N_total;j++)
595 int i, j, m, k, max, glo1, glo2;
608 for (j = 0; j < M; j++)
609 plan3D->
f_hat[0] += plan3D->
f[j];
613 for (j = 0; j < M; j++)
618 if (plan3D->
flags & NFSOFT_USE_NDFT)
620 nfft_adjoint_direct(&(plan3D->
p_nfft));
624 nfft_adjoint(&(plan3D->
p_nfft));
631 for (k = -N; k <= N; k++)
633 for (m = -N; m <= N; m++)
636 max = (ABS(m) > ABS(k) ? ABS(m) : ABS(k));
638 for (i = 1; i < 2* plan3D ->
N_total + 3; i++)
646 e2c(plan3D, ABS((k + m) % 2));
655 for (j = max; j <= N; j++)
657 if ((plan3D->
flags & NFSOFT_REPRESENT))
659 if ((k < 0) && (k % 2))
663 if ((m < 0) && (m % 2))
673 if ((plan3D->
flags & NFSOFT_NORMALIZED))
675 plan3D->
f_hat[glo1] = plan3D->
f_hat[glo1] * (1 / (2. * KPI)) * SQRT(
676 0.5 * (2. * (j) + 1.));
689 nfft_finalize(&plan->
p_nfft);
697 if (plan->
flags & NFSOFT_MALLOC_F_HAT)
704 if (plan->
flags & NFSOFT_MALLOC_F)
711 if (plan->
flags & NFSOFT_MALLOC_X)
718 int posN(
int n,
int m,
int B)
723 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