nfsoft/simple_test.c
/*
* Copyright (c) 2002, 2017 Jens Keiner, Stefan Kunis, Daniel Potts
*
* This program is free software; you can redistribute it and/or modify it under
* the terms of the GNU General Public License as published by the Free Software
* Foundation; either version 2 of the License, or (at your option) any later
* version.
*
* This program is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
* FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
* details.
*
* You should have received a copy of the GNU General Public License along with
* this program; if not, write to the Free Software Foundation, Inc., 51
* Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
/* Include standard C headers. */
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <complex.h>
/* Include NFFT3 library header. */
#include "nfft3.h"
static void simple_test_nfsoft(int bw, int M)
{
nfsoft_plan plan_nfsoft;
nfsoft_plan plan_ndsoft;
double t0, t1;
int j;
int k, m;
double d1, d2, d3;
double time, error;
unsigned int flags = NFSOFT_MALLOC_X | NFSOFT_MALLOC_F | NFSOFT_MALLOC_F_HAT;
k = 1000;
m = 5;
nfsoft_init_guru(&plan_ndsoft, bw, M, flags | NFSOFT_USE_NDFT
| NFSOFT_USE_DPT, PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT
| MALLOC_F | FFTW_INIT | FFT_OUT_OF_PLACE, m, k);
nfsoft_init_guru(&plan_nfsoft, bw, M, flags, PRE_PHI_HUT | PRE_PSI | MALLOC_X
| MALLOC_F_HAT | MALLOC_F | FFTW_INIT | FFT_OUT_OF_PLACE, m, k);
for (j = 0; j < plan_nfsoft.M_total; j++)
{
d1 = ((double) rand()) / RAND_MAX - 0.5;
d2 = 0.5 * ((double) rand()) / RAND_MAX;
d3 = ((double) rand()) / RAND_MAX - 0.5;
plan_nfsoft.x[3* j ] = d1;
plan_nfsoft.x[3* j + 1] = d2;
plan_nfsoft.x[3* j + 2] = d3;
plan_ndsoft.x[3* j ] = d1;
plan_ndsoft.x[3* j + 1] = d2;
plan_ndsoft.x[3* j + 2] = d3;
}
for (j = 0; j < (bw + 1) * (4* (bw +1)*(bw+1)-1)/3;j++)
{
d1=((double)rand())/RAND_MAX - 0.5;
d2=((double)rand())/RAND_MAX - 0.5;
plan_nfsoft.f_hat[j]=d1 + I*d2;
plan_ndsoft.f_hat[j]=d1 + I*d2;
}
if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20)
nfft_vpr_complex(plan_nfsoft.f_hat,(bw+1)*(4*(bw+1)*(bw+1)-1)/3,"randomly generated SO(3) Fourier coefficients");
else
nfft_vpr_complex(plan_ndsoft.f_hat,20,"1st-20th randomly generated SO(3) Fourier coefficient");
printf("\n---------------------------------------------\n");
nfsoft_precompute(&plan_nfsoft);
nfsoft_precompute(&plan_ndsoft);
t0 = nfft_clock_gettime_seconds();
nfsoft_trafo(&plan_nfsoft);
t1 = nfft_clock_gettime_seconds();
time = t1-t0;
if (plan_nfsoft.M_total<=20)
nfft_vpr_complex(plan_nfsoft.f,plan_nfsoft.M_total,"NFSOFT, function samples");
else
nfft_vpr_complex(plan_nfsoft.f,20,"NFSOFT, 1st-20th function sample");
printf(" computed in %11le seconds\n",time);
t0 = nfft_clock_gettime_seconds();
nfsoft_trafo(&plan_ndsoft);
t1 = nfft_clock_gettime_seconds();
time = t1-t0;
if (plan_ndsoft.M_total<=20)
nfft_vpr_complex(plan_ndsoft.f,plan_ndsoft.M_total,"NDSOFT, function samples");
else
nfft_vpr_complex(plan_ndsoft.f,20,"NDSOFT, 1st-20th function sample");
printf(" computed in %11le seconds\n",time);
error= nfft_error_l_infty_complex(plan_ndsoft.f,plan_nfsoft.f, plan_nfsoft.M_total);
printf("\n The NFSOFT of bandwidth=%d for %d rotations has infty-error %11le \n",bw, M,error);
printf("\n---------------------------------------------\n");
plan_nfsoft.f[0]=1.0;
plan_ndsoft.f[0]=1.0;
nfft_vpr_complex(plan_ndsoft.f,plan_ndsoft.M_total, "function samples to start adjoint trafo");
t0 = nfft_clock_gettime_seconds();
nfsoft_adjoint(&plan_nfsoft);
t1 = nfft_clock_gettime_seconds();
time = t1-t0;
if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20)
nfft_vpr_complex(plan_nfsoft.f_hat,(bw+1)*(4*(bw+1)*(bw+1)-1)/3,"SO(3) Fourier coefficients");
else
nfft_vpr_complex(plan_nfsoft.f_hat,20,"adjoint NFSOFT, 1st-20th Fourier coefficient");
printf(" computed in %11le seconds\n",time);
t0 = nfft_clock_gettime_seconds();
nfsoft_adjoint(&plan_ndsoft);
t1 = nfft_clock_gettime_seconds();
time = t1-t0;
if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20)
nfft_vpr_complex(plan_ndsoft.f_hat,(bw+1)*(4*(bw+1)*(bw+1)-1)/3,"SO(3) Fourier coefficients");
else
nfft_vpr_complex(plan_ndsoft.f_hat,20,"adjoint NDSOFT, 1st-20th Fourier coefficient");
printf(" computed in %11le seconds\n",time);
error=nfft_error_l_infty_complex(plan_ndsoft.f_hat,plan_nfsoft.f_hat, (bw+1)*(4*(bw+1)*(bw+1)-1)/3);
printf("\n The adjoint NFSOFT of bandwidth=%d for %d rotations has infty-error %11le \n",bw, M,error);
printf("\n---------------------------------------------\n");
nfsoft_finalize(&plan_ndsoft);
nfsoft_finalize(&plan_nfsoft);
}
int main(int argc, char **argv)
{
int N;
int M;
if (argc < 2)
{
printf(
"This test programm computes the NFSOFT with maximum polynomial degree N at M input rotations\n");
printf("Usage: simple_test N M \n");
printf("e.g.: simple_test 8 64\n");
exit(0);
}
N = atoi(argv[1]);
M = atoi(argv[2]);
printf(
"computing an NDSOFT, an NFSOFT, an adjoint NDSOFT, and an adjoint NFSOFT\n\n");
simple_test_nfsoft(N, M);
/* Exit the program. */
return EXIT_SUCCESS;
}
MALLOC_F_HAT
#define MALLOC_F_HAT
Definition nfft3.h:188
MALLOC_X
#define MALLOC_X
Definition nfft3.h:187
FFT_OUT_OF_PLACE
#define FFT_OUT_OF_PLACE
Definition nfft3.h:190
PRE_PSI
#define PRE_PSI
Definition nfft3.h:185
MALLOC_F
#define MALLOC_F
Definition nfft3.h:189
FFTW_INIT
#define FFTW_INIT
Definition nfft3.h:191
PRE_PHI_HUT
#define PRE_PHI_HUT
Definition nfft3.h:181
nfsoft_adjoint
void nfsoft_adjoint(nfsoft_plan *plan_nfsoft)
Definition nfsoft.c:545
nfsoft_precompute
void nfsoft_precompute(nfsoft_plan *plan)
Definition nfsoft.c:382
NFSOFT_USE_NDFT
#define NFSOFT_USE_NDFT
Definition nfft3.h:698
nfsoft_init_guru
void nfsoft_init_guru(nfsoft_plan *plan, int N, int M, unsigned int nfsoft_flags, unsigned int nfft_flags, int nfft_cutoff, int fpt_kappa)
Definition nfsoft.c:59
nfsoft_finalize
void nfsoft_finalize(nfsoft_plan *plan)
Definition nfsoft.c:645
NFSOFT_USE_DPT
#define NFSOFT_USE_DPT
Definition nfft3.h:699
NFSOFT_MALLOC_F_HAT
#define NFSOFT_MALLOC_F_HAT
Definition nfft3.h:702
NFSOFT_MALLOC_X
#define NFSOFT_MALLOC_X
Definition nfft3.h:700
NFSOFT_MALLOC_F
#define NFSOFT_MALLOC_F
Definition nfft3.h:703
nfsoft_trafo
void nfsoft_trafo(nfsoft_plan *plan_nfsoft)
Definition nfsoft.c:416
nfft3.h
Header file for the nfft3 library.
nfft_vpr_complex
void nfft_vpr_complex(fftw_complex *x, const NFFT_INT n, const char *text)
Print complex vector to standard output.
nfsoft_plan_::f
fftw_complex * f
Samples.
Definition nfft3.h:693
nfsoft_plan_::f_hat
fftw_complex * f_hat
Fourier coefficients.
Definition nfft3.h:693
nfsoft_plan_::M_total
NFFT_INT M_total
Total number of samples.
Definition nfft3.h:693
nfsoft_plan_::x
double * x
input nodes
Definition nfft3.h:693