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