NFFT Logo 3.2.3
nfst.c
1 /*
2  * Copyright (c) 2002, 2012 Jens Keiner, Stefan Kunis, Daniel Potts
3  *
4  * This program is free software; you can redistribute it and/or modify it under
5  * the terms of the GNU General Public License as published by the Free Software
6  * Foundation; either version 2 of the License, or (at your option) any later
7  * version.
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12  * details.
13  *
14  * You should have received a copy of the GNU General Public License along with
15  * this program; if not, write to the Free Software Foundation, Inc., 51
16  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  */
18 
19 /* $Id: nfst.c 3896 2012-10-10 12:19:26Z tovo $ */
20 
27 #include <stdio.h>
28 #include <math.h>
29 #include <string.h>
30 #include <stdlib.h>
31 
32 #include "nfft3util.h"
33 #include "nfft3.h"
34 #include "infft.h"
35 
39 #define NFST_DEFAULT_FLAGS PRE_PHI_HUT|\
40  PRE_PSI|\
41  MALLOC_X|\
42  MALLOC_F_HAT|\
43  MALLOC_F|\
44  FFTW_INIT|\
45  FFT_OUT_OF_PLACE
46 
47 #define FFTW_DEFAULT_FLAGS FFTW_ESTIMATE|\
48  FFTW_DESTROY_INPUT
49 
50 #define NFST_SUMMANDS ( 2 * ths->m + 2)
51 #define NODE(p,r) ( ths->x[(p) * ths->d + (r)])
52 
53 #define MACRO_ndst_init_result_trafo \
54  memset( f, 0, ths->M_total * sizeof( double));
55 #define MACRO_ndst_init_result_adjoint \
56  memset( f_hat, 0, ths->N_total * sizeof( double));
57 
58 
59 #define MACRO_nfst_D_init_result_A \
60  memset(g_hat, 0, nfst_prod_minus_a_int( ths->n, 0, ths->d) * sizeof( double));
61 #define MACRO_nfst_D_init_result_T \
62  memset(f_hat, 0, ths->N_total * sizeof( double));
63 
64 #define MACRO_nfst_B_init_result_A \
65  memset(f, 0, ths->M_total * sizeof( double));
66 #define MACRO_nfst_B_init_result_T \
67  memset(g, 0, nfst_prod_minus_a_int( ths->n, 0, ths->d) * sizeof( double));
68 
69 
70 #define NFST_PRE_WINFUN( d) ths->N[d] = 2 * ths->N[d]; \
71  ths->n[d] = nfst_fftw_2N( ths->n[d]);
72 
73 #define NFST_POST_WINFUN( d) ths->N[d] = (LRINT(0.5 * ths->N[d])); \
74  ths->n[d] = nfst_fftw_2N_rev( ths->n[d]);
75 
76 
77 #define NFST_WINDOW_HELP_INIT WINDOW_HELP_INIT
78 
79 
80 double nfst_phi_hut( nfst_plan *ths, int k, int d)
81 {
82  NFST_PRE_WINFUN( d);
83  double phi_hut_tmp = PHI_HUT( k, d);
84  NFST_POST_WINFUN( d);
85 
86  return phi_hut_tmp;
87 }
88 
89 double nfst_phi( nfst_plan *ths, double x, int d)
90 {
91  NFST_PRE_WINFUN( d);
92  double phi_tmp = PHI( x, d);
93  NFST_POST_WINFUN( d);
94 
95  return phi_tmp;
96 }
97 
98 int nfst_fftw_2N( int n)
99 {
100  return 2 * ( n + 1);
101 }
102 
103 int nfst_fftw_2N_rev( int n)
104 {
105  div_t n_div;
106 
107  n_div = div(n, 2);
108  return n_div.quot - 1;
109 }
110 
111 #define MACRO_with_sin_vec sin_vec[t][ka[t]]
112 #define MACRO_without_sin_vec sin( 2.0 * PI * (ka[t]+1) * NODE(j,t))
113 
114 
115 #define MACRO_with_PRE_PHI_HUT ths->c_phi_inv[t][kg[t]];
116 #define MACRO_compute_PHI_HUT_INV (1.0 / (nfst_phi_hut( ths, kg[t]+1, t)))
117 
118 #define MACRO_with_PRE_PSI ths->psi[(j * ths->d + t) * NFST_SUMMANDS + lc[t]];
119 #define MACRO_compute_PSI \
120  nfst_phi( ths, NODE(j,t) - (( double)(lc[t] + lb[t])) / nfst_fftw_2N( ths->n[t]), t)
121 
122 
123 
139 #define MACRO_ndst_malloc__sin_vec \
140  \
141  double **sin_vec; \
142  sin_vec = (double**)nfft_malloc( ths->d * sizeof( double*)); \
143  for( t = 0; t < ths->d; t++) \
144  sin_vec[t] = (double*)nfft_malloc( ( ths->N[t] - 1) * sizeof( double)); \
145 
146 
147 
148 
149 #define MACRO_ndst_free__sin_vec \
150 { \
151  /* free allocated memory */ \
152  for( t = 0; t < ths->d; t++) \
153  nfft_free( sin_vec[t]); \
154  nfft_free( sin_vec); \
155 }
156 
157 
158 
159 #define MACRO_ndst_init__sin_vec \
160 { \
161  for( t = 0; t < ths->d; t++) \
162  { \
163  cos_x[t] = cos( 2.0 * PI * NODE(j,t)); \
164  sin_vec[t][0] = sin( 2.0 * PI * NODE(j,t)); \
165  sin_vec[t][1] = sin( 4.0 * PI * NODE(j,t)); \
166  for( k = 2; k < ths->N[t] - 1; k++) \
167  sin_vec[t][k] = 2.0 * cos_x[t] * sin_vec[t][k-1] \
168  - sin_vec[t][k-2]; \
169  } \
170 }
171 
172 
173 #define MACRO_ndst_init__k__sin_k( which_one) \
174 { \
175  sin_k[0] = 1.0; \
176  for( t = 0; t < ths->d; t++) \
177  ka[t] = 0; \
178  \
179  for( t = 0; t < ths->d; t++) \
180  { \
181  sin_k[t+1] = sin_k[t] * MACRO_ ##which_one; \
182  } \
183 }
184 
185 
186 #define MACRO_ndst_count__k__sin_k( which_one) \
187 { \
188  ka[ths->d-1]++; \
189  i = ths->d - 1; \
190  while( ( ka[i] == ths->N[i] - 1) && ( i > 0)) \
191  { \
192  ka[i - 1]++; \
193  ka[i] = 0; \
194  \
195  i--; \
196  } \
197  for( t = i; t < ths->d; t++) \
198  sin_k[t+1] = sin_k[t] * MACRO_ ##which_one; \
199 }
200 
201 
202 #define MACRO_ndst_compute__trafo \
203 { \
204  f[j] += f_hat[k] * sin_k[ths->d]; \
205 }
206 
207 #define MACRO_ndst_compute__adjoint \
208 { \
209  f_hat[k] += f[j] * sin_k[ths->d]; \
210 }
211 
212 
213 /* slow (trafo) transform */
214 #define MACRO_ndst( which_one) \
215  void nfst_ ## which_one ## _direct ( nfst_plan *ths) \
216  { \
217  int j, k, t, i; \
218  int ka[ths->d]; \
219  double sin_k[ths->d+1]; \
220  double cos_x[ths->d]; \
221  \
222  double *f = ths->f; \
223  double *f_hat = ths->f_hat; \
224  \
225  MACRO_ndst_init_result_ ## which_one; \
226  \
227  if( ths->d == 1) \
228  for( j = 0; j < ths->M_total; j++) \
229  for( k = 0; k < ths->N_total; k++) \
230  { \
231  sin_k[ths->d] = sin( 2.0 * PI * (k+1) * NODE(j,0)); \
232  MACRO_ndst_compute__ ## which_one; \
233  } \
234  else \
235  if( 1 == 0) /*FIXME: remove slow slow ... */ \
236  /* slow ndst */ \
237  for( j = 0; j < ths->M_total; j++) \
238  { \
239  MACRO_ndst_init__k__sin_k(without_sin_vec); \
240  \
241  for( k = 0; k < ths->N_total; k++) \
242  { \
243  MACRO_ndst_compute__ ## which_one; \
244  \
245  MACRO_ndst_count__k__sin_k(without_sin_vec); \
246  } \
247  } \
248  else \
249  { \
250  /* fast nfst_trafo_direct */ \
251  MACRO_ndst_malloc__sin_vec; \
252  \
253  for( j = 0; j < ths->M_total; j++) \
254  { \
255  MACRO_ndst_init__sin_vec; \
256  \
257  MACRO_ndst_init__k__sin_k(with_sin_vec); \
258  \
259  for( k = 0; k < ths->N_total; k++) \
260  { \
261  MACRO_ndst_compute__ ## which_one; \
262  \
263  MACRO_ndst_count__k__sin_k(with_sin_vec); \
264  } \
265  } \
266  MACRO_ndst_free__sin_vec; \
267  } \
268  } /* ndst_{trafo, adjoint} */
269 
270 
271 MACRO_ndst(trafo)
272 MACRO_ndst(adjoint)
273 
274 
275 
276 
291 #define MACRO_nfst__lower_boundary( j,act_dim) \
292 { \
293  lb[(act_dim)] = \
294  (LRINT(NODE((j),(act_dim)) * nfst_fftw_2N( ths->n[(act_dim)]))) - ths->m; \
295 }
296 
297 #define MACRO_nfst_D_compute_A \
298 { \
299  g_hat[kg_plain[ths->d]] = f_hat[k_L] * c_phi_inv_k[ths->d]; \
300 }
301 
302 #define MACRO_nfst_D_compute_T \
303 { \
304  f_hat[k_L] = g_hat[kg_plain[ths->d]] * c_phi_inv_k[ths->d]; \
305 }
306 
307 
308 #define MACRO_init__kg \
309 { \
310  for( t = 0; t < ths->d; t++) \
311  kg[t] = 0; \
312  \
313  i = 0; \
314 }
315 
316 
317 #define MACRO_count__kg \
318 { \
319  kg[ths->d - 1]++; \
320  i = ths->d - 1; \
321  while( ( kg[i] == ths->N[i] - 1) && ( i > 0)) \
322  { \
323  kg[i - 1]++; \
324  kg[i] = 0; \
325  \
326  i--; \
327  } \
328 }
329 
330 
331 #define MACRO_update__c_phi_inv_k__lg_plain( which_one, which_phi) \
332 { \
333  for( t = i; t < ths->d; t++) { \
334  MACRO__c_phi_inv_k( which_phi); \
335  kg_plain[t+1] = kg_plain[t] * ths->n[t] + kg[t]; \
336  } \
337 }
338 
339 
340 #define MACRO__c_phi_inv_k( which_phi) \
341 { \
342  c_phi_inv_k[t+1] = 0.5 * c_phi_inv_k[t] * MACRO_ ## which_phi; \
343 }
344 
345 
346 #define MACRO_nfst_D(which_one) \
347 static inline void nfst_D_ ## which_one (nfst_plan *ths) \
348 { \
349  int k_L; \
350  \
351  int i, t; \
352  int kg[ths->d]; \
353  double c_phi_inv_k[ths->d+1]; \
354  int kg_plain[ths->d+1]; \
355  \
356  double *g_hat, *f_hat; \
357  \
358  g_hat = ths->g_hat; \
359  f_hat = ths->f_hat; \
360  \
361  MACRO_nfst_D_init_result_ ## which_one \
362  \
363  c_phi_inv_k[0] = 1; \
364  kg_plain[0] = 0; \
365  \
366  MACRO_init__kg; \
367  \
368  if( ths->nfst_flags & PRE_PHI_HUT) \
369  \
370  for( k_L = 0; k_L < ths->N_total; k_L++) \
371  { \
372  MACRO_update__c_phi_inv_k__lg_plain( which_one, with_PRE_PHI_HUT); \
373  \
374  MACRO_nfst_D_compute_ ## which_one; \
375  \
376  MACRO_count__kg; \
377  \
378  } /* for(k_L) */ \
379  \
380  else \
381  \
382  for( k_L = 0; k_L < ths->N_total; k_L++) \
383  { \
384  MACRO_update__c_phi_inv_k__lg_plain( which_one, compute_PHI_HUT_INV); \
385  \
386  MACRO_nfst_D_compute_ ## which_one; \
387  \
388  MACRO_count__kg \
389  \
390  } /* for(k_L) */ \
391 } /* nfst_D */
392 
393 MACRO_nfst_D(A)
394 MACRO_nfst_D(T)
395 
396 
397 
398 
399 
400 
401 
405 #define MACRO_nfst_B_PRE_FULL_PSI_compute_A \
406 { \
407  (*fj) += ths->psi[ix] * g[ths->psi_index_g[ix]]; \
408 }
409 
410 #define MACRO_nfst_B_PRE_FULL_PSI_compute_T \
411 { \
412  g[ths->psi_index_g[ix]] += ths->psi[ix] * (*fj); \
413 }
414 
415 
416 
417 #define MACRO_nfst_B_compute_A \
418 { \
419  (*fj) += phi_tilde[ths->d] * g[lg_plain[ths->d]]; \
420 }
421 
422 #define MACRO_nfst_B_compute_T \
423 { \
424  g[lg_plain[ths->d]] += phi_tilde[ths->d] * (*fj); \
425 }
426 
427 
428 
429 #define MACRO_compute_lg_offset__count_lg( i0) \
430 { \
431  /* determine index in g-array corresponding to lb[(i0)] */ \
432  if( lb[(i0)] < 0) \
433  { \
434  lg_offset[(i0)] = \
435  (lb[(i0)] % nfst_fftw_2N( ths->n[(i0)])) + nfst_fftw_2N( ths->n[(i0)]); \
436  } \
437  else \
438  { \
439  lg_offset[(i0)] = lb[(i0)] % nfst_fftw_2N( ths->n[(i0)]); \
440  } \
441  \
442  if( lg_offset[(i0)] > ths->n[(i0)]+1) \
443  lg_offset[(i0)] = -( nfst_fftw_2N( ths->n[(i0)]) - lg_offset[(i0)]); \
444 }
445 
446 
447 
448 #define MACRO_set__lg__to__lg_offset \
449 { \
450  if( lg_offset[i] <= 0) \
451  { \
452  lg[i] = -lg_offset[i]; \
453  count_lg[i] = -1; \
454  } \
455  else \
456  { \
457  lg[i] = +lg_offset[i]; \
458  count_lg[i] = +1; \
459  } \
460 }
461 
462 
463 
464 #define MACRO_count__lg(dim) \
465 { \
466  /* turn around when we hit one of the boundaries */ \
467  if( ((lg[(dim)] == 0) || (lg[(dim)] == (ths->n[(dim)] + 1))) ) \
468  count_lg[(dim)] *= -1; \
469  \
470  lg[(dim)] += count_lg[(dim)]; \
471 }
472 
473 
474 
475 #define MACRO_init_lb_lg_lc_phi_tilde_lg_plain( which_psi) \
476 { \
477  for( i = 0; i < ths->d; i++) \
478  { \
479  MACRO_nfst__lower_boundary( j, i); \
480  \
481  MACRO_compute_lg_offset__count_lg( i); \
482  MACRO_set__lg__to__lg_offset; \
483  \
484  /* counter for lg */ \
485  lc[i] = 0; \
486  } \
487  \
488  for( t = 0; t < ths->d; t++) \
489  { \
490  if( lg[t] == 0) \
491  { \
492  lg_plain[t+1] = lg_plain[t] * ths->n[t]; \
493  phi_tilde[t+1] = 0.0; \
494  } \
495  else \
496  if( lg[t] == ths->n[t]+1) \
497  { \
498  lg_plain[t+1] = lg_plain[t] * ths->n[t] + ths->n[t]-1; \
499  phi_tilde[t+1] = 0.0; \
500  } \
501  else \
502  { \
503  MACRO__phi_tilde( which_psi); \
504  lg_plain[t+1] = lg_plain[t] * ths->n[t] + lg[t]-1; \
505  } \
506  } \
507  \
508  i = 0; \
509 }
510 
511 
512 
513 #define MACRO_count__lg_lc \
514 { \
515  MACRO_count__lg( ths->d-1); \
516  \
517  lc[ths->d - 1]++; \
518  i = ths->d - 1; \
519  \
520  while( (lc[i] == NFST_SUMMANDS) && (i > 0)) \
521  { \
522  lc[i - 1]++; \
523  lc[i] = 0; \
524  \
525  /* ansonsten lg[i-1] verschieben */ \
526  MACRO_count__lg( i - 1); \
527  /* lg[i] = anfangswert */ \
528  MACRO_set__lg__to__lg_offset; \
529  \
530  i--; \
531  } \
532 }
533 
534 
535 #define MACRO_update__phi_tilde__lg_plain( which_psi) \
536 { \
537  for( t = i; t < ths->d; t++) \
538  { \
539  if( (lg[t] != 0) && (lg[t] != ths->n[t]+1)) \
540  { \
541  MACRO__phi_tilde( which_psi); \
542  lg_plain[t+1] = lg_plain[t] * ths->n[t] + lg[t]-1; \
543  } \
544  else \
545  phi_tilde[t+1] = 0.0; \
546  } \
547 }
548 
549 
550 
551 #define MACRO__phi_tilde( which_psi) \
552 { \
553  phi_tilde[t+1] = (double)count_lg[t] * phi_tilde[t] * MACRO_ ## which_psi; \
554 }
555 
556 
557 
558 
559 #define MACRO_nfst_B( which_one) \
560  static inline void nfst_B_ ## which_one ( nfst_plan *ths) \
561  { /* MACRO_nfst_B */ \
562  int lb[ths->d]; \
563  int j, t, i; \
564  int lprod, l_L, ix; \
565  int lc[ths->d]; \
566  int lg[ths->d]; \
567  int lg_offset[ths->d]; \
568  int count_lg[ths->d]; \
569  int lg_plain[ths->d+1]; \
570  double *f, *g; \
571  double phi_tilde[ths->d+1]; \
572  double *fj; \
573  \
574  f = ths->f; g = ths->g; \
575  \
576  MACRO_nfst_B_init_result_ ## which_one \
577  \
578  /* both flags are set */ \
579  if( (ths->nfst_flags & PRE_PSI) && (ths->nfst_flags & PRE_FULL_PSI)) \
580  { \
581  for( ix = 0, j = 0, fj = &f[0]; j < ths->M_total; j++, fj += 1) \
582  for( l_L = 0; l_L < ths->psi_index_f[j]; l_L++, ix++) \
583  { \
584  MACRO_nfst_B_PRE_FULL_PSI_compute_ ## which_one; \
585  } \
586  } \
587  else \
588  { \
589  phi_tilde[0] = 1; \
590  lg_plain[0] = 0; \
591  \
592  for( t = 0, lprod = 1; t < ths->d; t++) \
593  lprod *= NFST_SUMMANDS; \
594  \
595  /* PRE_PSI flag is set */ \
596  if( ths->nfst_flags & PRE_PSI) \
597  { \
598  for( j = 0, fj = &f[0]; j < ths->M_total; j++, fj += 1) \
599  { \
600  MACRO_init_lb_lg_lc_phi_tilde_lg_plain( with_PRE_PSI); \
601  \
602  for( l_L = 0; l_L < lprod; l_L++) \
603  { \
604  MACRO_update__phi_tilde__lg_plain( with_PRE_PSI); \
605  \
606  MACRO_nfst_B_compute_ ## which_one; \
607  \
608  MACRO_count__lg_lc; \
609  \
610  } /* for( l_L) */ \
611  } /* for( j) */ \
612  } /* if( PRE_PSI) */ \
613  \
614  /* no PSI flag is set */ \
615  else \
616  { \
617  for( j = 0, fj = &f[0]; j < ths->M_total; j++, fj += 1) \
618  { \
619  MACRO_init_lb_lg_lc_phi_tilde_lg_plain( compute_PSI); \
620  \
621  for( l_L = 0; l_L < lprod; l_L++) \
622  { \
623  MACRO_update__phi_tilde__lg_plain( compute_PSI); \
624  \
625  MACRO_nfst_B_compute_ ## which_one; \
626  \
627  MACRO_count__lg_lc; \
628  \
629  } /* for(l_L) */ \
630  } /* for(j) */ \
631  } /* else(PRE_PSI) */ \
632  }/* else( PRE_PRE && FULL_PRE_PSI) */ \
633 } /* nfst_B */
634 
635 MACRO_nfst_B(A)
636 MACRO_nfst_B(T)
637 
638 
639 
640 
641 
642 
647 void nfst_trafo( nfst_plan *ths)
648 {
653  ths->g_hat = ths->g1;
654  ths->g = ths->g2;
655 
656 
662  TIC(0)
663  nfst_D_A( ths);
664  TOC(0)
665 
666 
667 
673  TIC(1)
674  fftw_execute( ths->my_fftw_r2r_plan);
675  TOC(1)
676 
677 
678 
683  TIC(2)
684  nfst_B_A( ths);
685  TOC(2)
686 
687 } /* nfst_trafo */
688 
689 
690 
691 
692 void nfst_adjoint( nfst_plan *ths)
693 {
698  ths->g_hat = ths->g2;
699  ths->g = ths->g1;
700 
701 
707  TIC(2)
708  nfst_B_T( ths);
709  TOC(2)
710 
711 
718  TIC(1)
719  fftw_execute( ths->my_fftw_r2r_plan);
720  TOC(1)
721 
722 
728  TIC(0)
729  nfst_D_T( ths);
730  TOC(0)
731 
732 } /* nfst_adjoint */
733 
734 
735 
740 void nfst_precompute_phi_hut( nfst_plan *ths)
741 {
742  int kg[ths->d];
743  int t;
745  ths->c_phi_inv = (double**)nfft_malloc( ths->d * sizeof( double*));
746 
747  for( t = 0; t < ths->d; t++)
748  {
749  ths->c_phi_inv[t] = (double*)nfft_malloc( ( ths->N[t] - 1) * sizeof( double));
750 
751  for( kg[t] = 0; kg[t] < ths->N[t] - 1; kg[t]++)
752  {
753  ths->c_phi_inv[t][kg[t]] = MACRO_compute_PHI_HUT_INV;
754  }
755  }
756 } /* nfst_phi_hut */
757 
758 
759 
760 void nfst_precompute_psi( nfst_plan *ths)
761 {
762  int t;
763  int j;
764  int lc[ths->d];
765  int lb[ths->d];
767  for (t = 0; t < ths->d; t++)
768  {
769  for(j = 0; j < ths->M_total; j++)
770  {
771  MACRO_nfst__lower_boundary( j, t);
772 
773  for( lc[t] = 0; lc[t] < NFST_SUMMANDS; lc[t]++)
774  ths->psi[(j * ths->d + t) * NFST_SUMMANDS + lc[t]] = MACRO_compute_PSI;
775 
776  } /* for(j) */
777  } /* for(t) */
778 
779  /* full precomputation of psi */
780  if ( ths->nfst_flags & PRE_FULL_PSI)
781  nfst_full_psi( ths, ths->nfst_full_psi_eps);
782 
783 } /* nfst_precompute_psi */
784 
785 
786 
788 void nfst_full_psi(nfst_plan *ths, double eps)
789 {
790  int t, i;
791  int j;
792  int l_L;
793  int lc[ths->d];
794  int lg_plain[ths->d+1];
795  int count_lg[ths->d];
796  int lg_offset[ths->d];
797  int lg[ths->d];
798  int lprod;
799  int lb[ths->d];
801  double phi_tilde[ths->d+1];
802 
803  int *index_g, *index_f;
804  double *new_psi;
805  int ix, ix_old, size_psi;
806 
807  phi_tilde[0] = 1.0;
808  lg_plain[0] = 0;
809 
810  if(ths->nfst_flags & PRE_PSI)
811  {
812  size_psi = ths->M_total;
813  index_f = (int*)nfft_malloc( ths->M_total * sizeof( int));
814  index_g = (int*)nfft_malloc( size_psi * sizeof( int));
815  new_psi = (double*)nfft_malloc( size_psi * sizeof( double));
816 
817  for( t = 0,lprod = 1; t < ths->d; t++)
818  {
819  lprod *= NFST_SUMMANDS;
820  eps *= PHI( 0, t);
821  }
822 
823  for( ix = 0, ix_old = 0, j = 0; j < ths->M_total; j++)
824  {
825  MACRO_init_lb_lg_lc_phi_tilde_lg_plain( with_PRE_PSI);
826 
827  for( l_L = 0; l_L < lprod; l_L++)
828  {
829  MACRO_update__phi_tilde__lg_plain( with_PRE_PSI);
830 
831  if( fabs(phi_tilde[ths->d]) > eps)
832  {
833  index_g[ix] = lg_plain[ths->d];
834  new_psi[ix] = phi_tilde[ths->d];
835 
836  ix++;
837  if( ix == size_psi)
838  {
839  size_psi += ths->M_total;
840  index_g = (int*)realloc( index_g, size_psi * sizeof( int));
841  new_psi = (double*)realloc( new_psi, size_psi * sizeof( double));
842  }
843  }
844  MACRO_count__lg_lc;
845 
846  } /* for(l_L) */
847 
848  index_f[j] = ix - ix_old;
849  ix_old = ix;
850 
851  } /* for(j) */
852 
853  nfft_free( ths->psi);
854 
855  size_psi = ix;
856  ths->size_psi = size_psi;
857  index_g = (int*)realloc( index_g, size_psi * sizeof( int));
858  new_psi = (double*)realloc( new_psi, size_psi * sizeof( double));
859 
860  ths->psi = new_psi;
861  ths->psi_index_g = index_g;
862  ths->psi_index_f = index_f;
863 
864  } /* if(PRE_PSI) */
865 } /* nfst_full_psi */
866 
867 
868 
869 
870 void nfst_init_help( nfst_plan *ths)
871 {
872  int t;
874  ths->N_total = nfst_prod_minus_a_int( ths->N, 1, ths->d);
875 
876  ths->sigma = (double*)nfft_malloc( ths->d * sizeof( double));
877 
878  for( t = 0; t < ths->d; t++)
879  /* FIXME: n/N or (n+1)/N */
880  ths->sigma[t] = ((double)ths->n[t] + 1) / ths->N[t];
881 
882  /* assign r2r transform kinds for each dimension */
883  ths->r2r_kind = (fftw_r2r_kind*) nfft_malloc ( ths->d * sizeof( fftw_r2r_kind));
884  for (t = 0; t < ths->d; t++)
885  ths->r2r_kind[t] = FFTW_RODFT00;
886 
887 
888  WINDOW_HELP_INIT;
889 
890  if(ths->nfst_flags & MALLOC_X)
891  ths->x = (double*)nfft_malloc( ths->d * ths->M_total * sizeof( double));
892 
893  if(ths->nfst_flags & MALLOC_F_HAT)
894  ths->f_hat = (double*)nfft_malloc( ths->N_total * sizeof( double));
895 
896  if(ths->nfst_flags & MALLOC_F)
897  ths->f = (double*)nfft_malloc( ths->M_total * sizeof( double));
898 
899  if(ths->nfst_flags & PRE_PHI_HUT)
900  nfst_precompute_phi_hut( ths);
901 
902  /* NO FFTW_MALLOC HERE */
903  if(ths->nfst_flags & PRE_PSI)
904  {
905  ths->psi =
906  (double*)nfft_malloc( ths->M_total * ths->d * NFST_SUMMANDS * sizeof( double));
907 
911  ths->nfst_full_psi_eps = pow(10, -10);
912  }
913 
914  if(ths->nfst_flags & FFTW_INIT)
915  {
916  ths->g1 =
917  (double*)nfft_malloc( nfst_prod_minus_a_int( ths->n, 0, ths->d) * sizeof( double));
918 
919  if(ths->nfst_flags & FFT_OUT_OF_PLACE)
920  ths->g2 =
921  (double*)nfft_malloc( nfst_prod_minus_a_int( ths->n, 0, ths->d) * sizeof( double));
922  else
923  ths->g2 = ths->g1;
924 
925  ths->my_fftw_r2r_plan =
926  fftw_plan_r2r( ths->d, ths->n, ths->g1, ths->g2, ths->r2r_kind, ths->fftw_flags);
927  }
928 
929  ths->mv_trafo = (void (*) (void* ))nfst_trafo;
930  ths->mv_adjoint = (void (*) (void* ))nfst_adjoint;
931 }
932 
933 void nfst_init( nfst_plan *ths, int d, int *N, int M_total)
934 {
935  int t;
937  ths->d = d;
938 
939  ths->N = (int*)nfft_malloc( ths->d * sizeof( int));
940 
941  for(t = 0;t < d; t++)
942  ths->N[t] = N[t];
943 
944  ths->n = (int*)nfft_malloc( ths->d * sizeof( int));
945 
946  for( t = 0; t < d; t++)
947  ths->n[t] = 2 * X(next_power_of_2)( ths->N[t]) - 1;
948 
949  ths->M_total = M_total;
950 
951 /* Was soll dieser Ausdruck machen? Es handelt sich um eine Ganzzahl!
952 
953  WINDOW_HELP_ESTIMATE_m;
954 */
955 
956  ths->nfst_flags = NFST_DEFAULT_FLAGS;
957  ths->fftw_flags = FFTW_DEFAULT_FLAGS;
958 
959  nfst_init_help( ths);
960 }
961 
962 
963 void nfst_init_m( nfst_plan *ths, int d, int *N, int M_total, int m)
964 {
965  int t, n[d];
966 
967  for( t = 0; t < d; t++)
968  n[t] = nfst_fftw_2N( X(next_power_of_2)( N[t]));
969 
970  nfst_init_guru( ths, d, N, M_total, n, m, NFST_DEFAULT_FLAGS, FFTW_DEFAULT_FLAGS);
971 }
972 
973 
974 void nfst_init_guru( nfst_plan *ths, int d, int *N,
975  int M_total, int *n, int m,
976  unsigned nfst_flags, unsigned fftw_flags)
977 {
978  int t;
980  ths->d = d;
981  ths->M_total = M_total;
982 
983  ths->N = (int*)nfft_malloc( ths->d * sizeof( int));
984 
985  for( t = 0; t < d; t++)
986  ths->N[t] = N[t];
987 
988  ths->n = (int*)nfft_malloc( ths->d * sizeof( int));
989 
990  for( t = 0; t < d; t++)
991  ths->n[t] = n[t];
992 
993  ths->m = m;
994 
995  ths->nfst_flags = nfst_flags;
996  ths->fftw_flags = fftw_flags;
997 
998  nfst_init_help( ths);
999 }
1000 
1001 
1002 void nfst_init_1d( nfst_plan *ths, int N0, int M_total)
1003 {
1004  int N[1];
1005 
1006  N[0] = N0;
1007  nfst_init( ths, 1, N, M_total);
1008 }
1009 
1010 void nfst_init_2d( nfst_plan *ths, int N0, int N1, int M_total)
1011 {
1012  int N[2];
1013 
1014  N[0] = N0;
1015  N[1] = N1;
1016  nfst_init( ths, 2, N, M_total);
1017 }
1018 
1019 void nfst_init_3d( nfst_plan *ths, int N0, int N1, int N2, int M_total)
1020 {
1021  int N[3];
1022 
1023  N[0] = N0;
1024  N[1] = N1;
1025  N[2] = N2;
1026  nfst_init( ths, 3, N, M_total);
1027 }
1028 
1029 void nfst_finalize( nfst_plan *ths)
1030 {
1031  int t; /* index over dimensions */
1032 
1033  if( ths->nfst_flags & FFTW_INIT)
1034  {
1035  fftw_destroy_plan( ths->my_fftw_r2r_plan);
1036 
1037  if( ths->nfst_flags & FFT_OUT_OF_PLACE)
1038  nfft_free( ths->g2);
1039 
1040  nfft_free( ths->g1);
1041  }
1042 
1043  /* NO FFTW_FREE HERE */
1044  if( ths->nfst_flags & PRE_PSI)
1045  {
1046  if( ths->nfst_flags & PRE_FULL_PSI)
1047  {
1048  nfft_free( ths->psi_index_g);
1049  nfft_free( ths->psi_index_f);
1050  }
1051 
1052  nfft_free( ths->psi);
1053  }
1054 
1055  if( ths->nfst_flags & PRE_PHI_HUT) {
1056  for( t = 0; t < ths->d; t++)
1057  nfft_free( ths->c_phi_inv[t]);
1058  nfft_free( ths->c_phi_inv);
1059  }
1060 
1061  if( ths->nfst_flags & MALLOC_F)
1062  nfft_free( ths->f);
1063 
1064  if( ths->nfst_flags & MALLOC_F_HAT)
1065  nfft_free( ths->f_hat);
1066 
1067  if( ths->nfst_flags & MALLOC_X)
1068  nfft_free( ths->x);
1069 
1070  WINDOW_HELP_FINALIZE;
1071 
1072  nfft_free( ths->N);
1073  nfft_free( ths->n);
1074  nfft_free( ths->sigma);
1075 
1076  nfft_free(ths->r2r_kind);
1077 } /* nfst_finalize */
1078 
Generated on Tue Apr 30 2013 by Doxygen 1.8.1