NFFT  3.3.0
nfsft/simple_test_threads.c
1 /*
2  * Copyright (c) 2002, 2015 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: simple_test.c 3498 2010-05-07 18:46:08Z keiner $ */
20 
21 /* standard headers */
22 #include <stdio.h>
23 #include <math.h>
24 #include <string.h>
25 #include <stdlib.h>
26 /* It is important to include complex.h before nfft3.h. */
27 #include <complex.h>
28 #include <omp.h>
29 
30 #include "nfft3.h" /* NFFT3 header */
31 
32 #define __FES__ "E"
33 #define K(x) ((double) x)
34 
35 static void simple_test_nfsft(void)
36 {
37  const int N = 4; /* bandwidth/maximum degree */
38  const int M = 8; /* number of nodes */
39  nfsft_plan plan; /* transform plan */
40  int j, k, n; /* loop variables */
41 
42  /* precomputation (for fast polynomial transform) */
43  nfsft_precompute(N,1000.0,0U,0U);
44 
45  /* Initialize transform plan using the guru interface. All input and output
46  * arrays are allocated by nfsft_init_guru(). Computations are performed with
47  * respect to L^2-normalized spherical harmonics Y_k^n. The array of spherical
48  * Fourier coefficients is preserved during transformations. The NFFT uses a
49  * cut-off parameter m = 6. See the NFFT 3 manual for details.
50  */
51  nfsft_init_guru(&plan, N, M, NFSFT_MALLOC_X | NFSFT_MALLOC_F |
52  NFSFT_MALLOC_F_HAT | NFSFT_NORMALIZED | NFSFT_PRESERVE_F_HAT,
53  PRE_PHI_HUT | PRE_PSI | FFTW_INIT | FFT_OUT_OF_PLACE, 6);
54 
55  /* pseudo-random nodes */
56  for (j = 0; j < plan.M_total; j++)
57  {
58  plan.x[2*j]= nfft_drand48() - K(0.5);
59  plan.x[2*j+1]= K(0.5) * nfft_drand48();
60  }
61 
62  /* precomputation (for NFFT, node-dependent) */
63  nfsft_precompute_x(&plan);
64 
65  /* pseudo-random Fourier coefficients */
66  for (k = 0; k <= plan.N; k++)
67  for (n = -k; n <= k; n++)
68  plan.f_hat[NFSFT_INDEX(k,n,&plan)] =
69  nfft_drand48() - K(0.5) + _Complex_I*(nfft_drand48() - K(0.5));
70 
71  /* Direct transformation, display result. */
72  nfsft_trafo_direct(&plan);
73  printf("Vector f (NDSFT):\n");
74  for (j = 0; j < plan.M_total; j++)
75  printf("f[%+2d] = %+5.3" __FES__ " %+5.3" __FES__ "*I\n",j,
76  creal(plan.f[j]), cimag(plan.f[j]));
77 
78  printf("\n");
79 
80  /* Fast approximate transformation, display result. */
81  printf("Vector f (NDSFT):\n");
82  for (j = 0; j < plan.M_total; j++)
83  printf("f[%+2d] = %+5.3" __FES__ " %+5.3" __FES__ "*I\n",j,
84  creal(plan.f[j]), cimag(plan.f[j]));
85 
86  printf("\n");
87 
88  /* Direct adjoint transformation, display result. */
89  nfsft_adjoint_direct(&plan);
90  printf("Vector f_hat (NDSFT):\n");
91  for (k = 0; k <= plan.N; k++)
92  for (n = -k; n <= k; n++)
93  fprintf(stdout,"f_hat[%+2d,%+2d] = %+5.3" __FES__ " %+5.3" __FES__ "*I\n",k,n,
94  creal(plan.f_hat[NFSFT_INDEX(k,n,&plan)]),
95  cimag(plan.f_hat[NFSFT_INDEX(k,n,&plan)]));
96 
97  printf("\n");
98 
99  /* Fast approximate adjoint transformation, display result. */
100  nfsft_adjoint(&plan);
101  printf("Vector f_hat (NFSFT):\n");
102  for (k = 0; k <= plan.N; k++)
103  {
104  for (n = -k; n <= k; n++)
105  {
106  fprintf(stdout,"f_hat[%+2d,%+2d] = %+5.3" __FES__ " %+5.3" __FES__ "*I\n",k,n,
107  creal(plan.f_hat[NFSFT_INDEX(k,n,&plan)]),
108  cimag(plan.f_hat[NFSFT_INDEX(k,n,&plan)]));
109  }
110  }
111 
112  /* Finalize the plan. */
113  nfsft_finalize(&plan);
114 
115  /* Destroy data precomputed for fast polynomial transform. */
116  nfsft_forget();
117 }
118 
119 int main(void)
120 {
121  printf("nthreads = %d\n", nfft_get_num_threads());
122 
123  /* init */
124  fftw_init_threads();
125 
126  printf("Computing an NDSFT, an NFSFT, an adjoint NDSFT, and an adjoint NFSFT"
127  "...\n\n");
128  simple_test_nfsft();
129  return EXIT_SUCCESS;
130 }
double * x
the nodes for ,
Definition: nfft3.h:576
fftw_complex * f_hat
Fourier coefficients.
Definition: nfft3.h:576
int N
the bandwidth
Definition: nfft3.h:576
data structure for an NFSFT (nonequispaced fast spherical Fourier transform) plan with double precisi...
Definition: nfft3.h:576
NFFT_INT M_total
Total number of samples.
Definition: nfft3.h:576
fftw_complex * f
Samples.
Definition: nfft3.h:576