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