NFFT Logo 3.2.2
nfct.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: nfct.c 3896 2012-10-10 12:19:26Z tovo $ */
00020 
00021 /* Nonequispaced fast cosine transform
00022  * Author: Steffen Klatt 2004-2006, Jens Keiner 2010 */
00023 
00024 #include <stdio.h>
00025 #include <math.h>
00026 #include <string.h>
00027 #include <stdlib.h>
00028 
00029 #include "nfft3util.h"
00030 #include "nfft3.h"
00031 #include "infft.h"
00032 
00033 #undef X
00034 #if defined(NFFT_SINGLE)
00035 #define X(name) CONCAT(nfctf_, name)
00036 #elif defined(NFFT_LDOUBLE)
00037 #define X(name) CONCAT(nfctl_, name)
00038 #else
00039 #define X(name) CONCAT(nfct_, name)
00040 #endif
00041 
00042 static inline int fftw_2N(int n)
00043 {
00044   return 2 * (n - 1);
00045 }
00046 
00047 static inline int fftw_2N_rev(int n)
00048 {
00049   return (LRINT(K(0.5) * n) + 1);
00050 }
00051 
00052 static inline int prod_int(int *vec, int d)
00053 {
00054   int t, prod = 1;
00055 
00056   for (t = 0; t < d; t++)
00057     prod *= vec[t];
00058 
00059   return prod;
00060 }
00061 
00062 /* handy shortcuts */
00063 #define NFCT_DEFAULT_FLAGS PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | \
00064   MALLOC_F | FFTW_INIT | FFT_OUT_OF_PLACE
00065 
00066 #define FFTW_DEFAULT_FLAGS FFTW_ESTIMATE | FFTW_DESTROY_INPUT
00067 
00068 #define NFCT_SUMMANDS (2 * ths->m + 2)
00069 #define NODE(p,r) (ths->x[(p) * ths->d + (r)])
00070 
00071 #define MACRO_ndct_init_result_trafo \
00072   memset(f, 0, ths->M_total * sizeof(R));
00073 #define MACRO_ndct_init_result_adjoint \
00074   memset(f_hat, 0, ths->N_total * sizeof(R));
00075 
00076 #define MACRO_nfct_D_init_result_A \
00077   memset(g_hat, 0, prod_int(ths->n, ths->d) * sizeof(R));
00078 #define MACRO_nfct_D_init_result_T \
00079   memset(f_hat, 0, ths->N_total * sizeof(R));
00080 
00081 #define MACRO_nfct_B_init_result_A \
00082   memset(f, 0, ths->M_total * sizeof(R));
00083 #define MACRO_nfct_B_init_result_T \
00084   memset(g, 0, prod_int(ths->n, ths->d) * sizeof(R));
00085 
00086 #define NFCT_PRE_WINFUN(d) ths->N[d] = 2 * ths->N[d]; \
00087   ths->n[d] = fftw_2N(ths->n[d]);
00088 
00089 #define NFCT_POST_WINFUN(d) ths->N[d] = LRINT(K(0.5) * ths->N[d]); \
00090   ths->n[d] = fftw_2N_rev(ths->n[d]);
00091 
00092 #define NFCT_WINDOW_HELP_INIT WINDOW_HELP_INIT
00093 
00094 R X(phi_hut)(X(plan) *ths, int k, int d)
00095 {
00096   NFCT_PRE_WINFUN(d);
00097   R phi_hut_tmp = PHI_HUT(k, d);
00098   NFCT_POST_WINFUN(d);
00099 
00100   return phi_hut_tmp;
00101 }
00102 
00103 R X(phi)(X(plan) *ths, R x, int d)
00104 {
00105   NFCT_PRE_WINFUN(d);
00106   R phi_tmp = PHI(x, d);
00107   NFCT_POST_WINFUN(d);
00108 
00109   return phi_tmp;
00110 }
00111 
00112 #define MACRO_with_cos_vec cos_vec[t][ka[t]]
00113 #define MACRO_without_cos_vec COS(K(2.0) * KPI * ka[t] * NODE(j,t))
00114 
00115 #define MACRO_with_PRE_PHI_HUT ths->c_phi_inv[t][kg[t]];
00116 #define MACRO_compute_PHI_HUT_INV (K(1.0) / (X(phi_hut)(ths, kg[t], t)))
00117 
00118 #define MACRO_with_PRE_PSI ths->psi[(j * ths->d + t) * NFCT_SUMMANDS + lc[t]]
00119 #define MACRO_compute_PSI X(phi)(ths, NODE(j,t) - ((R)(lc[t] + lb[t])) / (K(2.0)*((R)(ths->n[t])-K(1.0))/*(R)(fftw_2N(ths->n[t]))*/), t)
00120 
00136 #define MACRO_ndct_malloc__cos_vec \
00137   R **cos_vec; \
00138   cos_vec = (R**)Y(malloc)(ths->d * sizeof(R*)); \
00139   for (t = 0; t < ths->d; t++) \
00140     cos_vec[t] = (R*)Y(malloc)(ths->N[t] * sizeof(R));
00141 
00142 #define MACRO_ndct_free__cos_vec \
00143 { \
00144   /* free allocated memory */ \
00145   for (t = 0; t < ths->d; t++) \
00146     Y(free)(cos_vec[t]); \
00147   Y(free)(cos_vec); \
00148 }
00149 
00150 #define MACRO_ndct_init__cos_vec \
00151 { \
00152   for(t = 0; t < ths->d; t++) \
00153   { \
00154     cos_vec[t][0] = K(1.0); \
00155     cos_vec[t][1] = COS(K(2.0) * KPI * NODE(j,t)); \
00156     for (k = 2; k < ths->N[t]; k++) \
00157       cos_vec[t][k] = K(2.0) * cos_vec[t][1] * cos_vec[t][k-1] - \
00158         cos_vec[t][k-2]; \
00159   } \
00160 }
00161 
00162 #define MACRO_ndct_init__k__cos_k(which_one) \
00163 { \
00164   cos_k[0] = K(1.0); \
00165   for (t = 0; t < ths->d; t++) \
00166     ka[t] = 0; \
00167 \
00168   for (t = 0; t < ths->d; t++) \
00169   { \
00170     cos_k[t+1] = cos_k[t] * MACRO_ ##which_one; \
00171   } \
00172 }
00173 
00174 #define MACRO_ndct_count__k__cos_k(which_one) \
00175 { \
00176   ka[ths->d-1]++; \
00177   i = ths->d - 1; \
00178   while ((ka[i] == ths->N[i]) && (i > 0)) \
00179   { \
00180     ka[i - 1]++; \
00181     ka[i] = 0; \
00182     i--; \
00183   } \
00184   for (t = i; t < ths->d; t++) \
00185     cos_k[t+1] = cos_k[t] * MACRO_ ##which_one; \
00186 }
00187 
00188 #define MACRO_ndct_compute__trafo \
00189 { \
00190   f[j] += f_hat[k] * cos_k[ths->d]; \
00191 }
00192 
00193 #define MACRO_ndct_compute__adjoint \
00194 { \
00195   f_hat[k] += f[j] * cos_k[ths->d]; \
00196 }
00197 
00198 /* slow (trafo) transform */
00199 #define MACRO_ndct(which_one) \
00200   void X(which_one ## _direct) (X(plan) *ths) \
00201   { \
00202     int j, k, t, i; \
00203     int ka[ths->d]; \
00204     R cos_k[ths->d+1]; \
00205     R *f = ths->f; \
00206     R *f_hat = ths->f_hat; \
00207 \
00208     MACRO_ndct_init_result_ ## which_one; \
00209     if (ths->d == 1) \
00210       for (j = 0; j < ths->M_total; j++) \
00211       { \
00212         for (k = 0; k < ths->N[0]; k++) \
00213         { \
00214           cos_k[ths->d] = COS(K(2.0) * KPI * k * NODE(j,0)); \
00215           MACRO_ndct_compute__ ## which_one; \
00216         } \
00217       } \
00218     else \
00219     { \
00220       /* fast nfct_trafo_direct */ \
00221       MACRO_ndct_malloc__cos_vec; \
00222 \
00223       for (j = 0; j < ths->M_total; j++) \
00224       { \
00225         MACRO_ndct_init__cos_vec; \
00226         MACRO_ndct_init__k__cos_k(with_cos_vec); \
00227 \
00228         for (k = 0; k < ths->N_total; k++) \
00229         { \
00230           MACRO_ndct_compute__ ## which_one; \
00231           MACRO_ndct_count__k__cos_k(with_cos_vec); \
00232         } \
00233       } \
00234       MACRO_ndct_free__cos_vec; \
00235     } \
00236 } /* ndct_{trafo, adjoint} */
00237 
00238 MACRO_ndct(trafo)
00239 MACRO_ndct(adjoint)
00240 
00255 #define MACRO_nfct__lower_boundary(j,act_dim) \
00256 { \
00257   lb[(act_dim)] = \
00258     LRINT(NODE((j),(act_dim)) * fftw_2N(ths->n[(act_dim)])) - ths->m; \
00259 }
00260 
00261 #define MACRO_nfct_D_compute_A \
00262 { \
00263   g_hat[kg_plain[ths->d]] = f_hat[k_L] * c_phi_inv_k[ths->d]; \
00264 }
00265 
00266 #define MACRO_nfct_D_compute_T \
00267 { \
00268   f_hat[k_L] = g_hat[kg_plain[ths->d]] * c_phi_inv_k[ths->d]; \
00269 }
00270 
00271 #define MACRO_init__kg \
00272 { \
00273   for (t = 0; t < ths->d; t++) \
00274     kg[t] = 0; \
00275   i = 0; \
00276 }
00277 
00278 #define MACRO_count__kg \
00279 { \
00280 \
00281   kg[ths->d - 1]++; \
00282   i = ths->d - 1; \
00283   while ((kg[i] == ths->N[i]) && (i > 0)) \
00284   { \
00285     kg[i - 1]++; \
00286     kg[i] = 0; \
00287     i--; \
00288   } \
00289 }
00290 
00291 #define MACRO_update__phi_inv_k__kg_plain(what_kind, which_phi) \
00292 { \
00293   for (t = i; t < ths->d; t++) \
00294   { \
00295     MACRO__c_phi_inv_k__ ## what_kind(which_phi); \
00296     kg_plain[t+1] = kg_plain[t] * ths->n[t] + kg[t]; \
00297   } \
00298 }
00299 
00300 #define MACRO__c_phi_inv_k__A(which_phi) \
00301 { \
00302   if (kg[t] == 0) \
00303   { \
00304     c_phi_inv_k[t+1] = c_phi_inv_k[t] * MACRO_ ## which_phi; \
00305   } \
00306   else \
00307   { \
00308     c_phi_inv_k[t+1] = K(0.5) * c_phi_inv_k[t] * MACRO_ ## which_phi; \
00309   } \
00310 }
00311 
00312 #define MACRO__c_phi_inv_k__T(which_phi) \
00313 { \
00314   c_phi_inv_k[t+1] = c_phi_inv_k[t] * MACRO_ ## which_phi; \
00315 }
00316 
00317 #define MACRO_nfct_D(which_one) \
00318 static inline void D_ ## which_one (X(plan) *ths) \
00319 { \
00320   int k_L; /* plain index */ \
00321   int i, t; \
00322   int kg[ths->d]; /* multi index in g_hat,c_phi */ \
00323   R c_phi_inv_k[ths->d+1]; /* postfix product of PHI_HUT_INV */ \
00324   int kg_plain[ths->d+1]; /* postfix plain index */ \
00325   R *g_hat, *f_hat; /* local copy */ \
00326 \
00327   g_hat = ths->g_hat; \
00328   f_hat = ths->f_hat; \
00329 \
00330   MACRO_nfct_D_init_result_ ## which_one \
00331 \
00332   c_phi_inv_k[0] = K(1.0); \
00333   kg_plain[0] = 0; \
00334 \
00335   MACRO_init__kg; \
00336 \
00337   if (ths->nfct_flags & PRE_PHI_HUT) \
00338     for (k_L = 0; k_L < ths->N_total; k_L++) \
00339     { \
00340       MACRO_update__phi_inv_k__kg_plain(which_one, with_PRE_PHI_HUT); \
00341       MACRO_nfct_D_compute_ ## which_one; \
00342       MACRO_count__kg; \
00343     } /* for (k_L) */ \
00344   else \
00345     for (k_L = 0; k_L < ths->N_total; k_L++) \
00346     { \
00347       MACRO_update__phi_inv_k__kg_plain(which_one,compute_PHI_HUT_INV); \
00348       MACRO_nfct_D_compute_ ## which_one; \
00349       MACRO_count__kg; \
00350     } /* for(k_L) */ \
00351 } /* nfct_D */
00352 
00353 MACRO_nfct_D(A)
00354 MACRO_nfct_D(T)
00355 
00359 #define MACRO_nfct_B_PRE_FULL_PSI_compute_A \
00360 { \
00361   (*fj) += ths->psi[ix] * g[ths->psi_index_g[ix]]; \
00362 }
00363 
00364 #define MACRO_nfct_B_PRE_FULL_PSI_compute_T \
00365 { \
00366   g[ths->psi_index_g[ix]] += ths->psi[ix] * (*fj); \
00367 }
00368 
00369 #define MACRO_nfct_B_compute_A \
00370 { \
00371   (*fj) += phi_tilde[ths->d] * g[lg_plain[ths->d]]; \
00372 }
00373 
00374 #define MACRO_nfct_B_compute_T \
00375 { \
00376   g[lg_plain[ths->d]] += phi_tilde[ths->d] * (*fj); \
00377 }
00378 
00379 #define MACRO_compute_lg_offset__count_lg(i0) \
00380 { \
00381   /* determine index in g-array corresponding to lb[(i0)] */ \
00382   if (lb[(i0)] < 0) \
00383     lg_offset[(i0)] = \
00384       (lb[(i0)] % fftw_2N(ths->n[(i0)])) + fftw_2N(ths->n[(i0)]); \
00385   else \
00386     lg_offset[(i0)] = lb[(i0)] % (fftw_2N(ths->n[(i0)])); \
00387     if (lg_offset[(i0)] >= ths->n[(i0)]) \
00388       lg_offset[(i0)] = -(fftw_2N(ths->n[(i0)]) - lg_offset[(i0)]); \
00389 }
00390 
00391 #define MACRO_set__lg__to__lg_offset \
00392 { \
00393   if (lg_offset[i] <= 0) \
00394   { \
00395     lg[i] = -lg_offset[i]; \
00396     count_lg[i] = -1; \
00397   } \
00398   else \
00399   { \
00400     lg[i] = +lg_offset[i]; \
00401     count_lg[i] = +1; \
00402   } \
00403 }
00404 
00405 #define MACRO_count__lg(dim) \
00406 { \
00407   /* turn around if we hit one of the boundaries */ \
00408   if ((lg[(dim)] == 0) || (lg[(dim)] == ths->n[(dim)]-1)) \
00409     count_lg[(dim)] *= -1; \
00410   /* move array index */ \
00411   lg[(dim)] += count_lg[(dim)]; \
00412 }
00413 
00414 #define MACRO_init_lb_lg_lc \
00415 { \
00416   for (i = 0; i < ths->d; i++) \
00417   { \
00418     MACRO_nfct__lower_boundary(j, i); \
00419     MACRO_compute_lg_offset__count_lg(i); \
00420     MACRO_set__lg__to__lg_offset; \
00421     /* counter for lg */ \
00422     lc[i] = 0; \
00423    } \
00424    i = 0; \
00425 }
00426 
00427 #define MACRO_count__lg_lc \
00428 { \
00429   MACRO_count__lg(ths->d-1); \
00430   lc[ths->d - 1]++; \
00431   i = ths->d - 1; \
00432   while ((lc[i] == NFCT_SUMMANDS) && (i > 0)) \
00433   { \
00434     lc[i - 1]++; \
00435     lc[i] = 0; \
00436     /* ansonsten lg[i-1] verschieben */ \
00437     MACRO_count__lg(i - 1); \
00438     /* lg[i] = anfangswert */ \
00439     MACRO_set__lg__to__lg_offset; \
00440     i--; \
00441   } \
00442 }
00443 
00444 #define  MACRO_update_phi_tilde_lg_plain(which_one, which_psi) \
00445 { \
00446   for (t = i; t < ths->d; t++) \
00447   { \
00448     MACRO__phi_tilde__ ## which_one(which_psi); \
00449     lg_plain[t+1]  = lg_plain[t]  * ths->n[t] + lg[t]; \
00450   } \
00451 }
00452 
00453 #define MACRO__phi_tilde__A(which_psi) \
00454 { \
00455   phi_tilde[t+1] = phi_tilde[t] * MACRO_ ## which_psi; \
00456 }
00457 
00458 #define MACRO__phi_tilde__T(which_psi) \
00459 { \
00460   if(lg[t] == 0 || lg[t] == ths->n[t] - 1) \
00461   { \
00462     phi_tilde[t+1] = phi_tilde[t] * MACRO_ ## which_psi; \
00463   } \
00464   else \
00465   { \
00466     phi_tilde[t+1] = K(0.5) * phi_tilde[t] * MACRO_ ## which_psi; \
00467   } \
00468 }
00469 
00470 #define MACRO_nfct_B(which_one) \
00471 static inline void B_ ## which_one (nfct_plan *ths) \
00472 { /* MACRO_nfct_B */ \
00473   int lb[ths->d]; /* multi band with respect to x_j */ \
00474   int j, t, i; /* index nodes, help vars */ \
00475   int lprod, l_L, ix; /* index one row of B */ \
00476   int lc[ths->d]; /* multi index 0<=lc<2m+2 */ \
00477   int lg[ths->d]; /* real index of g in array */ \
00478   int lg_offset[ths->d]; /* offset in g according to u */ \
00479   int count_lg[ths->d]; /* count summands (2m+2) */ \
00480   int lg_plain[ths->d+1]; /* index of g in multi_array */ \
00481   R *f, *g; /* local copy */ \
00482   R phi_tilde[ths->d+1]; /* holds values for psi */ \
00483   R *fj; /* pointer to final result */ \
00484 \
00485   f = ths->f; g = ths->g; \
00486 \
00487   MACRO_nfct_B_init_result_ ## which_one \
00488 \
00489   /* both flags are set */ \
00490   if ((ths->nfct_flags & PRE_PSI)&&(ths->nfct_flags & PRE_FULL_PSI)) \
00491   { \
00492     for (ix = 0, j = 0, fj = &f[0]; j < ths->M_total; j++, fj += 1) \
00493       for (l_L = 0; l_L < ths->psi_index_f[j]; l_L++, ix++) \
00494       { \
00495         MACRO_nfct_B_PRE_FULL_PSI_compute_ ## which_one; \
00496       } \
00497   } \
00498   else \
00499   { \
00500     phi_tilde[0] = K(1.0); \
00501     lg_plain[0]  = 0; \
00502 \
00503     for (t = 0, lprod = 1; t < ths->d; t++) \
00504       lprod *= NFCT_SUMMANDS; \
00505 \
00506     /* PRE_PSI flag is set */ \
00507     if (ths->nfct_flags & PRE_PSI) \
00508       for (j = 0, fj = &f[0]; j < ths->M_total; j++, fj += 1) \
00509         { \
00510           MACRO_init_lb_lg_lc; \
00511           for (l_L = 0; l_L < lprod; l_L++) \
00512           { \
00513             MACRO_update_phi_tilde_lg_plain(which_one, with_PRE_PSI); \
00514             MACRO_nfct_B_compute_ ## which_one; \
00515             MACRO_count__lg_lc; \
00516           } /* for(l_L) */ \
00517         } /* for(j) */ \
00518 \
00519     /* no PSI flag is set */ \
00520     else \
00521       for (j = 0, fj = &f[0]; j < ths->M_total; j++, fj += 1) \
00522       { \
00523         MACRO_init_lb_lg_lc; \
00524         for (l_L = 0; l_L < lprod; l_L++) \
00525         { \
00526           MACRO_update_phi_tilde_lg_plain(which_one,compute_PSI); \
00527           MACRO_nfct_B_compute_ ## which_one; \
00528           MACRO_count__lg_lc; \
00529         } /* for (l_L) */ \
00530       } /* for (j) */ \
00531   } /* else (PRE_PSI && FULL_PRE_PSI) */ \
00532 } /* nfct_B */
00533 
00534 MACRO_nfct_B(A)
00535 MACRO_nfct_B(T)
00536 
00537 /* more memory, but faster */
00538 #define MACRO_nfct_full_psi(which_one) \
00539 static inline void full_psi__ ## which_one(nfct_plan *ths) \
00540 { \
00541   int t, i; /* index over all dimensions */ \
00542   int j; /* node index */ \
00543   int l_L; /* plain index 0 <= l_L < lprod */ \
00544   int lc[ths->d]; /* multi index 0<=lj<u+o+1 */ \
00545   int lg_plain[ths->d+1]; /* postfix plain index */ \
00546   int count_lg[ths->d]; \
00547   int lg_offset[ths->d]; \
00548   int lg[ths->d]; \
00549   int lprod; /* 'bandwidth' of matrix B */ \
00550   int lb[ths->d]; /* depends on x_j */ \
00551 \
00552   R phi_tilde[ths->d+1]; \
00553   R eps = ths->nfct_full_psi_eps; \
00554 \
00555   int *index_g, *index_f; \
00556   R *new_psi; \
00557   int ix, ix_old, size_psi; \
00558 \
00559   phi_tilde[0] = K(1.0); \
00560   lg_plain[0]  =   0; \
00561  \
00562   if (ths->nfct_flags & PRE_PSI) \
00563   { \
00564     size_psi = ths->M_total; \
00565     index_f = (int*)Y(malloc)(ths->M_total  * sizeof(int)); \
00566     index_g = (int*)Y(malloc)(size_psi * sizeof(int)); \
00567     new_psi = (R*)Y(malloc)(size_psi * sizeof(R)); \
00568 \
00569     for (t = 0,lprod = 1; t < ths->d; t++) \
00570     { \
00571       lprod *= NFCT_SUMMANDS; \
00572       eps *= nfct_phi(ths, K(0.0), t); \
00573     } \
00574 \
00575     for (ix = 0, ix_old = 0, j = 0; j < ths->M_total; j++) \
00576     { \
00577       MACRO_init_lb_lg_lc; \
00578 \
00579       for (l_L = 0; l_L < lprod; l_L++) \
00580       { \
00581         MACRO_update_phi_tilde_lg_plain(which_one, with_PRE_PSI); \
00582 \
00583         if (phi_tilde[ths->d] > eps) \
00584         { \
00585           index_g[ix] = lg_plain[ths->d]; \
00586           new_psi[ix] = phi_tilde[ths->d]; \
00587 \
00588           ix++; \
00589           if (ix == size_psi) \
00590           { \
00591             size_psi += ths->M_total; \
00592             index_g = (int*)realloc(index_g, size_psi * sizeof(int)); \
00593             new_psi = (R*)realloc(new_psi, size_psi * sizeof(R)); \
00594           } \
00595         } \
00596 \
00597         MACRO_count__lg_lc; \
00598 \
00599       } /* for (l_L) */ \
00600 \
00601       index_f[j] = ix - ix_old; \
00602       ix_old = ix; \
00603 \
00604     } /* for(j) */ \
00605 \
00606     Y(free)(ths->psi); \
00607     size_psi = ix; \
00608     ths->size_psi = size_psi; \
00609 \
00610     index_g = (int*)realloc(index_g, size_psi * sizeof(int)); \
00611     new_psi = (R*)realloc(new_psi, size_psi * sizeof(R)); \
00612 \
00613     ths->psi = new_psi; \
00614     ths->psi_index_g = index_g; \
00615     ths->psi_index_f = index_f; \
00616 \
00617   } /* if(PRE_PSI) */ \
00618 }
00619 
00620 MACRO_nfct_full_psi(A)
00621 MACRO_nfct_full_psi(T)
00622 
00623 /* user routines */
00624 
00625 void X(trafo)(X(plan) *ths)
00626 {
00627   /* use ths->my_fftw_r2r_plan */
00628   ths->g_hat = ths->g1;
00629   ths->g = ths->g2;
00630 
00631   /* form \f$ \hat g_k = \frac{\hat f_k}{c_k\left(\phi\right)} \text{ for }
00632    * k \in I_N \f$ */
00633   TIC(0)
00634   D_A(ths);
00635   TOC(0)
00636 
00637   /* Compute by d-variate discrete Fourier transform
00638    * \f$ g_l = \sum_{k \in I_N} \hat g_k {\rm e}^{-2\pi {\rm i} \frac{kl}{n}}
00639    * \text{ for } l \in I_n \f$ */
00640   TIC(1)
00641   Z(execute)(ths->my_fftw_r2r_plan);
00642   TOC(1)
00643 
00644   if (ths->nfct_flags & PRE_FULL_PSI)
00645     full_psi__A(ths);
00646 
00647   /* Set \f$ f_j = \sum_{l \in I_n,m(x_j)} g_l \psi\left(x_j-\frac{l}{n}\right)
00648    * \text{ for } j=0,\hdots,M-1 \f$ */
00649   TIC(2)
00650   B_A(ths);
00651   TOC(2)
00652 
00653   if (ths->nfct_flags & PRE_FULL_PSI)
00654   {
00655     Y(free)(ths->psi_index_g);
00656     Y(free)(ths->psi_index_f);
00657   }
00658 } /* nfct_trafo */
00659 
00660 void X(adjoint)(X(plan) *ths)
00661 {
00662   /* use ths->my_fftw_plan */
00663   ths->g_hat = ths->g2;
00664   ths->g = ths->g1;
00665 
00666   if (ths->nfct_flags & PRE_FULL_PSI)
00667     full_psi__T(ths);
00668 
00669   /* Set \f$ g_l = \sum_{j=0}^{M-1} f_j \psi\left(x_j-\frac{l}{n}\right)
00670    * \text{ for } l \in I_n,m(x_j) \f$ */
00671   TIC(2)
00672   B_T(ths);
00673   TOC(2)
00674 
00675   if (ths->nfct_flags & PRE_FULL_PSI)
00676   {
00677     Y(free)(ths->psi_index_g);
00678     Y(free)(ths->psi_index_f);
00679   }
00680 
00681   /* Compute by d-variate discrete cosine transform
00682    * \f$ \hat g_k = \sum_{l \in I_n} g_l {\rm e}^{-2\pi {\rm i} \frac{kl}{n}}
00683    * \text{ for }  k \in I_N\f$ */
00684   TIC(1)
00685   Z(execute)(ths->my_fftw_r2r_plan);
00686   TOC(1)
00687 
00688   /* Form \f$ \hat f_k = \frac{\hat g_k}{c_k\left(\phi\right)} \text{ for }
00689    * k \in I_N \f$ */
00690   TIC(0)
00691   D_T(ths);
00692   TOC(0)
00693 
00694 } /* nfct_adjoint */
00695 
00696 static inline void precompute_phi_hut(X(plan) *ths)
00697 {
00698   int kg[ths->d]; /* index over all frequencies */
00699   int t; /* index over all dimensions */
00700 
00701   ths->c_phi_inv = (R**)Y(malloc)(ths->d * sizeof(R*));
00702 
00703   for (t = 0; t < ths->d; t++)
00704   {
00705     ths->c_phi_inv[t] = (R*)Y(malloc)(ths->N[t] * sizeof(R));
00706 
00707     for (kg[t] = 0; kg[t] < ths->N[t]; kg[t]++)
00708     {
00709       ths->c_phi_inv[t][kg[t]] = MACRO_compute_PHI_HUT_INV;
00710     }
00711   }
00712 } /* nfct_phi_hut */
00713 
00714 void X(precompute_psi)(X(plan) *ths)
00715 {
00716   int t; /* index over all dimensions */
00717   int j; /* index over all nodes */
00718   int lc[ths->d]; /* index 0<=lj<u+o+1 */
00719   int lb[ths->d]; /* depends on x_j */
00720 
00721   for (t = 0; t < ths->d; t++)
00722   {
00723     for (j = 0; j < ths->M_total; j++)
00724     {
00725       MACRO_nfct__lower_boundary(j, t);
00726       for(lc[t] = 0; lc[t] < NFCT_SUMMANDS; lc[t]++)
00727         ths->psi[(j * ths->d + t) * NFCT_SUMMANDS + lc[t]] = MACRO_compute_PSI;
00728     } /* for (j) */
00729   } /* for (t) */
00730 } /* nfct_precompute_psi */
00731 
00732 static inline void init_help(X(plan) *ths)
00733 {
00734   int t; /* index over all dimensions */
00735 
00736   ths->N_total = prod_int(ths->N, ths->d);
00737   ths->sigma = (R*)Y(malloc)(ths->d * sizeof(R));
00738 
00739   for (t = 0; t < ths->d; t++)
00740     ths->sigma[t] = ((R)(ths->n[t] - 1)) / ths->N[t];
00741 
00742   /* Assign r2r transform kinds for each dimension */
00743   ths->r2r_kind = (Z(r2r_kind)*)Y(malloc)(ths->d * sizeof (Z(r2r_kind)));
00744   for (t = 0; t < ths->d; t++)
00745     ths->r2r_kind[t] = FFTW_REDFT00;
00746 
00747   NFCT_WINDOW_HELP_INIT;
00748 
00749   if (ths->nfct_flags & MALLOC_X)
00750     ths->x = (R*)Y(malloc)(ths->d * ths->M_total * sizeof(R));
00751 
00752   if (ths->nfct_flags & MALLOC_F_HAT)
00753     ths->f_hat = (R*)Y(malloc)(ths->N_total * sizeof(R));
00754 
00755   if (ths->nfct_flags & MALLOC_F)
00756     ths->f = (R*)Y(malloc)(ths->M_total * sizeof(R));
00757 
00758   if (ths->nfct_flags & PRE_PHI_HUT)
00759     precompute_phi_hut(ths);
00760 
00761   /* NO FFTW_MALLOC HERE */
00762   if (ths->nfct_flags & PRE_PSI)
00763   {
00764     ths->psi =
00765       (R*)Y(malloc)(ths->M_total * ths->d * NFCT_SUMMANDS * sizeof(R));
00766 
00767     /* Set default for full_psi_eps */
00768     ths->nfct_full_psi_eps = POW(K(10.0), K(-10.0));
00769   }
00770 
00771   if (ths->nfct_flags & FFTW_INIT)
00772   {
00773     ths->g1 =
00774       (R*)Y(malloc)(prod_int(ths->n, ths->d) * sizeof(R));
00775 
00776     if (ths->nfct_flags & FFT_OUT_OF_PLACE)
00777       ths->g2 =
00778         (R*) Y(malloc)(prod_int(ths->n, ths->d) * sizeof(R));
00779     else
00780       ths->g2 = ths->g1;
00781 
00782     ths->my_fftw_r2r_plan =
00783       Z(plan_r2r)(ths->d, ths->n, ths->g1, ths->g2, ths->r2r_kind,
00784         ths->fftw_flags);
00785   }
00786 
00787   ths->mv_trafo = (void (*) (void* ))X(trafo);
00788   ths->mv_adjoint = (void (*) (void* ))X(adjoint);
00789 }
00790 
00791 void X(init)(X(plan) *ths, int d, int *N, int M_total)
00792 {
00793   int t;
00794 
00795   ths->d = d;
00796   ths->M_total = M_total;
00797   ths->N = (int*) Y(malloc)(ths->d * sizeof(int));
00798 
00799   for (t = 0;t < d; t++)
00800     ths->N[t] = N[t];
00801 
00802   ths->n = (int*) Y(malloc)(ths->d * sizeof(int));
00803 
00804   for (t = 0; t < d; t++)
00805     ths->n[t] = fftw_2N(Y(next_power_of_2)(ths->N[t]));
00806 
00807 /* Was soll dieser Ausdruck machen? Es handelt sich um eine Ganzzahl!
00808 
00809   WINDOW_HELP_ESTIMATE_m;
00810 */
00811   ths->nfct_flags = NFCT_DEFAULT_FLAGS;
00812   ths->fftw_flags = FFTW_DEFAULT_FLAGS;
00813 
00814   init_help(ths);
00815 }
00816 
00817 /* Was macht diese Funktion. Wird sie gebraucht? Bei NFST ist sie auch in
00818  * nfft3.h deklariert.
00819 void nfct_init_m(nfct_plan *ths, int d, int *N, int M_total, int m)
00820 {
00821   int t, n[d];
00822 
00823   for(t = 0; t < d; t++)
00824     n[t] = fftw_2N(X(next_power_of_2)(N[t]));
00825 
00826   nfct_init_guru(ths, d, N, M_total, n, m, NFCT_DEFAULT_FLAGS, FFTW_DEFAULT_FLAGS);
00827 }
00828 */
00829 
00830 void X(init_guru)(X(plan) *ths, int d, int *N, int M_total, int *n, int m,
00831   unsigned nfct_flags, unsigned fftw_flags)
00832 {
00833   int t; /* index over all dimensions */
00834 
00835   ths->d = d;
00836   ths->M_total = M_total;
00837 
00838   ths->N = (int*)Y(malloc)(ths->d * sizeof(int));
00839 
00840   for (t = 0; t < d; t++)
00841     ths->N[t] = N[t];
00842 
00843   ths->n = (int*)Y(malloc)(ths->d * sizeof(int));
00844 
00845   for (t = 0; t < d; t++)
00846     ths->n[t] = n[t];
00847 
00848   ths->m = m;
00849 
00850   ths->nfct_flags = nfct_flags;
00851   ths->fftw_flags = fftw_flags;
00852 
00853   init_help(ths);
00854 }
00855 
00856 void X(init_1d)(X(plan) *ths, int N0, int M_total)
00857 {
00858   int N[1];
00859 
00860   N[0] = N0;
00861   X(init)(ths, 1, N, M_total);
00862 }
00863 
00864 void X(init_2d)(X(plan) *ths, int N0, int N1, int M_total)
00865 {
00866   int N[2];
00867 
00868   N[0] = N0;
00869   N[1] = N1;
00870   X(init)(ths, 2, N, M_total);
00871 }
00872 
00873 void X(init_3d)(X(plan) *ths, int N0, int N1, int N2, int M_total)
00874 {
00875   int N[3];
00876 
00877   N[0] = N0;
00878   N[1] = N1;
00879   N[2] = N2;
00880   X(init)(ths, 3, N, M_total);
00881 }
00882 
00883 void X(finalize)(X(plan) *ths)
00884 {
00885   int t; /* dimension index*/
00886 
00887   if (ths->nfct_flags & FFTW_INIT)
00888   {
00889     Z(destroy_plan)(ths->my_fftw_r2r_plan);
00890 
00891     if (ths->nfct_flags & FFT_OUT_OF_PLACE)
00892       Y(free)(ths->g2);
00893 
00894     Y(free)(ths->g1);
00895   }
00896 
00897   /* NO FFTW_FREE HERE */
00898   if (ths->nfct_flags & PRE_PSI)
00899   {
00900     Y(free)(ths->psi);
00901   }
00902 
00903   if (ths->nfct_flags & PRE_PHI_HUT)
00904   {
00905     for (t = 0; t < ths->d; t++)
00906       Y(free)(ths->c_phi_inv[t]);
00907     Y(free)(ths->c_phi_inv);
00908   }
00909 
00910   if (ths->nfct_flags & MALLOC_F)
00911     Y(free)(ths->f);
00912 
00913   if(ths->nfct_flags & MALLOC_F_HAT)
00914     Y(free)(ths->f_hat);
00915 
00916   if (ths->nfct_flags & MALLOC_X)
00917     Y(free)(ths->x);
00918 
00919   WINDOW_HELP_FINALIZE;
00920 
00921   Y(free)(ths->N);
00922   Y(free)(ths->n);
00923   Y(free)(ths->sigma);
00924 
00925   Y(free)(ths->r2r_kind);
00926 } /* nfct_finalize */

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409