NFFT Logo 3.2.2
nfsoft.c
00001 /*
00002  * Copyright (c) 2002, 2012 Jens Keiner, Stefan Kunis, Daniel Potts
00003  *
00004  * This program is free software; you can redistribute it and/or modify it under
00005  * the terms of the GNU General Public License as published by the Free Software
00006  * Foundation; either version 2 of the License, or (at your option) any later
00007  * version.
00008  *
00009  * This program is distributed in the hope that it will be useful, but WITHOUT
00010  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00011  * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
00012  * details.
00013  *
00014  * You should have received a copy of the GNU General Public License along with
00015  * this program; if not, write to the Free Software Foundation, Inc., 51
00016  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
00017  */
00018 
00019 /* $Id: nfsoft.c 3896 2012-10-10 12:19:26Z tovo $ */
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     int rr;
00269     for (rr = 0; rr < N + 2; rr++)
00270       fprintf(stderr, "a[%4d] = %10e b[%4d] = %10e c[%4d] = %10e\n",rr,alpha[rr],rr,beta[rr],rr,gamma[rr]);
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       //fprintf(stdout,"\n k= %d, m= %d \n",k,m);
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         //fprintf(stdout,"%f \t", plan3D->nfft_plan.f_hat[NFSOFT_INDEX(k,m,i-N-1,N)-1]);
00537         //fprintf(stdout,"another index: %d for k=%d,m=%d,l=%d,N=%d \n", NFSOFT_INDEX(k,m,i-N-1,N)-1,k,m,i-N-1,N);
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   //nfft_vpr_complex(my_plan->cheby,N+1,"chebychev");
00565 
00566 
00567       if (even>0)
00568       {
00569         //my_plan->aux[N-1]= -1/(2*I)* my_plan->cheby[N-2];
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 //nfft_vpr_complex(my_plan->wig_coeffs,my_plan->N_total,"chebychev ");
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   //nothing much to be done for polynomial degree 0
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   //nfft_vpr_complex(plan3D->nfft_plan.f_hat,plan3D->nfft_plan.N_total,"all results");
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       //fprintf(stdout,"k=%d,m=%d \n",k,m);
00650       //nfft_vpr_complex(plan3D->cheby,2*plan3D->N_total+2,"euler");
00651       e2c(plan3D, ABS((k + m) % 2));
00652 
00653       //nfft_vpr_complex(plan3D->wig_coeffs,plan3D->N_total+1,"chebys");
00654       SO3_fpt_transposed(plan3D->wig_coeffs, plan3D->internal_fpt_set, N, k, m,
00655           plan3D->flags);
00656       //nfft_vpr_complex(plan3D->wig_coeffs,plan3D->N_total+1,"wigners");
00657       //  SO3_fpt_transposed(plan3D->wig_coeffs,N,k,m,plan3D->flags,plan3D->fpt_kappa);
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   /* Finalise the nfft plan. */
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     //fprintf(stderr,"deallocating f_hat\n");
00705     free(plan->f_hat);
00706   }
00707 
00708   /* De-allocate memory for samples, if neccesary. */
00709   if (plan->flags & NFSOFT_MALLOC_F)
00710   {
00711     //fprintf(stderr,"deallocating f\n");
00712     free(plan->f);
00713   }
00714 
00715   /* De-allocate memory for nodes, if neccesary. */
00716   if (plan->flags & NFSOFT_MALLOC_X)
00717   {
00718     //fprintf(stderr,"deallocating x\n");
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   //(n > -B? pos=posN(n-1,m,B)+B+1-MAX(ABS(m),ABS(n-1)): pos= 0)
00732   return pos;
00733 }
00734 

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409