NFFT Logo 3.2.2
nfst.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: nfst.c 3896 2012-10-10 12:19:26Z tovo $ */
00020 
00027 #include <stdio.h>
00028 #include <math.h>
00029 #include <string.h>
00030 #include <stdlib.h>
00031 
00032 #include "nfft3util.h"
00033 #include "nfft3.h"
00034 #include "infft.h"
00035 
00039 #define NFST_DEFAULT_FLAGS   PRE_PHI_HUT|\
00040                              PRE_PSI|\
00041                              MALLOC_X|\
00042                              MALLOC_F_HAT|\
00043                              MALLOC_F|\
00044                              FFTW_INIT|\
00045                              FFT_OUT_OF_PLACE
00046 
00047 #define FFTW_DEFAULT_FLAGS   FFTW_ESTIMATE|\
00048                              FFTW_DESTROY_INPUT
00049 
00050 #define NFST_SUMMANDS ( 2 * ths->m + 2)
00051 #define NODE(p,r) ( ths->x[(p) * ths->d + (r)])
00052 
00053 #define MACRO_ndst_init_result_trafo      \
00054   memset( f, 0, ths->M_total * sizeof( double));
00055 #define MACRO_ndst_init_result_adjoint \
00056   memset( f_hat, 0, ths->N_total * sizeof( double));
00057 
00058 
00059 #define MACRO_nfst_D_init_result_A        \
00060   memset(g_hat, 0, nfst_prod_minus_a_int( ths->n, 0, ths->d) * sizeof( double));
00061 #define MACRO_nfst_D_init_result_T        \
00062   memset(f_hat, 0, ths->N_total * sizeof( double));
00063 
00064 #define MACRO_nfst_B_init_result_A        \
00065   memset(f, 0,     ths->M_total * sizeof( double));
00066 #define MACRO_nfst_B_init_result_T        \
00067   memset(g, 0, nfst_prod_minus_a_int( ths->n, 0, ths->d) * sizeof( double));
00068 
00069 
00070 #define NFST_PRE_WINFUN( d)  ths->N[d] = 2 * ths->N[d];         \
00071                              ths->n[d] = nfst_fftw_2N( ths->n[d]);
00072 
00073 #define NFST_POST_WINFUN( d) ths->N[d] = (LRINT(0.5 * ths->N[d]));   \
00074                              ths->n[d] = nfst_fftw_2N_rev( ths->n[d]);
00075 
00076 
00077 #define NFST_WINDOW_HELP_INIT  WINDOW_HELP_INIT
00078 
00079 
00080 double nfst_phi_hut( nfst_plan *ths, int k, int d)
00081 {
00082   NFST_PRE_WINFUN( d);
00083   double phi_hut_tmp = PHI_HUT( k, d);
00084   NFST_POST_WINFUN( d);
00085 
00086   return phi_hut_tmp;
00087 }
00088 
00089 double nfst_phi( nfst_plan *ths, double x, int d)
00090 {
00091   NFST_PRE_WINFUN( d);
00092   double phi_tmp = PHI( x, d);
00093   NFST_POST_WINFUN( d);
00094 
00095   return phi_tmp;
00096 }
00097 
00098 int nfst_fftw_2N( int n)
00099 {
00100   return 2 * ( n + 1);
00101 }
00102 
00103 int nfst_fftw_2N_rev( int n)
00104 {
00105   div_t n_div;
00106 
00107   n_div = div(n, 2);
00108   return n_div.quot - 1;
00109 }
00110 
00111 #define MACRO_with_sin_vec     sin_vec[t][ka[t]]
00112 #define MACRO_without_sin_vec  sin( 2.0 * PI * (ka[t]+1) * NODE(j,t))
00113 
00114 
00115 #define MACRO_with_PRE_PHI_HUT     ths->c_phi_inv[t][kg[t]];
00116 #define MACRO_compute_PHI_HUT_INV  (1.0 / (nfst_phi_hut( ths, kg[t]+1, t)))
00117 
00118 #define MACRO_with_PRE_PSI     ths->psi[(j * ths->d + t) * NFST_SUMMANDS + lc[t]];
00119 #define MACRO_compute_PSI      \
00120   nfst_phi( ths, NODE(j,t) - (( double)(lc[t] + lb[t])) / nfst_fftw_2N( ths->n[t]), t)
00121 
00122 
00123 
00139 #define MACRO_ndst_malloc__sin_vec                                              \
00140                                                                                 \
00141   double **sin_vec;                                                             \
00142   sin_vec = (double**)nfft_malloc( ths->d * sizeof( double*));                       \
00143   for( t = 0; t < ths->d; t++)                                                  \
00144     sin_vec[t] = (double*)nfft_malloc( ( ths->N[t] - 1) * sizeof( double));          \
00145 
00146 
00147 
00148 
00149 #define MACRO_ndst_free__sin_vec                                                \
00150 {                                                                               \
00151   /* free allocated memory */                                                   \
00152   for( t = 0; t < ths->d; t++)                                                  \
00153   nfft_free( sin_vec[t]);                                                          \
00154   nfft_free( sin_vec);                                                               \
00155 }
00156 
00157 
00158 
00159 #define MACRO_ndst_init__sin_vec                                                \
00160 {                                                                               \
00161   for( t = 0; t < ths->d; t++)                                                  \
00162   {                                                                             \
00163     cos_x[t] = cos( 2.0 * PI * NODE(j,t));                                      \
00164     sin_vec[t][0] = sin( 2.0 * PI * NODE(j,t));                                 \
00165     sin_vec[t][1] = sin( 4.0 * PI * NODE(j,t));                                 \
00166     for( k = 2; k < ths->N[t] - 1; k++)                                         \
00167       sin_vec[t][k] = 2.0 * cos_x[t] * sin_vec[t][k-1]                          \
00168                       - sin_vec[t][k-2];                                        \
00169   }                                                                             \
00170 }
00171 
00172 
00173 #define MACRO_ndst_init__k__sin_k( which_one)                                   \
00174 {                                                                               \
00175   sin_k[0] = 1.0;                                                               \
00176   for( t = 0; t < ths->d; t++)                                                  \
00177     ka[t] = 0;                                                                  \
00178                                                                                 \
00179   for( t = 0; t < ths->d; t++)                                                  \
00180   {                                                                             \
00181     sin_k[t+1] = sin_k[t] * MACRO_ ##which_one;                                 \
00182   }                                                                             \
00183 }
00184 
00185 
00186 #define MACRO_ndst_count__k__sin_k( which_one)                                  \
00187 {                                                                               \
00188   ka[ths->d-1]++;                                                               \
00189   i = ths->d - 1;                                                               \
00190   while( ( ka[i] == ths->N[i] - 1) && ( i > 0))                                 \
00191   {                                                                             \
00192     ka[i - 1]++;                                                                \
00193     ka[i] = 0;                                                                  \
00194                                                                                 \
00195     i--;                                                                        \
00196   }                                                                             \
00197   for( t = i; t < ths->d; t++)                                                  \
00198     sin_k[t+1] = sin_k[t] * MACRO_ ##which_one;                                 \
00199 }
00200 
00201 
00202 #define MACRO_ndst_compute__trafo                                               \
00203 {                                                                               \
00204   f[j] += f_hat[k] * sin_k[ths->d];                                             \
00205 }
00206 
00207 #define MACRO_ndst_compute__adjoint                                             \
00208 {                                                                               \
00209   f_hat[k] += f[j] * sin_k[ths->d];                                             \
00210 }
00211 
00212 
00213 /* slow (trafo) transform */
00214 #define MACRO_ndst( which_one)                                                  \
00215   void nfst_ ## which_one ## _direct ( nfst_plan *ths)                                     \
00216   {                                                                             \
00217     int j, k, t, i;                                                             \
00218     int ka[ths->d];                                                             \
00219     double sin_k[ths->d+1];                                                     \
00220     double cos_x[ths->d];                                                       \
00221                                                                                 \
00222     double *f     = ths->f;                                                     \
00223     double *f_hat = ths->f_hat;                                                 \
00224                                                                                 \
00225     MACRO_ndst_init_result_ ## which_one;                                       \
00226                                                                                 \
00227     if( ths->d == 1)                                                            \
00228       for( j = 0; j < ths->M_total; j++)                                        \
00229         for( k = 0; k < ths->N_total; k++)                                      \
00230         {                                                                       \
00231           sin_k[ths->d] = sin( 2.0 * PI * (k+1) * NODE(j,0));                   \
00232           MACRO_ndst_compute__ ## which_one;                                    \
00233         }                                                                       \
00234     else                                                                        \
00235       if( 1 == 0) /*FIXME: remove slow slow ... */                              \
00236         /* slow ndst */                                                         \
00237         for( j = 0; j < ths->M_total; j++)                                      \
00238         {                                                                       \
00239           MACRO_ndst_init__k__sin_k(without_sin_vec);                           \
00240                                                                                 \
00241           for( k = 0; k < ths->N_total; k++)                                    \
00242           {                                                                     \
00243             MACRO_ndst_compute__ ## which_one;                                  \
00244                                                                                 \
00245             MACRO_ndst_count__k__sin_k(without_sin_vec);                        \
00246           }                                                                     \
00247         }                                                                       \
00248       else                                                                      \
00249       {                                                                         \
00250         /* fast nfst_trafo_direct */                                                   \
00251         MACRO_ndst_malloc__sin_vec;                                             \
00252                                                                                 \
00253         for( j = 0; j < ths->M_total; j++)                                      \
00254         {                                                                       \
00255           MACRO_ndst_init__sin_vec;                                             \
00256                                                                                 \
00257           MACRO_ndst_init__k__sin_k(with_sin_vec);                              \
00258                                                                                 \
00259           for( k = 0; k < ths->N_total; k++)                                    \
00260           {                                                                     \
00261             MACRO_ndst_compute__ ## which_one;                                  \
00262                                                                                 \
00263             MACRO_ndst_count__k__sin_k(with_sin_vec);                           \
00264           }                                                                     \
00265         }                                                                       \
00266         MACRO_ndst_free__sin_vec;                                               \
00267       }                                                                         \
00268   } /* ndst_{trafo, adjoint} */
00269 
00270 
00271 MACRO_ndst(trafo)
00272 MACRO_ndst(adjoint)
00273 
00274 
00275 
00276 
00291 #define MACRO_nfst__lower_boundary( j,act_dim)                                  \
00292 {                                                                               \
00293   lb[(act_dim)] =                                                               \
00294     (LRINT(NODE((j),(act_dim)) * nfst_fftw_2N( ths->n[(act_dim)]))) - ths->m;   \
00295 }
00296 
00297 #define MACRO_nfst_D_compute_A                                                  \
00298 {                                                                               \
00299   g_hat[kg_plain[ths->d]] = f_hat[k_L] * c_phi_inv_k[ths->d];                   \
00300 }
00301 
00302 #define MACRO_nfst_D_compute_T                                                  \
00303 {                                                                               \
00304   f_hat[k_L] = g_hat[kg_plain[ths->d]] * c_phi_inv_k[ths->d];                   \
00305 }
00306 
00307 
00308 #define MACRO_init__kg                                                          \
00309 {                                                                               \
00310   for( t = 0; t < ths->d; t++)                                                  \
00311     kg[t] = 0;                                                                  \
00312                                                                                 \
00313   i = 0;                                                                        \
00314 }
00315 
00316 
00317 #define MACRO_count__kg                                                         \
00318 {                                                                               \
00319   kg[ths->d - 1]++;                                                             \
00320   i = ths->d - 1;                                                               \
00321   while( ( kg[i] == ths->N[i] - 1) && ( i > 0))                                 \
00322   {                                                                             \
00323     kg[i - 1]++;                                                                \
00324     kg[i] = 0;                                                                  \
00325                                                                                 \
00326     i--;                                                                        \
00327   }                                                                             \
00328 }
00329 
00330 
00331 #define MACRO_update__c_phi_inv_k__lg_plain( which_one, which_phi)              \
00332 {                                                                               \
00333   for( t = i; t < ths->d; t++) {                                                \
00334     MACRO__c_phi_inv_k( which_phi);                                             \
00335     kg_plain[t+1] = kg_plain[t] * ths->n[t] + kg[t];                            \
00336   }                                                                             \
00337 }
00338 
00339 
00340 #define MACRO__c_phi_inv_k( which_phi)                                          \
00341 {                                                                               \
00342   c_phi_inv_k[t+1] = 0.5 * c_phi_inv_k[t] * MACRO_ ## which_phi;                \
00343 }
00344 
00345 
00346 #define MACRO_nfst_D(which_one)                                                 \
00347 static inline void nfst_D_ ## which_one (nfst_plan *ths)                               \
00348 {                                                                               \
00349   int k_L;                                    \
00350                                                                                 \
00351   int i, t;                                                                     \
00352   int kg[ths->d];                             \
00353   double c_phi_inv_k[ths->d+1];               \
00354   int kg_plain[ths->d+1];                     \
00355                                                                                 \
00356   double *g_hat, *f_hat;                      \
00357                                                                                 \
00358   g_hat = ths->g_hat;                                                           \
00359   f_hat = ths->f_hat;                                                           \
00360                                                                                 \
00361   MACRO_nfst_D_init_result_ ## which_one                                        \
00362                                                                                 \
00363   c_phi_inv_k[0] = 1;                                                           \
00364   kg_plain[0]    = 0;                                                           \
00365                                                                                 \
00366   MACRO_init__kg;                                                               \
00367                                                                                 \
00368   if( ths->nfst_flags & PRE_PHI_HUT)                                            \
00369                                                                                 \
00370     for( k_L = 0; k_L < ths->N_total; k_L++)                                    \
00371     {                                                                           \
00372       MACRO_update__c_phi_inv_k__lg_plain( which_one, with_PRE_PHI_HUT);        \
00373                                                                                 \
00374       MACRO_nfst_D_compute_ ## which_one;                                       \
00375                                                                                 \
00376       MACRO_count__kg;                                                          \
00377                                                                                 \
00378     } /* for(k_L) */                                                            \
00379                                                                                 \
00380   else                                                                          \
00381                                                                                 \
00382     for( k_L = 0; k_L < ths->N_total; k_L++)                                    \
00383     {                                                                           \
00384       MACRO_update__c_phi_inv_k__lg_plain( which_one, compute_PHI_HUT_INV);     \
00385                                                                                 \
00386       MACRO_nfst_D_compute_ ## which_one;                                       \
00387                                                                                 \
00388       MACRO_count__kg                                                           \
00389                                                                                 \
00390     } /* for(k_L) */                                                            \
00391 } /* nfst_D */
00392 
00393 MACRO_nfst_D(A)
00394 MACRO_nfst_D(T)
00395 
00396 
00397 
00398 
00399 
00400 
00401 
00405 #define MACRO_nfst_B_PRE_FULL_PSI_compute_A                                     \
00406 {                                                                               \
00407   (*fj) += ths->psi[ix] * g[ths->psi_index_g[ix]];                              \
00408 }
00409 
00410 #define MACRO_nfst_B_PRE_FULL_PSI_compute_T                                     \
00411 {                                                                               \
00412   g[ths->psi_index_g[ix]] += ths->psi[ix] * (*fj);                              \
00413 }
00414 
00415 
00416 
00417 #define MACRO_nfst_B_compute_A                                                  \
00418 {                                                                               \
00419   (*fj) += phi_tilde[ths->d] * g[lg_plain[ths->d]];                             \
00420 }
00421 
00422 #define MACRO_nfst_B_compute_T                                                  \
00423 {                                                                               \
00424   g[lg_plain[ths->d]] += phi_tilde[ths->d] * (*fj);                             \
00425 }
00426 
00427 
00428 
00429 #define MACRO_compute_lg_offset__count_lg( i0)                                  \
00430 {                                                                               \
00431   /* determine index in g-array corresponding to lb[(i0)] */                    \
00432   if( lb[(i0)] < 0)                                                             \
00433   {                                                                             \
00434     lg_offset[(i0)] =                                                           \
00435       (lb[(i0)] % nfst_fftw_2N( ths->n[(i0)])) + nfst_fftw_2N( ths->n[(i0)]);   \
00436   }                                                                             \
00437   else                                                                          \
00438   {                                                                             \
00439     lg_offset[(i0)] = lb[(i0)] % nfst_fftw_2N( ths->n[(i0)]);                   \
00440   }                                                                             \
00441                                                                                 \
00442   if( lg_offset[(i0)] > ths->n[(i0)]+1)                                         \
00443       lg_offset[(i0)] = -( nfst_fftw_2N( ths->n[(i0)]) - lg_offset[(i0)]);      \
00444 }
00445 
00446 
00447 
00448 #define MACRO_set__lg__to__lg_offset                                            \
00449 {                                                                               \
00450   if( lg_offset[i] <= 0)                                                        \
00451   {                                                                             \
00452     lg[i] = -lg_offset[i];                                                      \
00453     count_lg[i] = -1;                                                           \
00454   }                                                                             \
00455   else                                                                          \
00456   {                                                                             \
00457     lg[i] = +lg_offset[i];                                                      \
00458     count_lg[i] = +1;                                                           \
00459   }                                                                             \
00460 }
00461 
00462 
00463 
00464 #define MACRO_count__lg(dim)                                                    \
00465 {                                                                               \
00466   /* turn around when we hit one of the boundaries */                           \
00467   if( ((lg[(dim)] == 0) || (lg[(dim)] == (ths->n[(dim)] + 1))) )                \
00468     count_lg[(dim)] *= -1;                                                      \
00469                                                                                 \
00470   lg[(dim)] += count_lg[(dim)];                                                 \
00471 }
00472 
00473 
00474 
00475 #define MACRO_init_lb_lg_lc_phi_tilde_lg_plain( which_psi)                      \
00476 {                                                                               \
00477   for( i = 0; i < ths->d; i++)                                                  \
00478   {                                                                             \
00479     MACRO_nfst__lower_boundary( j, i);                                          \
00480                                                                                 \
00481     MACRO_compute_lg_offset__count_lg( i);                                      \
00482     MACRO_set__lg__to__lg_offset;                                               \
00483                                                                                 \
00484     /* counter for lg */                                                        \
00485     lc[i] = 0;                                                                  \
00486   }                                                                             \
00487                                                                                 \
00488   for( t = 0; t < ths->d; t++)                                                  \
00489   {                                                                             \
00490     if( lg[t] == 0)                                                             \
00491     {                                                                           \
00492       lg_plain[t+1]  = lg_plain[t] * ths->n[t];                                 \
00493       phi_tilde[t+1] = 0.0;                                                     \
00494     }                                                                           \
00495     else                                                                        \
00496       if( lg[t] == ths->n[t]+1)                                                 \
00497       {                                                                         \
00498         lg_plain[t+1]  = lg_plain[t] * ths->n[t] + ths->n[t]-1;                 \
00499         phi_tilde[t+1] = 0.0;                                                   \
00500       }                                                                         \
00501       else                                                                      \
00502       {                                                                         \
00503         MACRO__phi_tilde( which_psi);                                           \
00504         lg_plain[t+1]  = lg_plain[t] * ths->n[t] + lg[t]-1;                     \
00505       }                                                                         \
00506   }                                                                             \
00507                                                                                 \
00508   i = 0;                                                                        \
00509 }
00510 
00511 
00512 
00513 #define MACRO_count__lg_lc                                                      \
00514 {                                                                               \
00515   MACRO_count__lg( ths->d-1);                                                   \
00516                                                                                 \
00517   lc[ths->d - 1]++;                                                             \
00518   i = ths->d - 1;                                                               \
00519                                                                                 \
00520   while( (lc[i] == NFST_SUMMANDS) && (i > 0))                                   \
00521   {                                                                             \
00522     lc[i - 1]++;                                                                \
00523     lc[i] = 0;                                                                  \
00524                                                                                 \
00525     /* ansonsten lg[i-1] verschieben */                                         \
00526     MACRO_count__lg( i - 1);                                                    \
00527     /* lg[i] = anfangswert */                                                   \
00528     MACRO_set__lg__to__lg_offset;                                               \
00529                                                                                 \
00530     i--;                                                                        \
00531   }                                                                             \
00532 }
00533 
00534 
00535 #define MACRO_update__phi_tilde__lg_plain( which_psi)                           \
00536 {                                                                               \
00537   for( t = i; t < ths->d; t++)                                                  \
00538   {                                                                             \
00539     if( (lg[t] != 0) && (lg[t] != ths->n[t]+1))                                 \
00540     {                                                                           \
00541       MACRO__phi_tilde( which_psi);                                             \
00542       lg_plain[t+1] = lg_plain[t] * ths->n[t] + lg[t]-1;                        \
00543     }                                                                           \
00544     else                                                                        \
00545       phi_tilde[t+1] = 0.0;                                                     \
00546   }                                                                             \
00547 }
00548 
00549 
00550 
00551 #define MACRO__phi_tilde( which_psi)                                            \
00552 {                                                                               \
00553   phi_tilde[t+1] = (double)count_lg[t] * phi_tilde[t] * MACRO_ ## which_psi;    \
00554 }
00555 
00556 
00557 
00558 
00559 #define MACRO_nfst_B( which_one)                                                \
00560   static inline void nfst_B_ ## which_one ( nfst_plan *ths)                            \
00561   { /* MACRO_nfst_B */                                                          \
00562     int lb[ths->d];                       \
00563     int j, t, i;                          \
00564     int lprod, l_L, ix;                   \
00565     int lc[ths->d];                       \
00566     int lg[ths->d];                       \
00567     int lg_offset[ths->d];                \
00568     int count_lg[ths->d];                 \
00569     int lg_plain[ths->d+1];               \
00570     double *f, *g;                        \
00571     double phi_tilde[ths->d+1];           \
00572     double *fj;                           \
00573                                                                                 \
00574     f = ths->f; g = ths->g;                                                     \
00575                                                                                 \
00576     MACRO_nfst_B_init_result_ ## which_one                                      \
00577                                                                                 \
00578     /* both flags are set */                                                    \
00579     if( (ths->nfst_flags & PRE_PSI) && (ths->nfst_flags & PRE_FULL_PSI))        \
00580     {                                                                           \
00581       for( ix = 0, j = 0, fj = &f[0]; j < ths->M_total; j++, fj += 1)           \
00582         for( l_L = 0; l_L < ths->psi_index_f[j]; l_L++, ix++)                   \
00583         {                                                                       \
00584           MACRO_nfst_B_PRE_FULL_PSI_compute_ ## which_one;                      \
00585         }                                                                       \
00586     }                                                                           \
00587     else                                                                        \
00588     {                                                                           \
00589       phi_tilde[0] = 1;                                                         \
00590       lg_plain[0]  = 0;                                                         \
00591                                                                                 \
00592       for( t = 0, lprod = 1; t < ths->d; t++)                                   \
00593         lprod *= NFST_SUMMANDS;                                                 \
00594                                                                                 \
00595       /* PRE_PSI flag is set */                                                 \
00596       if( ths->nfst_flags & PRE_PSI)                                            \
00597       {                                                                         \
00598         for( j = 0, fj = &f[0]; j < ths->M_total; j++, fj += 1)                 \
00599         {                                                                       \
00600           MACRO_init_lb_lg_lc_phi_tilde_lg_plain( with_PRE_PSI);                \
00601                                                                                 \
00602           for( l_L = 0; l_L < lprod; l_L++)                                     \
00603           {                                                                     \
00604             MACRO_update__phi_tilde__lg_plain( with_PRE_PSI);                   \
00605                                                                                 \
00606             MACRO_nfst_B_compute_ ## which_one;                                 \
00607                                                                                 \
00608             MACRO_count__lg_lc;                                                 \
00609                                                                                 \
00610           } /* for( l_L) */                                                     \
00611         } /* for( j) */                                                         \
00612       } /* if( PRE_PSI) */                                                      \
00613                                                                                 \
00614       /* no PSI flag is set */                                                  \
00615       else                                                                      \
00616       {                                                                         \
00617         for( j = 0, fj = &f[0]; j < ths->M_total; j++, fj += 1)                 \
00618         {                                                                       \
00619           MACRO_init_lb_lg_lc_phi_tilde_lg_plain( compute_PSI);                 \
00620                                                                                 \
00621           for( l_L = 0; l_L < lprod; l_L++)                                     \
00622           {                                                                     \
00623             MACRO_update__phi_tilde__lg_plain( compute_PSI);                    \
00624                                                                                 \
00625             MACRO_nfst_B_compute_ ## which_one;                                 \
00626                                                                                 \
00627             MACRO_count__lg_lc;                                                 \
00628                                                                                 \
00629           } /* for(l_L) */                                                      \
00630         } /* for(j) */                                                          \
00631       } /* else(PRE_PSI) */                                                     \
00632     }/* else( PRE_PRE && FULL_PRE_PSI) */                                       \
00633 } /* nfst_B */
00634 
00635 MACRO_nfst_B(A)
00636 MACRO_nfst_B(T)
00637 
00638 
00639 
00640 
00641 
00642 
00647 void nfst_trafo( nfst_plan *ths)
00648 {
00653   ths->g_hat = ths->g1;
00654   ths->g     = ths->g2;
00655 
00656 
00662   TIC(0)
00663   nfst_D_A( ths);
00664   TOC(0)
00665 
00666 
00667   
00673   TIC(1)
00674   fftw_execute( ths->my_fftw_r2r_plan);
00675   TOC(1)
00676 
00677 
00678   
00683   TIC(2)
00684   nfst_B_A( ths);
00685   TOC(2)
00686 
00687 } /* nfst_trafo */
00688 
00689 
00690 
00691 
00692 void nfst_adjoint( nfst_plan *ths)
00693 {
00698   ths->g_hat = ths->g2;
00699   ths->g     = ths->g1;
00700 
00701 
00707   TIC(2)
00708   nfst_B_T( ths);
00709   TOC(2)
00710 
00711 
00718   TIC(1)
00719   fftw_execute( ths->my_fftw_r2r_plan);
00720   TOC(1)
00721 
00722 
00728   TIC(0)
00729   nfst_D_T( ths);
00730   TOC(0)
00731 
00732 } /* nfst_adjoint */
00733 
00734 
00735 
00740 void nfst_precompute_phi_hut( nfst_plan *ths)
00741 {
00742   int kg[ths->d];                      
00743   int t;                               
00745   ths->c_phi_inv = (double**)nfft_malloc( ths->d * sizeof( double*));
00746 
00747   for( t = 0; t < ths->d; t++)
00748   {
00749     ths->c_phi_inv[t] = (double*)nfft_malloc( ( ths->N[t] - 1) * sizeof( double));
00750 
00751     for( kg[t] = 0; kg[t] < ths->N[t] - 1; kg[t]++)
00752     {
00753       ths->c_phi_inv[t][kg[t]] = MACRO_compute_PHI_HUT_INV;
00754     }
00755   }
00756 } /* nfst_phi_hut */
00757 
00758 
00759 
00760 void nfst_precompute_psi( nfst_plan *ths)
00761 {
00762   int t;                                
00763   int j;                                
00764   int lc[ths->d];                       
00765   int lb[ths->d];                       
00767   for (t = 0; t < ths->d; t++)
00768   {
00769     for(j = 0; j < ths->M_total; j++)
00770     {
00771       MACRO_nfst__lower_boundary( j, t);
00772 
00773       for( lc[t] = 0; lc[t] < NFST_SUMMANDS; lc[t]++)
00774         ths->psi[(j * ths->d + t) * NFST_SUMMANDS + lc[t]] = MACRO_compute_PSI;
00775 
00776     } /* for(j) */
00777   }  /* for(t) */
00778 
00779   /* full precomputation of psi */
00780   if ( ths->nfst_flags & PRE_FULL_PSI)
00781     nfst_full_psi( ths, ths->nfst_full_psi_eps);
00782 
00783 } /* nfst_precompute_psi */
00784 
00785 
00786 
00788 void nfst_full_psi(nfst_plan *ths, double eps)
00789 {
00790   int t, i;                             
00791   int j;                                
00792   int l_L;                              
00793   int lc[ths->d];                       
00794   int lg_plain[ths->d+1];               
00795   int count_lg[ths->d];
00796   int lg_offset[ths->d];
00797   int lg[ths->d];
00798   int lprod;                            
00799   int lb[ths->d];                       
00801   double phi_tilde[ths->d+1];
00802 
00803   int *index_g, *index_f;
00804   double *new_psi;
00805   int ix, ix_old, size_psi;
00806 
00807   phi_tilde[0] = 1.0;
00808   lg_plain[0]  =   0;
00809 
00810   if(ths->nfst_flags & PRE_PSI)
00811   {
00812     size_psi = ths->M_total;
00813     index_f  =    (int*)nfft_malloc( ths->M_total  * sizeof( int));
00814     index_g  =    (int*)nfft_malloc( size_psi * sizeof( int));
00815     new_psi  = (double*)nfft_malloc( size_psi * sizeof( double));
00816 
00817     for( t = 0,lprod = 1; t < ths->d; t++)
00818     {
00819       lprod *= NFST_SUMMANDS;
00820       eps *= PHI( 0, t);
00821     }
00822 
00823     for( ix = 0, ix_old = 0, j = 0; j < ths->M_total; j++)
00824     {
00825       MACRO_init_lb_lg_lc_phi_tilde_lg_plain( with_PRE_PSI);
00826 
00827       for( l_L = 0; l_L < lprod; l_L++)
00828       {
00829         MACRO_update__phi_tilde__lg_plain( with_PRE_PSI);
00830 
00831         if( fabs(phi_tilde[ths->d]) > eps)
00832         {
00833           index_g[ix] =  lg_plain[ths->d];
00834           new_psi[ix] = phi_tilde[ths->d];
00835 
00836           ix++;
00837           if( ix == size_psi)
00838           {
00839             size_psi += ths->M_total;
00840             index_g   =    (int*)realloc( index_g, size_psi * sizeof( int));
00841             new_psi   = (double*)realloc( new_psi, size_psi * sizeof( double));
00842           }
00843         }
00844         MACRO_count__lg_lc;
00845 
00846       } /* for(l_L) */
00847 
00848       index_f[j] = ix - ix_old;
00849       ix_old     = ix;
00850 
00851     } /* for(j) */
00852 
00853     nfft_free( ths->psi);
00854 
00855     size_psi      = ix;
00856     ths->size_psi = size_psi;
00857     index_g       = (int*)realloc( index_g, size_psi * sizeof( int));
00858     new_psi       = (double*)realloc( new_psi, size_psi * sizeof( double));
00859 
00860     ths->psi         = new_psi;
00861     ths->psi_index_g = index_g;
00862     ths->psi_index_f = index_f;
00863 
00864   } /* if(PRE_PSI) */
00865 } /* nfst_full_psi */
00866 
00867 
00868 
00869 
00870 void nfst_init_help( nfst_plan *ths)
00871 {
00872   int t;                                
00874   ths->N_total = nfst_prod_minus_a_int( ths->N, 1, ths->d);
00875 
00876   ths->sigma   = (double*)nfft_malloc( ths->d * sizeof( double));
00877 
00878   for( t = 0; t < ths->d; t++)
00879     /* FIXME: n/N or (n+1)/N */
00880     ths->sigma[t] = ((double)ths->n[t] + 1) / ths->N[t];
00881 
00882   /* assign r2r transform kinds for each dimension */
00883   ths->r2r_kind = (fftw_r2r_kind*) nfft_malloc ( ths->d * sizeof( fftw_r2r_kind));
00884   for (t = 0; t < ths->d; t++)
00885     ths->r2r_kind[t] = FFTW_RODFT00;
00886 
00887 
00888   WINDOW_HELP_INIT;
00889 
00890   if(ths->nfst_flags & MALLOC_X)
00891     ths->x = (double*)nfft_malloc( ths->d * ths->M_total * sizeof( double));
00892 
00893   if(ths->nfst_flags & MALLOC_F_HAT)
00894     ths->f_hat = (double*)nfft_malloc( ths->N_total * sizeof( double));
00895 
00896   if(ths->nfst_flags & MALLOC_F)
00897     ths->f = (double*)nfft_malloc( ths->M_total * sizeof( double));
00898 
00899   if(ths->nfst_flags & PRE_PHI_HUT)
00900     nfst_precompute_phi_hut( ths);
00901 
00902   /* NO FFTW_MALLOC HERE */
00903   if(ths->nfst_flags & PRE_PSI)
00904   {
00905     ths->psi =
00906       (double*)nfft_malloc( ths->M_total * ths->d * NFST_SUMMANDS * sizeof( double));
00907 
00911     ths->nfst_full_psi_eps = pow(10, -10);
00912   }
00913 
00914   if(ths->nfst_flags & FFTW_INIT)
00915   {
00916       ths->g1 =
00917         (double*)nfft_malloc( nfst_prod_minus_a_int( ths->n, 0, ths->d) * sizeof( double));
00918 
00919       if(ths->nfst_flags & FFT_OUT_OF_PLACE)
00920         ths->g2 =
00921           (double*)nfft_malloc( nfst_prod_minus_a_int( ths->n, 0, ths->d) * sizeof( double));
00922       else
00923         ths->g2 = ths->g1;
00924 
00925       ths->my_fftw_r2r_plan =
00926         fftw_plan_r2r( ths->d, ths->n, ths->g1, ths->g2, ths->r2r_kind, ths->fftw_flags);
00927   }
00928 
00929   ths->mv_trafo = (void (*) (void* ))nfst_trafo;
00930   ths->mv_adjoint = (void (*) (void* ))nfst_adjoint;
00931 }
00932 
00933 void nfst_init( nfst_plan *ths, int d, int *N, int M_total)
00934 {
00935   int t;                                
00937   ths->d = d;
00938 
00939   ths->N      = (int*)nfft_malloc( ths->d * sizeof( int));
00940 
00941   for(t = 0;t < d; t++)
00942     ths->N[t] = N[t];
00943 
00944   ths->n      = (int*)nfft_malloc( ths->d * sizeof( int));
00945 
00946   for( t = 0; t < d; t++)
00947     ths->n[t] = 2 * X(next_power_of_2)( ths->N[t]) - 1;
00948 
00949   ths->M_total = M_total;
00950 
00951 /* Was soll dieser Ausdruck machen? Es handelt sich um eine Ganzzahl!
00952 
00953   WINDOW_HELP_ESTIMATE_m;
00954 */  
00955 
00956   ths->nfst_flags = NFST_DEFAULT_FLAGS;
00957   ths->fftw_flags = FFTW_DEFAULT_FLAGS;
00958 
00959   nfst_init_help( ths);
00960 }
00961 
00962 
00963 void nfst_init_m( nfst_plan *ths, int d, int *N, int M_total, int m)
00964 {
00965   int t, n[d];
00966 
00967   for( t = 0; t < d; t++)
00968     n[t] = nfst_fftw_2N( X(next_power_of_2)( N[t]));
00969 
00970   nfst_init_guru( ths, d, N, M_total, n, m, NFST_DEFAULT_FLAGS, FFTW_DEFAULT_FLAGS);
00971 }
00972 
00973 
00974 void nfst_init_guru( nfst_plan *ths, int d, int *N,
00975      int M_total, int *n, int m,
00976      unsigned nfst_flags, unsigned fftw_flags)
00977 {
00978   int t;             
00980   ths->d = d;
00981   ths->M_total = M_total;
00982 
00983   ths->N      = (int*)nfft_malloc( ths->d * sizeof( int));
00984 
00985   for( t = 0; t < d; t++)
00986     ths->N[t]      = N[t];
00987 
00988   ths->n      = (int*)nfft_malloc( ths->d * sizeof( int));
00989 
00990   for( t = 0; t < d; t++)
00991     ths->n[t]      = n[t];
00992 
00993   ths->m = m;
00994 
00995   ths->nfst_flags = nfst_flags;
00996   ths->fftw_flags = fftw_flags;
00997 
00998   nfst_init_help( ths);
00999 }
01000 
01001 
01002 void nfst_init_1d( nfst_plan *ths, int N0, int M_total)
01003 {
01004   int N[1];
01005 
01006   N[0] = N0;
01007   nfst_init( ths, 1, N, M_total);
01008 }
01009 
01010 void nfst_init_2d( nfst_plan *ths, int N0, int N1, int M_total)
01011 {
01012   int N[2];
01013 
01014   N[0] = N0;
01015   N[1] = N1;
01016   nfst_init( ths, 2, N, M_total);
01017 }
01018 
01019 void nfst_init_3d( nfst_plan *ths, int N0, int N1, int N2, int M_total)
01020 {
01021   int N[3];
01022 
01023   N[0] = N0;
01024   N[1] = N1;
01025   N[2] = N2;
01026   nfst_init( ths, 3, N, M_total);
01027 }
01028 
01029 void nfst_finalize( nfst_plan *ths)
01030 {
01031   int t; /* index over dimensions */
01032 
01033   if( ths->nfst_flags & FFTW_INIT)
01034   {
01035     fftw_destroy_plan( ths->my_fftw_r2r_plan);
01036 
01037     if( ths->nfst_flags & FFT_OUT_OF_PLACE)
01038       nfft_free( ths->g2);
01039 
01040     nfft_free( ths->g1);
01041   }
01042 
01043   /* NO FFTW_FREE HERE */
01044   if( ths->nfst_flags & PRE_PSI)
01045   {
01046     if( ths->nfst_flags & PRE_FULL_PSI)
01047     {
01048       nfft_free( ths->psi_index_g);
01049       nfft_free( ths->psi_index_f);
01050     }
01051 
01052     nfft_free( ths->psi);
01053   }
01054 
01055   if( ths->nfst_flags & PRE_PHI_HUT) {
01056     for( t = 0; t < ths->d; t++)
01057       nfft_free( ths->c_phi_inv[t]);
01058     nfft_free( ths->c_phi_inv);
01059   }
01060 
01061   if( ths->nfst_flags & MALLOC_F)
01062     nfft_free( ths->f);
01063 
01064   if( ths->nfst_flags & MALLOC_F_HAT)
01065     nfft_free( ths->f_hat);
01066 
01067   if( ths->nfst_flags & MALLOC_X)
01068     nfft_free( ths->x);
01069 
01070   WINDOW_HELP_FINALIZE;
01071 
01072   nfft_free( ths->N);
01073   nfft_free( ths->n);
01074   nfft_free( ths->sigma);
01075 
01076   nfft_free(ths->r2r_kind);
01077 } /* nfst_finalize */
01078 

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409