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