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