NFFT  3.3.0
fpt/simple_test.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$ */
20 
21 /* standard headers */
22 #include <stdlib.h>
23 #include <stdio.h>
24 #include <math.h>
25 
26 /* It is important to include complex.h before nfft3.h and fftw3.h. */
27 #include <complex.h>
28 
29 #include <fftw3.h>
30 
31 /* NFFT3 header */
32 #include "nfft3.h"
33 
34 int main(void)
35 {
36  /* This example shows the use of the fast polynomial transform to evaluate a
37  * finite expansion in Legendre polynomials,
38  *
39  * f(x) = a_0 P_0(x) + a_1 P_1(x) + ... + a_N P_N(x) (1)
40  *
41  * at the Chebyshev nodes x_j = cos(j*pi/N), j=0,1,...,N. */
42  const int N = 8;
43 
44  /* An fpt_set is a data structure that contains precomputed data for a number
45  * of different polynomial transforms. Here, we need only one transform. the
46  * second parameter (t) is the exponent of the maximum transform size desired
47  * (2^t), i.e., t = 3 means that N in (1) can be at most N = 8. */
48  fpt_set set = fpt_init(1,lrint(ceil(log2((double)N))),0U);
49 
50  /* Three-term recurrence coefficients for Legendre polynomials */
51  double *alpha = malloc((N+2)*sizeof(double)),
52  *beta = malloc((N+2)*sizeof(double)),
53  *gamma = malloc((N+2)*sizeof(double));
54 
55  /* alpha[0] and beta[0] are not referenced. */
56  alpha[0] = beta[0] = 0.0;
57  /* gamma[0] contains the value of P_0(x) (which is a constant). */
58  gamma[0] = 1.0;
59 
60  /* Actual three-term recurrence coefficients for Legendre polynomials */
61  {
62  int k;
63  for (k = 0; k <= N; k++)
64  {
65  alpha[k+1] = ((double)(2*k+1))/((double)(k+1));
66  beta[k+1] = 0.0;
67  gamma[k+1] = -((double)(k))/((double)(k+1));
68  }
69  }
70 
71  printf(
72  "Computing a fast polynomial transform (FPT) and a fast discrete cosine \n"
73  "transform (DCT) to evaluate\n\n"
74  " f_j = a_0 P_0(x_j) + a_1 P_1(x_j) + ... + a_N P_N(x_j), j=0,1,...,N,\n\n"
75  "with N=%d, x_j = cos(j*pi/N), j=0,1,...N, the Chebyshev nodes, a_k,\n"
76  "k=0,1,...,N, random Fourier coefficients in [-1,1]x[-1,1]*I, and P_k,\n"
77  "k=0,1,...,N, the Legendre polynomials.",N
78  );
79 
80  /* Random seed, makes things reproducible. */
81  nfft_srand48(314);
82 
83  /* The function fpt_repcompute actually does the precomputation for a single
84  * transform. It needs arrays alpha, beta, and gamma, containing the three-
85  * term recurrence coefficients, here of the Legendre polynomials. The format
86  * is explained above. The sixth parameter (k_start) is where the index in the
87  * linear combination (1) starts, here k_start=0. The seventh parameter
88  * (kappa) is the threshold which has an influence on the accuracy of the fast
89  * polynomial transform. Usually, kappa = 1000 is a good choice. */
90  fpt_precompute(set,0,alpha,beta,gamma,0,1000.0);
91 
92 
93  {
94  /* Arrays for Fourier coefficients and function values. */
95  double _Complex *a = malloc((N+1)*sizeof(double _Complex));
96  double _Complex *b = malloc((N+1)*sizeof(double _Complex));
97  double *f = malloc((N+1)*sizeof(double _Complex));
98 
99  /* Plan for discrete cosine transform */
100  const int NP1 = N + 1;
101  fftw_r2r_kind kind = FFTW_REDFT00;
102  fftw_plan p = fftw_plan_many_r2r(1, &NP1, 1, (double*)b, NULL, 2, 1,
103  (double*)f, NULL, 1, 1, &kind, 0U);
104 
105  /* random Fourier coefficients */
106  {
107  int k;
108  printf("\n2) Random Fourier coefficients a_k, k=0,1,...,N:\n");
109  for (k = 0; k <= N; k++)
110  {
111  a[k] = 2.0*nfft_drand48() - 1.0; /* for debugging: use k+1 */
112  printf(" a_%-2d = %+5.3lE\n",k,creal(a[k]));
113  }
114  }
115 
116  /* fast polynomial transform */
117  fpt_trafo(set,0,a,b,N,0U);
118 
119  /* Renormalize coefficients b_j, j=1,2,...,N-1 owing to how FFTW defines a
120  * DCT-I; see
121  * http://www.fftw.org/fftw3_doc/1d-Real_002deven-DFTs-_0028DCTs_0029.html
122  * for details */
123  {
124  int j;
125  for (j = 1; j < N; j++)
126  b[j] *= 0.5;
127  }
128 
129  /* discrete cosine transform */
130  fftw_execute(p);
131 
132  {
133  int j;
134  printf("\n3) Function values f_j, j=1,1,...,M:\n");
135  for (j = 0; j <= N; j++)
136  printf(" f_%-2d = %+5.3lE\n",j,f[j]);
137  }
138 
139  /* cleanup */
140  free(a);
141  free(b);
142  free(f);
143 
144  /* cleanup */
145  fftw_destroy_plan(p);
146  }
147 
148  /* cleanup */
149  fpt_finalize(set);
150  free(alpha);
151  free(beta);
152  free(gamma);
153 
154  return EXIT_SUCCESS;
155 }
double * gamma
Precomputed recursion coefficients /f$^n/f$ for /f$k = 0,/ldots, N_{{max}}; n=-k,/ldots,k/f$ of associated Legendre-functions /f$P_k^n/f$.
Holds data for a set of cascade summations.
Definition: fpt.c:96
double * alpha
Precomputed recursion coefficients /f$^n/f$ for /f$k = 0,/ldots, N_{{max}}; n=-k,/ldots,k/f$ of associated Legendre-functions /f$P_k^n/f$.
double * beta
Precomputed recursion coefficients /f$^n/f$ for /f$k = 0,/ldots, N_{{max}}; n=-k,/ldots,k/f$ of associated Legendre-functions /f$P_k^n/f$.