00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "config.h"
00022
00023
00024 #include <stdio.h>
00025 #include <math.h>
00026 #include <string.h>
00027 #include <stdlib.h>
00028 #ifdef HAVE_COMPLEX_H
00029 #include <complex.h>
00030 #endif
00031
00032
00033 #include "nfft3util.h"
00034
00035 #include "nfft3.h"
00036 #include "infft.h"
00037
00038 static void simple_test_nfsoft(int bw, int M)
00039 {
00040 nfsoft_plan plan_nfsoft;
00041 nfsoft_plan plan_ndsoft;
00043 ticks t0, t1;
00044 int j;
00045 int k, m;
00046 double d1, d2, d3;
00047 double time, error;
00048 unsigned int flags = NFSOFT_MALLOC_X | NFSOFT_MALLOC_F | NFSOFT_MALLOC_F_HAT;
00051 k = 1000;
00052 m = 5;
00061 nfsoft_init_guru(&plan_ndsoft, bw, M, flags | NFSOFT_USE_NDFT
00062 | NFSOFT_USE_DPT, PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT
00063 | MALLOC_F | FFTW_INIT | FFT_OUT_OF_PLACE, m, k);
00064
00065 nfsoft_init_guru(&plan_nfsoft, bw, M, flags, PRE_PHI_HUT | PRE_PSI | MALLOC_X
00066 | MALLOC_F_HAT | MALLOC_F | FFTW_INIT | FFT_OUT_OF_PLACE, m, k);
00067
00069 for (j = 0; j < plan_nfsoft.M_total; j++)
00070 {
00071 d1 = ((R) rand()) / RAND_MAX - 0.5;
00072 d2 = 0.5 * ((R) rand()) / RAND_MAX;
00073 d3 = ((R) rand()) / RAND_MAX - 0.5;
00074
00075 plan_nfsoft.x[3* j ] = d1;
00076 plan_nfsoft.x[3* j + 1] = d2;
00077 plan_nfsoft.x[3* j + 2] = d3;
00079 plan_ndsoft.x[3* j ] = d1;
00080 plan_ndsoft.x[3* j + 1] = d2;
00081 plan_ndsoft.x[3* j + 2] = d3;
00082 }
00083
00085 for (j = 0; j < (bw + 1) * (4* (bw +1)*(bw+1)-1)/3;j++)
00086 {
00087 d1=((R)rand())/RAND_MAX - 0.5;
00088 d2=((R)rand())/RAND_MAX - 0.5;
00089 plan_nfsoft.f_hat[j]=d1 + I*d2;
00090 plan_ndsoft.f_hat[j]=d1 + I*d2;
00091 }
00092
00093 if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20)
00094 nfft_vpr_complex(plan_nfsoft.f_hat,(bw+1)*(4*(bw+1)*(bw+1)-1)/3,"randomly generated SO(3) Fourier coefficients");
00095 else
00096 nfft_vpr_complex(plan_ndsoft.f_hat,20,"1st-20th randomly generated SO(3) Fourier coefficient");
00097
00098 printf("\n---------------------------------------------\n");
00099
00101 nfsoft_precompute(&plan_nfsoft);
00102 nfsoft_precompute(&plan_ndsoft);
00103
00104
00106 t0 = getticks();
00107 nfsoft_trafo(&plan_nfsoft);
00108 t1 = getticks();
00109 time = nfft_elapsed_seconds(t1,t0);
00110 if (plan_nfsoft.M_total<=20)
00111 nfft_vpr_complex(plan_nfsoft.f,plan_nfsoft.M_total,"NFSOFT, function samples");
00112 else
00113 nfft_vpr_complex(plan_nfsoft.f,20,"NFSOFT, 1st-20th function sample");
00114 printf(" computed in %11le seconds\n",time);
00115
00117 t0 = getticks();
00118 nfsoft_trafo(&plan_ndsoft);
00119 t1 = getticks();
00120 time = nfft_elapsed_seconds(t1,t0);
00121 if (plan_ndsoft.M_total<=20)
00122 nfft_vpr_complex(plan_ndsoft.f,plan_ndsoft.M_total,"NDSOFT, function samples");
00123 else
00124 nfft_vpr_complex(plan_ndsoft.f,20,"NDSOFT, 1st-20th function sample");
00125 printf(" computed in %11le seconds\n",time);
00126
00128 error= X(error_l_infty_complex)(plan_ndsoft.f,plan_nfsoft.f, plan_nfsoft.M_total);
00129 printf("\n The NFSOFT of bandwidth=%d for %d rotations has infty-error %11le \n",bw, M,error);
00130
00131 printf("\n---------------------------------------------\n");
00132
00133 plan_nfsoft.f[0]=1.0;
00134 plan_ndsoft.f[0]=1.0;
00135 nfft_vpr_complex(plan_ndsoft.f,plan_ndsoft.M_total, "function samples to start adjoint trafo");
00136
00138 t0 = getticks();
00139 nfsoft_adjoint(&plan_nfsoft);
00140 t1 = getticks();
00141 time = nfft_elapsed_seconds(t1,t0);
00142 if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20)
00143 nfft_vpr_complex(plan_nfsoft.f_hat,(bw+1)*(4*(bw+1)*(bw+1)-1)/3,"SO(3) Fourier coefficients");
00144 else
00145 nfft_vpr_complex(plan_nfsoft.f_hat,20,"adjoint NFSOFT, 1st-20th Fourier coefficient");
00146 printf(" computed in %11le seconds\n",time);
00147
00149 t0 = getticks();
00150 nfsoft_adjoint(&plan_ndsoft);
00151 t1 = getticks();
00152 time = nfft_elapsed_seconds(t1,t0);
00153 if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20)
00154 nfft_vpr_complex(plan_ndsoft.f_hat,(bw+1)*(4*(bw+1)*(bw+1)-1)/3,"SO(3) Fourier coefficients");
00155 else
00156 nfft_vpr_complex(plan_ndsoft.f_hat,20,"adjoint NDSOFT, 1st-20th Fourier coefficient");
00157 printf(" computed in %11le seconds\n",time);
00158
00159
00161 error=X(error_l_infty_complex)(plan_ndsoft.f_hat,plan_nfsoft.f_hat, (bw+1)*(4*(bw+1)*(bw+1)-1)/3);
00162 printf("\n The adjoint NFSOFT of bandwidth=%d for %d rotations has infty-error %11le \n",bw, M,error);
00163
00164 printf("\n---------------------------------------------\n");
00165
00167 nfsoft_finalize(&plan_ndsoft);
00168 nfsoft_finalize(&plan_nfsoft);
00169 }
00170
00187 int main(int argc, char **argv)
00188 {
00189 int N;
00190 int M;
00192 if (argc < 2)
00193 {
00194 printf(
00195 "This test programm computes the NFSOFT with maximum polynomial degree N at M input rotations\n");
00196 printf("Usage: simple_test N M \n");
00197 printf("e.g.: simple_test 8 64\n");
00198 exit(0);
00199 }
00200
00202 N = atoi(argv[1]);
00203 M = atoi(argv[2]);
00204
00205 printf(
00206 "computing an NDSOFT, an NFSOFT, an adjoint NDSOFT, and an adjoint NFSOFT\n\n");
00207
00208 simple_test_nfsoft(N, M);
00209
00210
00211
00212
00213 return EXIT_SUCCESS;
00214
00215 }