![]() |
NFFT
3.3.1
|
00001 /* 00002 * Copyright (c) 2002, 2016 Jens Keiner, Stefan Kunis, Daniel Potts 00003 * 00004 * This program is free software; you can redistribute it and/or modify it under 00005 * the terms of the GNU General Public License as published by the Free Software 00006 * Foundation; either version 2 of the License, or (at your option) any later 00007 * version. 00008 * 00009 * This program is distributed in the hope that it will be useful, but WITHOUT 00010 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 00011 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 00012 * details. 00013 * 00014 * You should have received a copy of the GNU General Public License along with 00015 * this program; if not, write to the Free Software Foundation, Inc., 51 00016 * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 00017 */ 00018 00019 /* Include standard C headers. */ 00020 #include <stdio.h> 00021 #include <math.h> 00022 #include <string.h> 00023 #include <stdlib.h> 00024 #include <complex.h> 00025 00026 /* Include NFFT3 library header. */ 00027 #include "nfft3.h" 00028 00029 static void simple_test_nfsoft(int bw, int M) 00030 { 00031 nfsoft_plan plan_nfsoft; 00032 nfsoft_plan plan_ndsoft; 00034 double t0, t1; 00035 int j; 00036 int k, m; 00037 double d1, d2, d3; 00038 double time, error; 00039 unsigned int flags = NFSOFT_MALLOC_X | NFSOFT_MALLOC_F | NFSOFT_MALLOC_F_HAT; 00042 k = 1000; 00043 m = 5; 00052 nfsoft_init_guru(&plan_ndsoft, bw, M, flags | NFSOFT_USE_NDFT 00053 | NFSOFT_USE_DPT, PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT 00054 | MALLOC_F | FFTW_INIT | FFT_OUT_OF_PLACE, m, k); 00055 00056 nfsoft_init_guru(&plan_nfsoft, bw, M, flags, PRE_PHI_HUT | PRE_PSI | MALLOC_X 00057 | MALLOC_F_HAT | MALLOC_F | FFTW_INIT | FFT_OUT_OF_PLACE, m, k); 00058 00060 for (j = 0; j < plan_nfsoft.M_total; j++) 00061 { 00062 d1 = ((double) rand()) / RAND_MAX - 0.5; 00063 d2 = 0.5 * ((double) rand()) / RAND_MAX; 00064 d3 = ((double) rand()) / RAND_MAX - 0.5; 00065 00066 plan_nfsoft.x[3* j ] = d1; 00067 plan_nfsoft.x[3* j + 1] = d2; 00068 plan_nfsoft.x[3* j + 2] = d3; 00070 plan_ndsoft.x[3* j ] = d1; 00071 plan_ndsoft.x[3* j + 1] = d2; 00072 plan_ndsoft.x[3* j + 2] = d3; 00073 } 00074 00076 for (j = 0; j < (bw + 1) * (4* (bw +1)*(bw+1)-1)/3;j++) 00077 { 00078 d1=((double)rand())/RAND_MAX - 0.5; 00079 d2=((double)rand())/RAND_MAX - 0.5; 00080 plan_nfsoft.f_hat[j]=d1 + I*d2; 00081 plan_ndsoft.f_hat[j]=d1 + I*d2; 00082 } 00083 00084 if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20) 00085 nfft_vpr_complex(plan_nfsoft.f_hat,(bw+1)*(4*(bw+1)*(bw+1)-1)/3,"randomly generated SO(3) Fourier coefficients"); 00086 else 00087 nfft_vpr_complex(plan_ndsoft.f_hat,20,"1st-20th randomly generated SO(3) Fourier coefficient"); 00088 00089 printf("\n---------------------------------------------\n"); 00090 00092 nfsoft_precompute(&plan_nfsoft); 00093 nfsoft_precompute(&plan_ndsoft); 00094 00095 00097 t0 = nfft_clock_gettime_seconds(); 00098 nfsoft_trafo(&plan_nfsoft); 00099 t1 = nfft_clock_gettime_seconds(); 00100 time = t1-t0; 00101 if (plan_nfsoft.M_total<=20) 00102 nfft_vpr_complex(plan_nfsoft.f,plan_nfsoft.M_total,"NFSOFT, function samples"); 00103 else 00104 nfft_vpr_complex(plan_nfsoft.f,20,"NFSOFT, 1st-20th function sample"); 00105 printf(" computed in %11le seconds\n",time); 00106 00108 t0 = nfft_clock_gettime_seconds(); 00109 nfsoft_trafo(&plan_ndsoft); 00110 t1 = nfft_clock_gettime_seconds(); 00111 time = t1-t0; 00112 if (plan_ndsoft.M_total<=20) 00113 nfft_vpr_complex(plan_ndsoft.f,plan_ndsoft.M_total,"NDSOFT, function samples"); 00114 else 00115 nfft_vpr_complex(plan_ndsoft.f,20,"NDSOFT, 1st-20th function sample"); 00116 printf(" computed in %11le seconds\n",time); 00117 00119 error= nfft_error_l_infty_complex(plan_ndsoft.f,plan_nfsoft.f, plan_nfsoft.M_total); 00120 printf("\n The NFSOFT of bandwidth=%d for %d rotations has infty-error %11le \n",bw, M,error); 00121 00122 printf("\n---------------------------------------------\n"); 00123 00124 plan_nfsoft.f[0]=1.0; 00125 plan_ndsoft.f[0]=1.0; 00126 nfft_vpr_complex(plan_ndsoft.f,plan_ndsoft.M_total, "function samples to start adjoint trafo"); 00127 00129 t0 = nfft_clock_gettime_seconds(); 00130 nfsoft_adjoint(&plan_nfsoft); 00131 t1 = nfft_clock_gettime_seconds(); 00132 time = t1-t0; 00133 if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20) 00134 nfft_vpr_complex(plan_nfsoft.f_hat,(bw+1)*(4*(bw+1)*(bw+1)-1)/3,"SO(3) Fourier coefficients"); 00135 else 00136 nfft_vpr_complex(plan_nfsoft.f_hat,20,"adjoint NFSOFT, 1st-20th Fourier coefficient"); 00137 printf(" computed in %11le seconds\n",time); 00138 00140 t0 = nfft_clock_gettime_seconds(); 00141 nfsoft_adjoint(&plan_ndsoft); 00142 t1 = nfft_clock_gettime_seconds(); 00143 time = t1-t0; 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 /* Exit the program. */ 00204 return EXIT_SUCCESS; 00205 00206 }