1%TEST_NFSFT Example program:  of fast algorithm
2%Use simple_test.m for the basic usage exampe.
3
4% Copyright (c) 2002, 2017 Jens Keiner, Stefan Kunis, Daniel Potts
5%
6% This program is free software; you can redistribute it and/or modify it under
7% the terms of the GNU General Public License as published by the Free Software
8% Foundation; either version 2 of the License, or (at your option) any later
9% version.
10%
11% This program is distributed in the hope that it will be useful, but WITHOUT
12% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
13% FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
14% details.
15%
16% You should have received a copy of the GNU General Public License along with
17% this program; if not, write to the Free Software Foundation, Inc., 51
18% Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
19%
20fprintf('Number of threads: %d\n', nfsft_get_num_threads());
21
22% maximum degree (bandwidth) of spherical harmonics expansions
23N = 32;
24direct = 1;
25
26% threshold (affects accuracy, 1000 is the default)
27kappa = 10000;
28
29% random Fourier coefficients
30fh = f_hat(rand((N+1)*(N+1),1)./(1:(N+1)*(N+1))');
31
32% precomputation
33tic
34nfsft_precompute(N,kappa,0,0*FPT_NO_FAST_ALGORITHM);
35fprintf('Time of precomputation:    %g seconds\n', toc);
36
37% number of nodes
38M = (2*N+2) * (N+1);
39
40% nodes
41ph=(-N-1:N)/(2*N+2)*2*pi;
42th=(0:N)/(2*N+2)*2*pi;
43[ph,th]=meshgrid(ph,th);
44X=[ph(:)';th(:)'];
45
46fprintf('Performing a fast and direct NFSFT with bandwidth N = %d and M = %d nodes\n', N,M);
47
48% Create plan.
49plan = nfsft_init_advanced(N,M,bitor(NFSFT_NORMALIZED,NFSFT_PRESERVE_F_HAT));
50plan2 = nfsft_init_advanced(N,M,bitor(NFSFT_NORMALIZED,NFSFT_PRESERVE_F_HAT)+NFSFT_EQUISPACED);
51
52% Set nodes.
53nfsft_set_x(plan,X);
54
55% set Fourier coefficients
56nfsft_set_f_hat(plan,double(fh));
57nfsft_set_f_hat(plan2,double(fh));
58
59% direct (slow) transform
60if(direct)
61tic
62ndsft_trafo(plan);
63fprintf('Time of direct transform:  %g seconds\n', toc);
64
65% get function values
66f_direct = nfsft_get_f(plan);
67end
68
69% fast NFSFT transform
70tic
71nfsft_trafo(plan);
72fprintf('Time of fast transform:    %g seconds\n', toc);
73% get function values
74f_fast = nfsft_get_f(plan);
75
76% fast FSFT transform
77tic
78nfsft_trafo(plan2);
79fprintf('Time of FSFT transform:    %g seconds\n', toc);
80% get function values
81f_fsft = nfsft_get_f(plan2);
82
83
84% random function values for adjoint
85f = rand(M,1);
86
87% adjoint transform
88if(direct)
89nfsft_set_f(plan,f)
90tic
91ndsft_adjoint(plan);
92fprintf('Time of direct adjoint:    %g seconds\n', toc);
93
94% get function values
95fh2_direct = f_hat(nfsft_get_f_hat(plan));
96end
97
98% Fast adjoint transform
99nfsft_set_f(plan,f)
100tic
101nfsft_adjoint(plan);
102fprintf('Time of fast adjoint:      %g seconds\n', toc);
103fh2_fast = f_hat(nfsft_get_f_hat(plan));
104
105% FSFT adjoint transform
106nfsft_set_f(plan2,f)
107tic
108nfsft_adjoint(plan2);
109fprintf('Time of FSFT adjoint:      %g seconds\n', toc);
110fh2_fsft = f_hat(nfsft_get_f_hat(plan2));
111
112% finalize plan
113nfsft_finalize(plan);
114nfsft_finalize(plan2);
115
116% display error
117fprintf('Relative error trafo fast vs FSFT:       %g\n', norm(f_fast-f_fsft)/norm(f_fast));
118if(direct)
119fprintf('Relative error trafo FSFT vs direct:     %g\n', norm(f_fsft-f_direct)/norm(f_direct));
120end
121fprintf('Relative error adjoint fast vs FSFT:     %g\n', norm(fh2_fast-fh2_fsft)/norm(fh2_fast));
122if(direct)
123fprintf('Relative error adjoint FSFT vs direct:   %g\n', norm(fh2_fsft-fh2_direct)/norm(fh2_direct));
124end
125
126% destroy all precomputations
127% nfsft_forget()
128