00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "config.h"
00022
00023 #include <stdio.h>
00024 #include <math.h>
00025 #include <string.h>
00026 #include <stdlib.h>
00027 #ifdef HAVE_COMPLEX_H
00028 #include <complex.h>
00029 #endif
00030 #include "nfft3.h"
00031 #include "nfft3util.h"
00032 #include "infft.h"
00033 #include "wigner.h"
00034
00035 #define DEFAULT_NFFT_CUTOFF 6
00036 #define FPT_THRESHOLD 1000
00037
00038 void nfsoft_init(nfsoft_plan *plan, int N, int M)
00039 {
00040 nfsoft_init_advanced(plan, N, M, NFSOFT_MALLOC_X | NFSOFT_MALLOC_F
00041 | NFSOFT_MALLOC_F_HAT);
00042 }
00043
00044 void nfsoft_init_advanced(nfsoft_plan *plan, int N, int M,
00045 unsigned int nfsoft_flags)
00046 {
00047 nfsoft_init_guru(plan, N, M, nfsoft_flags, PRE_PHI_HUT | PRE_PSI | MALLOC_X
00048 | MALLOC_F_HAT | MALLOC_F | FFTW_INIT | FFT_OUT_OF_PLACE,
00049 DEFAULT_NFFT_CUTOFF, FPT_THRESHOLD);
00050 }
00051
00052 void nfsoft_init_guru(nfsoft_plan *plan, int B, int M,
00053 unsigned int nfsoft_flags, unsigned int nfft_flags, int nfft_cutoff,
00054 int fpt_kappa)
00055 {
00056 int N[3];
00057 int n[3];
00058
00059 N[0] = 2* B + 2;
00060 N[1] = 2* B + 2;
00061 N[2] = 2* B + 2;
00062
00063 n[0] = 8* B ;
00064 n[1] = 8* B ;
00065 n[2] = 8* B ;
00066
00067 nfft_init_guru(&plan->p_nfft, 3, N, M, n, nfft_cutoff, nfft_flags,
00068 FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
00069
00070 if ((plan->p_nfft).nfft_flags & PRE_LIN_PSI)
00071 {
00072 nfft_precompute_lin_psi(&(plan->p_nfft));
00073 }
00074
00075 plan->N_total = B;
00076 plan->M_total = M;
00077 plan->fpt_kappa = fpt_kappa;
00078 plan->flags = nfsoft_flags;
00079
00080 if (plan->flags & NFSOFT_MALLOC_F_HAT)
00081 {
00082 plan->f_hat = (C*) nfft_malloc((B + 1) * (4* (B +1)*(B+1)-1)/3*sizeof(C));
00083 if (plan->f_hat == NULL ) printf("Allocation failed!\n");
00084 }
00085
00086 if (plan->flags & NFSOFT_MALLOC_X)
00087 {
00088 plan->x = (R*) nfft_malloc(plan->M_total*3*sizeof(R));
00089 if (plan->x == NULL ) printf("Allocation failed!\n");
00090 }
00091 if (plan->flags & NFSOFT_MALLOC_F)
00092 {
00093 plan->f = (C*) nfft_malloc(plan->M_total*sizeof(C));
00094 if (plan->f == NULL ) printf("Allocation failed!\n");
00095 }
00096
00097 plan->wig_coeffs = (C*) nfft_malloc((X(next_power_of_2)(B)+1)*sizeof(C));
00098 plan->cheby = (C*) nfft_malloc((2*B+2)*sizeof(C));
00099 plan->aux = (C*) nfft_malloc((2*B+4)*sizeof(C));
00100
00101 if (plan->wig_coeffs == NULL ) printf("Allocation failed!\n");
00102 if (plan->cheby == NULL ) printf("Allocation failed!\n");
00103 if (plan->aux == NULL ) printf("Allocation failed!\n");
00104
00105 plan->mv_trafo = (void (*) (void* ))nfsoft_trafo;
00106 plan->mv_adjoint = (void (*) (void* ))nfsoft_adjoint;
00107
00108 plan->internal_fpt_set = 0;
00109 }
00110
00111 static void c2e(nfsoft_plan *my_plan, int even)
00112 {
00113 int j, N;
00114
00116 N = 2* (my_plan ->N_total+1);
00117
00119 my_plan->cheby[my_plan->N_total+1] = my_plan->wig_coeffs[0];
00120 my_plan->cheby[0]=0.0;
00121
00122 for (j=1;j<my_plan->N_total+1;j++)
00123 {
00124 my_plan->cheby[my_plan->N_total+1+j]=0.5* my_plan->wig_coeffs[j];
00125 my_plan->cheby[my_plan->N_total+1-j]=0.5* my_plan->wig_coeffs[j];
00126 }
00127
00128 C *aux= (C*) nfft_malloc((N+2)*sizeof(C));
00129
00130 for(j=1;j<N;j++)
00131 aux[j]=my_plan->cheby[j];
00132
00133 aux[0]=0.;
00134 aux[N]=0.;
00135
00136 if (even>0)
00137 {
00138 my_plan->cheby[0]=(C) (-1.)/(2.*_Complex_I) * aux[1];
00139 for (j=1;j<N;j++)
00140 {
00141 my_plan->cheby[j]=(1./(2.*_Complex_I)*(aux[j+1]-aux[j-1]));
00142 }
00143
00144 }
00145 free(aux);
00146 aux = NULL;
00147 }
00148
00149
00150 static fpt_set SO3_fpt_init(int l, fpt_set set, unsigned int flags, int kappa)
00151 {
00152 int N, t, k_start, k_end, k, m;
00153 int glo = 0;
00154 R *alpha, *beta, *gamma;
00155
00157 if (flags & NFSOFT_USE_DPT)
00158 {
00159 if (l < 2)
00160 N = 2;
00161 else
00162 N = l;
00163
00164 t = (int) log2(X(next_power_of_2)(N));
00165
00166 }
00167 else
00168 {
00170 if (l < 2)
00171 N = 2;
00172 else
00173 N = X(next_power_of_2)(l);
00174
00175 t = (int) log2(N);
00176 }
00177
00179 alpha = (R*) nfft_malloc((N + 2) * sizeof(R));
00180 beta = (R*) nfft_malloc((N + 2) * sizeof(R));
00181 gamma = (R*) nfft_malloc((N + 2) * sizeof(R));
00182
00184 if (flags & NFSOFT_NO_STABILIZATION)
00185 {
00186 set = fpt_init((2* N + 1) * (2* N + 1), t, 0U | FPT_NO_STABILIZATION);
00187 }
00188 else
00189 {
00190 set = fpt_init((2* N + 1) * (2* N + 1), t, 0U);
00191 }
00192
00193 for (k = -N; k <= N; k++)
00194 for (m = -N; m <= N; m++)
00195 {
00197 k_start = (ABS(k) >= ABS(m)) ? ABS(k) : ABS(m);
00198 k_end = N;
00199
00200 SO3_alpha_row(alpha, N, k, m);
00201 SO3_beta_row(beta, N, k, m);
00202 SO3_gamma_row(gamma, N, k, m);
00203
00204 fpt_precompute(set, glo, alpha, beta, gamma, k_start, kappa);
00205 glo++;
00206 }
00207
00208 free(alpha);
00209 free(beta);
00210 free(gamma);
00211 alpha = NULL;
00212 beta = NULL;
00213 gamma = NULL;
00214
00215 return set;
00216 }
00217
00218 static fpt_set SO3_single_fpt_init(int l, int k, int m, unsigned int flags, int kappa)
00219 {
00220 int N, t, k_start, k_end;
00221 R *alpha, *beta, *gamma;
00222 fpt_set set = 0;
00223
00225 if (flags & NFSOFT_USE_DPT)
00226 {
00227 if (l < 2)
00228 N = 2;
00229 else
00230 N = l;
00231
00232 t = (int) log2(X(next_power_of_2)(N));
00233
00234 }
00235 else
00236 {
00238 if (l < 2)
00239 N = 2;
00240 else
00241 N = X(next_power_of_2)(l);
00242
00243 t = (int) log2(N);
00244 }
00245
00247 alpha = (R*) nfft_malloc((N + 2) * sizeof(R));
00248 beta = (R*) nfft_malloc((N + 2) * sizeof(R));
00249 gamma = (R*) nfft_malloc((N + 2) * sizeof(R));
00250
00252 {
00253 unsigned int fptflags = 0U
00254 | IF(flags & NFSOFT_USE_DPT,FPT_NO_FAST_ALGORITHM,IF(t > 1,FPT_NO_DIRECT_ALGORITHM,0U))
00255 | IF(flags & NFSOFT_NO_STABILIZATION,FPT_NO_STABILIZATION,0U);
00256 set = fpt_init(1, t, fptflags);
00257 }
00258
00260 k_start = (ABS(k) >= ABS(m)) ? ABS(k) : ABS(m);
00261 k_end = N;
00262
00263 SO3_alpha_row(alpha, N, k, m);
00264 SO3_beta_row(beta, N, k, m);
00265 SO3_gamma_row(gamma, N, k, m);
00266
00267
00268
00269
00270
00271
00272
00273 fpt_precompute(set, 0, alpha, beta, gamma, k_start, kappa);
00274
00275 free(alpha);
00276 free(beta);
00277 free(gamma);
00278 alpha = NULL;
00279 beta = NULL;
00280 gamma = NULL;
00281
00282 return set;
00283 }
00284
00285 void SO3_fpt(C *coeffs, fpt_set set, int l, int k, int m, unsigned int flags)
00286 {
00287 int N;
00289 C* x;
00291 C* y;
00292
00293 int trafo_nr;
00294 int k_start, k_end, j;
00295 int function_values = 0;
00296
00298 if (flags & NFSOFT_USE_DPT)
00299 {
00300 N = l;
00301 if (l < 2)
00302 N = 2;
00303 }
00304 else
00305 {
00306 if (l < 2)
00307 N = 2;
00308 else
00309 N = X(next_power_of_2)(l);
00310 }
00311
00313 k_start = (ABS(k) >= ABS(m)) ? ABS(k) : ABS(m);
00314 k_end = N;
00315 trafo_nr = (N + k) * (2* N + 1) + (m + N);
00316
00318 x = (C*) nfft_malloc((k_end + 1) * sizeof(C));
00319
00320 for (j = 0; j <= k_end; j++)
00321 x[j] = K(0.0);
00322
00323
00324 for (j = 0; j <= l - k_start; j++)
00325 {
00326 x[j + k_start] = coeffs[j];
00327 }
00328 for (j = l - k_start + 1; j <= k_end - k_start; j++)
00329 {
00330 x[j + k_start] = K(0.0);
00331 }
00332
00334 y = (C*) nfft_malloc((k_end + 1) * sizeof(C));
00335
00336 if (flags & NFSOFT_USE_DPT)
00337 {
00338 fpt_trafo_direct(set, trafo_nr, &x[k_start], y, k_end, 0U
00339 | (function_values ? FPT_FUNCTION_VALUES : 0U));
00340 }
00341 else
00342 {
00343 fpt_trafo(set, trafo_nr, &x[k_start], y, k_end, 0U
00344 | (function_values ? FPT_FUNCTION_VALUES : 0U));
00345 }
00346
00348 for (j = 0; j <= l; j++)
00349 {
00350 coeffs[j] = y[j];
00351 }
00352
00355 free(x);
00356 free(y);
00357 x = NULL;
00358 y = NULL;
00359 }
00360
00361 void SO3_fpt_transposed(C *coeffs, fpt_set set, int l, int k, int m,
00362 unsigned int flags)
00363 {
00364 int N, k_start, k_end, j;
00365 int trafo_nr;
00366 int function_values = 0;
00368 C* x;
00370 C* y;
00371
00374 if (flags & NFSOFT_USE_DPT)
00375 {
00376 N = l;
00377 if (l < 2)
00378 N = 2;
00379 }
00380 else
00381 {
00382 if (l < 2)
00383 N = 2;
00384 else
00385 N = X(next_power_of_2)(l);
00386 }
00387
00389 k_start = (ABS(k) >= ABS(m)) ? ABS(k) : ABS(m);
00390 k_end = N;
00391 trafo_nr = (N + k) * (2* N + 1) + (m + N);
00392
00394 y = (C*) nfft_malloc((k_end + 1) * sizeof(C));
00396 x = (C*) nfft_malloc((k_end + 1) * sizeof(C));
00397
00398 for (j = 0; j <= l; j++)
00399 {
00400 y[j] = coeffs[j];
00401 }
00402 for (j = l + 1; j <= k_end; j++)
00403 {
00404 y[j] = K(0.0);
00405 }
00406
00407 if (flags & NFSOFT_USE_DPT)
00408 {
00409 fpt_transposed_direct(set, trafo_nr, &x[k_start], y, k_end, 0U
00410 | (function_values ? FPT_FUNCTION_VALUES : 0U));
00411 }
00412 else
00413 {
00414 fpt_transposed(set, trafo_nr, &x[k_start], y, k_end, 0U
00415 | (function_values ? FPT_FUNCTION_VALUES : 0U));
00416 }
00417
00418 for (j = 0; j <= l; j++)
00419 {
00420 coeffs[j] = x[j];
00421 }
00422
00424 free(x);
00425 free(y);
00426 x = NULL;
00427 y = NULL;
00428 }
00429
00430 void nfsoft_precompute(nfsoft_plan *plan3D)
00431 {
00432 int j;
00433 int N = plan3D->N_total;
00434 int M = plan3D->M_total;
00435
00438 for (j = 0; j < M; j++)
00439 {
00440 plan3D->p_nfft.x[3* j ] = plan3D->x[3* j + 2];
00441 plan3D->p_nfft.x[3* j + 1] = plan3D->x[3* j ];
00442 plan3D->p_nfft.x[3* j + 2] = plan3D->x[3* j + 1];
00443 }
00444
00445 for (j = 0; j < 3* plan3D ->p_nfft.M_total; j++)
00446 {
00447 plan3D->p_nfft.x[j] = plan3D->p_nfft.x[j] * (1 / (2* PI ));
00448 }
00449
00450 if ((plan3D->p_nfft).nfft_flags & FG_PSI)
00451 {
00452 nfft_precompute_one_psi(&(plan3D->p_nfft));
00453 }
00454 if ((plan3D->p_nfft).nfft_flags & PRE_PSI)
00455 {
00456 nfft_precompute_one_psi(&(plan3D->p_nfft));
00457 }
00458
00460 plan3D->internal_fpt_set = SO3_fpt_init(N, plan3D->internal_fpt_set,
00461 plan3D->flags, plan3D->fpt_kappa);
00462
00463 if ((plan3D->p_nfft).nfft_flags & MALLOC_F_HAT)
00464 for (j = 0; j < plan3D->p_nfft.N_total; j++)
00465 plan3D->p_nfft.f_hat[j] = 0.0;
00466
00467 }
00468
00469 void nfsoft_trafo(nfsoft_plan *plan3D)
00470 {
00471 int i, j, m, k, max, glo1, glo2;
00472
00473 i = 0;
00474 glo1 = 0;
00475 glo2 = 0;
00476
00477 int N = plan3D->N_total;
00478 int M = plan3D->M_total;
00479
00481 if (N == 0)
00482 {
00483 for (j = 0; j < M; j++)
00484 plan3D->f[j] = plan3D->f_hat[0];
00485 return;
00486 }
00487
00488 for (j = 0; j < plan3D->p_nfft.N_total; j++)
00489 plan3D->p_nfft.f_hat[j] = 0.0;
00490
00491 for (k = -N; k <= N; k++)
00492 {
00493 for (m = -N; m <= N; m++)
00494 {
00495
00496 max = (ABS(m) > ABS(k) ? ABS(m) : ABS(k));
00497
00498 for (j = 0; j <= N - max; j++)
00499 {
00500 plan3D->wig_coeffs[j] = plan3D->f_hat[glo1];
00501
00502 if ((plan3D->flags & NFSOFT_NORMALIZED))
00503 {
00504 plan3D->wig_coeffs[j] = plan3D->wig_coeffs[j] * (1. / (2. * PI))
00505 * SQRT(0.5 * (2. * (max + j) + 1.));
00506 }
00507
00508 if ((plan3D->flags & NFSOFT_REPRESENT))
00509 {
00510 if ((k < 0) && (k % 2))
00511 {
00512 plan3D->wig_coeffs[j] = plan3D->wig_coeffs[j] * (-1);
00513 }
00514 if ((m < 0) && (m % 2))
00515 plan3D->wig_coeffs[j] = plan3D->wig_coeffs[j] * (-1);
00516
00517 if ((m + k) % 2)
00518 plan3D->wig_coeffs[j] = plan3D->wig_coeffs[j] * (-1);
00519
00520 }
00521
00522 glo1++;
00523 }
00524
00525 for (j = N - max + 1; j < X(next_power_of_2)(N) + 1; j++)
00526 plan3D->wig_coeffs[j] = 0.0;
00527
00528 SO3_fpt(plan3D->wig_coeffs, plan3D->internal_fpt_set, N, k, m, plan3D->flags);
00529
00530 c2e(plan3D, ABS((k + m) % 2));
00531
00532 for (i = 1; i <= 2* plan3D ->N_total + 2; i++)
00533 {
00534 plan3D->p_nfft.f_hat[NFSOFT_INDEX(k, m, i - N - 1, N) - 1]
00535 = plan3D->cheby[i - 1];
00536
00537
00538 }
00539
00540 }
00541 }
00542
00543 if (plan3D->flags & NFSOFT_USE_NDFT)
00544 {
00545 nfft_trafo_direct(&(plan3D->p_nfft));
00546 }
00547 else
00548 {
00549 nfft_trafo(&(plan3D->p_nfft));
00550 }
00551
00552 for (j = 0; j < plan3D->M_total; j++)
00553 plan3D->f[j] = plan3D->p_nfft.f[j];
00554
00555 }
00556
00557 static void e2c(nfsoft_plan *my_plan, int even)
00558 {
00559 int N;
00560 int j;
00561
00563 N = 2* (my_plan ->N_total+1);
00564
00565
00566
00567 if (even>0)
00568 {
00569
00570 my_plan->aux[0]= 1/(2*_Complex_I)*my_plan->cheby[1];
00571
00572 for(j=1;j<N-1;j++)
00573 {
00574 my_plan->aux[j]=1/(2*_Complex_I)*(my_plan->cheby[j+1]-my_plan->cheby[j-1]);
00575 }
00576 my_plan->aux[N-1]=1/(2*_Complex_I)*(-my_plan->cheby[j-1]);
00577
00578
00579 for(j=0;j<N;j++)
00580 {
00581 my_plan->cheby[j]= my_plan->aux[j];
00582 }
00583 }
00584
00585 my_plan->wig_coeffs[0]=my_plan->cheby[my_plan->N_total+1];
00586
00587 for(j=1;j<=my_plan->N_total;j++)
00588 {
00589 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]);
00590 }
00591
00592
00593
00594
00595
00596 }
00597
00598 void nfsoft_adjoint(nfsoft_plan *plan3D)
00599 {
00600 int i, j, m, k, max, glo1, glo2;
00601
00602 i = 0;
00603 glo1 = 0;
00604 glo2 = 0;
00605
00606 int N = plan3D->N_total;
00607 int M = plan3D->M_total;
00608
00609
00610 if (N == 0)
00611 {
00612 plan3D->f_hat[0]=0;
00613 for (j = 0; j < M; j++)
00614 plan3D->f_hat[0] += plan3D->f[j];
00615 return;
00616 }
00617
00618 for (j = 0; j < M; j++)
00619 {
00620 plan3D->p_nfft.f[j] = plan3D->f[j];
00621 }
00622
00623 if (plan3D->flags & NFSOFT_USE_NDFT)
00624 {
00625 nfft_adjoint_direct(&(plan3D->p_nfft));
00626 }
00627 else
00628 {
00629 nfft_adjoint(&(plan3D->p_nfft));
00630 }
00631
00632
00633
00634 glo1 = 0;
00635
00636 for (k = -N; k <= N; k++)
00637 {
00638 for (m = -N; m <= N; m++)
00639 {
00640
00641 max = (ABS(m) > ABS(k) ? ABS(m) : ABS(k));
00642
00643 for (i = 1; i < 2* plan3D ->N_total + 3; i++)
00644 {
00645 plan3D->cheby[i - 1] = plan3D->p_nfft.f_hat[NFSOFT_INDEX(k, m, i - N
00646 - 1, N) - 1];
00647 }
00648
00649
00650
00651 e2c(plan3D, ABS((k + m) % 2));
00652
00653
00654 SO3_fpt_transposed(plan3D->wig_coeffs, plan3D->internal_fpt_set, N, k, m,
00655 plan3D->flags);
00656
00657
00658
00659
00660 for (j = max; j <= N; j++)
00661 {
00662 if ((plan3D->flags & NFSOFT_REPRESENT))
00663 {
00664 if ((k < 0) && (k % 2))
00665 {
00666 plan3D->wig_coeffs[j] = -plan3D->wig_coeffs[j];
00667 }
00668 if ((m < 0) && (m % 2))
00669 plan3D->wig_coeffs[j] = -plan3D->wig_coeffs[j];
00670
00671 if ((m + k) % 2)
00672 plan3D->wig_coeffs[j] = plan3D->wig_coeffs[j] * (-1);
00673
00674 }
00675
00676 plan3D->f_hat[glo1] = plan3D->wig_coeffs[j];
00677
00678 if ((plan3D->flags & NFSOFT_NORMALIZED))
00679 {
00680 plan3D->f_hat[glo1] = plan3D->f_hat[glo1] * (1 / (2. * PI)) * SQRT(
00681 0.5 * (2. * (j) + 1.));
00682 }
00683
00684 glo1++;
00685 }
00686
00687 }
00688 }
00689 }
00690
00691 void nfsoft_finalize(nfsoft_plan *plan)
00692 {
00693
00694 nfft_finalize(&plan->p_nfft);
00695 free(plan->wig_coeffs);
00696 free(plan->cheby);
00697 free(plan->aux);
00698
00699 fpt_finalize(plan->internal_fpt_set);
00700 plan->internal_fpt_set = NULL;
00701
00702 if (plan->flags & NFSOFT_MALLOC_F_HAT)
00703 {
00704
00705 free(plan->f_hat);
00706 }
00707
00708
00709 if (plan->flags & NFSOFT_MALLOC_F)
00710 {
00711
00712 free(plan->f);
00713 }
00714
00715
00716 if (plan->flags & NFSOFT_MALLOC_X)
00717 {
00718
00719 free(plan->x);
00720 }
00721 }
00722
00723 int posN(int n, int m, int B)
00724 {
00725 int pos;
00726
00727 if (n > -B)
00728 pos = posN(n - 1, m, B) + B + 1 - MAX(ABS(m), ABS(n - 1));
00729 else
00730 pos = 0;
00731
00732 return pos;
00733 }
00734