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