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