NFFT Logo 3.2.2
nfsft/simple_test_threads.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: simple_test.c 3498 2010-05-07 18:46:08Z keiner $ */
00020 
00021 #include "config.h"
00022 
00023 /* standard headers */
00024 #include <stdio.h>
00025 #include <math.h>
00026 #include <string.h>
00027 #include <stdlib.h>
00028 /* It is important to include complex.h before nfft3.h. */
00029 #ifdef HAVE_COMPLEX_H
00030 #include <complex.h>
00031 #endif
00032 #include <omp.h>
00033 
00034 #include "nfft3.h" /* NFFT3 header */
00035 #include "nfft3util.h" /* NFFT3 utilities header*/
00036 #include "infft.h" /* NFFT3 internal header */
00037 
00038 static void simple_test_nfsft(void)
00039 {
00040   const int N = 4; /* bandwidth/maximum degree */
00041   const int M = 8; /* number of nodes */
00042   nfsft_plan plan; /* transform plan */
00043   int j, k, n; /* loop variables */
00044 
00045   /* precomputation (for fast polynomial transform) */
00046   nfsft_precompute(N,1000.0,0U,0U);
00047 
00048   /* Initialize transform plan using the guru interface. All input and output
00049    * arrays are allocated by nfsft_init_guru(). Computations are performed with
00050    * respect to L^2-normalized spherical harmonics Y_k^n. The array of spherical
00051    * Fourier coefficients is preserved during transformations. The NFFT uses a
00052    * cut-off parameter m = 6. See the NFFT 3 manual for details.
00053    */
00054   nfsft_init_guru(&plan, N, M, NFSFT_MALLOC_X | NFSFT_MALLOC_F |
00055     NFSFT_MALLOC_F_HAT | NFSFT_NORMALIZED | NFSFT_PRESERVE_F_HAT,
00056     PRE_PHI_HUT | PRE_PSI | FFTW_INIT | FFT_OUT_OF_PLACE, 6);
00057 
00058   /* pseudo-random nodes */
00059   for (j = 0; j < plan.M_total; j++)
00060   {
00061     plan.x[2*j]= nfft_drand48() - K(0.5);
00062     plan.x[2*j+1]= K(0.5) * nfft_drand48();
00063   }
00064 
00065   /* precomputation (for NFFT, node-dependent) */
00066   nfsft_precompute_x(&plan);
00067 
00068   /* pseudo-random Fourier coefficients */
00069   for (k = 0; k <= plan.N; k++)
00070     for (n = -k; n <= k; n++)
00071       plan.f_hat[NFSFT_INDEX(k,n,&plan)] =
00072           nfft_drand48() - K(0.5) + _Complex_I*(nfft_drand48() - K(0.5));
00073 
00074   /* Direct transformation, display result. */
00075   nfsft_trafo_direct(&plan);
00076   printf("Vector f (NDSFT):\n");
00077   for (j = 0; j < plan.M_total; j++)
00078     printf("f[%+2d] = %+5.3" FE " %+5.3" FE "*I\n",j,
00079       creal(plan.f[j]), cimag(plan.f[j]));
00080 
00081   printf("\n");
00082 
00083   /* Fast approximate transformation, display result. */
00084   printf("Vector f (NDSFT):\n");
00085   for (j = 0; j < plan.M_total; j++)
00086     printf("f[%+2d] = %+5.3" FE " %+5.3" FE "*I\n",j,
00087       creal(plan.f[j]), cimag(plan.f[j]));
00088 
00089   printf("\n");
00090 
00091   /* Direct adjoint transformation, display result. */
00092   nfsft_adjoint_direct(&plan);
00093   printf("Vector f_hat (NDSFT):\n");
00094   for (k = 0; k <= plan.N; k++)
00095     for (n = -k; n <= k; n++)
00096       fprintf(stdout,"f_hat[%+2d,%+2d] = %+5.3" FE " %+5.3" FE "*I\n",k,n,
00097         creal(plan.f_hat[NFSFT_INDEX(k,n,&plan)]),
00098         cimag(plan.f_hat[NFSFT_INDEX(k,n,&plan)]));
00099 
00100   printf("\n");
00101 
00102   /* Fast approximate adjoint transformation, display result. */
00103   nfsft_adjoint(&plan);
00104   printf("Vector f_hat (NFSFT):\n");
00105   for (k = 0; k <= plan.N; k++)
00106   {
00107     for (n = -k; n <= k; n++)
00108     {
00109       fprintf(stdout,"f_hat[%+2d,%+2d] = %+5.3" FE " %+5.3" FE "*I\n",k,n,
00110         creal(plan.f_hat[NFSFT_INDEX(k,n,&plan)]),
00111         cimag(plan.f_hat[NFSFT_INDEX(k,n,&plan)]));
00112     }
00113   }
00114 
00115   /* Finalize the plan. */
00116   nfsft_finalize(&plan);
00117 
00118   /* Destroy data precomputed for fast polynomial transform. */
00119   nfsft_forget();
00120 }
00121 
00122 int main(void)
00123 {
00124   printf("nthreads = %d\n", nfft_get_omp_num_threads());
00125 
00126   /* init */
00127   fftw_init_threads();
00128 
00129   printf("Computing an NDSFT, an NFSFT, an adjoint NDSFT, and an adjoint NFSFT"
00130     "...\n\n");
00131   simple_test_nfsft();
00132   return EXIT_SUCCESS;
00133 }

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409