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);
37 void nfsoft_init(nfsoft_plan *plan,
int N,
int M)
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)
81 plan->f_hat = (C*) nfft_malloc((B + 1) * (4* (B +1)*(B+1)-1)/3*
sizeof(C));
82 if (plan->f_hat == NULL ) printf(
"Allocation failed!\n");
85 if (plan->flags & NFSOFT_MALLOC_X)
87 plan->x = (R*) nfft_malloc(plan->M_total*3*
sizeof(R));
88 if (plan->x == NULL ) printf(
"Allocation failed!\n");
90 if (plan->flags & NFSOFT_MALLOC_F)
92 plan->f = (C*) nfft_malloc(plan->M_total*
sizeof(C));
93 if (plan->f == NULL ) printf(
"Allocation failed!\n");
96 plan->wig_coeffs = (C*) nfft_malloc((
X(next_power_of_2)(B)+1)*
sizeof(C));
97 plan->cheby = (C*) nfft_malloc((2*B+2)*
sizeof(C));
98 plan->aux = (C*) nfft_malloc((2*B+4)*
sizeof(C));
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;
107 plan->internal_fpt_set = SO3_fpt_init(plan->N_total, plan->flags, plan->fpt_kappa);
111 static void c2e(nfsoft_plan *my_plan,
int even)
116 N = 2* (my_plan ->N_total+1);
119 my_plan->cheby[my_plan->N_total+1] = my_plan->wig_coeffs[0];
120 my_plan->cheby[0]=0.0;
122 for (j=1;j<my_plan->N_total+1;j++)
124 my_plan->cheby[my_plan->N_total+1+j]=0.5* my_plan->wig_coeffs[j];
125 my_plan->cheby[my_plan->N_total+1-j]=0.5* my_plan->wig_coeffs[j];
128 C *aux= (C*) nfft_malloc((N+2)*
sizeof(C));
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;
155 R *alpha, *beta, *gamma;
158 if (flags & NFSOFT_USE_DPT)
165 t = (int) log2(
X(next_power_of_2)(N));
174 N =
X(next_power_of_2)(l);
180 alpha = (R*) nfft_malloc((N + 2) *
sizeof(R));
181 beta = (R*) nfft_malloc((N + 2) *
sizeof(R));
182 gamma = (R*) nfft_malloc((N + 2) *
sizeof(R));
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;
222 R *alpha, *beta, *gamma;
226 if (flags & NFSOFT_USE_DPT)
233 t = (int) log2(
X(next_power_of_2)(N));
242 N =
X(next_power_of_2)(l);
248 alpha = (R*) nfft_malloc((N + 2) *
sizeof(R));
249 beta = (R*) nfft_malloc((N + 2) *
sizeof(R));
250 gamma = (R*) nfft_malloc((N + 2) *
sizeof(R));
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);
319 x = (C*) nfft_malloc((k_end + 1) *
sizeof(C));
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);
335 y = (C*) nfft_malloc((k_end + 1) *
sizeof(C));
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);
395 y = (C*) nfft_malloc((k_end + 1) *
sizeof(C));
397 x = (C*) nfft_malloc((k_end + 1) *
sizeof(C));
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++)
431 void nfsoft_precompute(nfsoft_plan *plan3D)
434 int N = plan3D->N_total;
435 int M = plan3D->M_total;
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];
446 for (j = 0; j < 3* plan3D ->p_nfft.M_total; j++)
448 plan3D->p_nfft.x[j] = plan3D->p_nfft.x[j] * (1 / (2* KPI ));
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));
462 void nfsoft_trafo(nfsoft_plan *plan3D)
464 int i, j, m, k, max, glo1, glo2;
470 int N = plan3D->N_total;
471 int M = plan3D->M_total;
476 for (j = 0; j < M; j++)
477 plan3D->f[j] = plan3D->f_hat[0];
481 for (j = 0; j < plan3D->p_nfft.N_total; j++)
482 plan3D->p_nfft.f_hat[j] = 0.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++)
493 plan3D->wig_coeffs[j] = plan3D->f_hat[glo1];
495 if ((plan3D->flags & NFSOFT_NORMALIZED))
497 plan3D->wig_coeffs[j] = plan3D->wig_coeffs[j] * (1. / (2. * KPI))
498 * SQRT(0.5 * (2. * (max + j) + 1.));
501 if ((plan3D->flags & NFSOFT_REPRESENT))
503 if ((k < 0) && (k % 2))
505 plan3D->wig_coeffs[j] = plan3D->wig_coeffs[j] * (-1);
507 if ((m < 0) && (m % 2))
508 plan3D->wig_coeffs[j] = plan3D->wig_coeffs[j] * (-1);
511 plan3D->wig_coeffs[j] = plan3D->wig_coeffs[j] * (-1);
518 for (j = N - max + 1; j <
X(next_power_of_2)(N) + 1; j++)
519 plan3D->wig_coeffs[j] = 0.0;
521 SO3_fpt(plan3D->wig_coeffs, plan3D->internal_fpt_set, N, k, m, plan3D->flags);
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++)
546 plan3D->f[j] = plan3D->p_nfft.f[j];
550 static void e2c(nfsoft_plan *my_plan,
int even)
556 N = 2* (my_plan ->N_total+1);
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]);
574 my_plan->cheby[j]= my_plan->aux[j];
578 my_plan->wig_coeffs[0]=my_plan->cheby[my_plan->N_total+1];
580 for(j=1;j<=my_plan->N_total;j++)
582 my_plan->wig_coeffs[j]=0.5*(my_plan->cheby[my_plan->N_total+j+1]+my_plan->cheby[my_plan->N_total+1-j]);
591 void nfsoft_adjoint(nfsoft_plan *plan3D)
593 int i, j, m, k, max, glo1, glo2;
599 int N = plan3D->N_total;
600 int M = plan3D->M_total;
606 for (j = 0; j < M; j++)
607 plan3D->f_hat[0] += plan3D->f[j];
611 for (j = 0; j < M; j++)
613 plan3D->p_nfft.f[j] = plan3D->f[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++)
638 plan3D->cheby[i - 1] = plan3D->p_nfft.f_hat[NFSOFT_INDEX(k, m, i - N
644 e2c(plan3D, ABS((k + m) % 2));
647 SO3_fpt_transposed(plan3D->wig_coeffs, plan3D->internal_fpt_set, N, k, m,
653 for (j = max; j <= N; j++)
655 if ((plan3D->flags & NFSOFT_REPRESENT))
657 if ((k < 0) && (k % 2))
659 plan3D->wig_coeffs[j] = -plan3D->wig_coeffs[j];
661 if ((m < 0) && (m % 2))
662 plan3D->wig_coeffs[j] = -plan3D->wig_coeffs[j];
665 plan3D->wig_coeffs[j] = plan3D->wig_coeffs[j] * (-1);
669 plan3D->f_hat[glo1] = plan3D->wig_coeffs[j];
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.));
684 void nfsoft_finalize(nfsoft_plan *plan)
687 nfft_finalize(&plan->p_nfft);
688 free(plan->wig_coeffs);
692 fpt_finalize(plan->internal_fpt_set);
693 plan->internal_fpt_set = NULL;
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));
#define X(name)
Include header for C99 complex datatype.
static void c2e(nfsft_plan *plan)
Converts coefficients with , from a linear combination of Chebyshev polynomials.
Header file for functions related to Wigner-d/D functions.
void SO3_beta_row(double *beta, int N, int m, int n)
Compute three-term-recurrence coefficients of Wigner-d functions for all degrees .
void SO3_gamma_row(double *gamma, int N, int m, int n)
Compute three-term-recurrence coefficients of Wigner-d functions for all degrees .
void SO3_alpha_row(double *alpha, int N, int m, int n)
Compute three-term-recurrence coefficients of Wigner-d functions for all degrees .