![]() |
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 <stdlib.h> 00020 #include <math.h> 00021 #include <complex.h> 00022 00023 #include "nfft3.h" 00024 00026 #define CSWAP(x,y) {double _Complex * NFFT_SWAP_temp__; \ 00027 NFFT_SWAP_temp__=(x); (x)=(y); (y)=NFFT_SWAP_temp__;} 00028 00029 void simple_test_nnfft_1d(void) 00030 { 00031 int j,k; 00032 nnfft_plan my_plan; 00034 int N[1]; 00035 N[0]=10; 00036 00038 nnfft_init(&my_plan, 1, 3, 19, N); 00039 00041 for(j=0;j<my_plan.M_total;j++) 00042 { 00043 my_plan.x[j]=((double)rand())/((double)RAND_MAX)-0.5; 00044 } 00046 for(j=0;j<my_plan.N_total;j++) 00047 { 00048 my_plan.v[j]=((double)rand())/((double)RAND_MAX)-0.5; 00049 } 00050 00052 /* if(my_plan.nnfft_flags & PRE_PSI) 00053 nnfft_precompute_psi(&my_plan); 00054 00055 if(my_plan.nnfft_flags & PRE_FULL_PSI) 00056 nnfft_precompute_full_psi(&my_plan); 00057 if(my_plan.nnfft_flags & PRE_LIN_PSI) 00058 nnfft_precompute_lin_psi(&my_plan); 00060 /* if(my_plan.nnfft_flags & PRE_PHI_HUT) 00061 nnfft_precompute_phi_hut(&my_plan); 00062 */ 00063 00064 nnfft_precompute_one_psi(&my_plan); 00065 00066 00068 for(k=0;k<my_plan.N_total;k++) 00069 my_plan.f_hat[k] = ((double)rand())/((double)RAND_MAX) + _Complex_I*((double)rand())/((double)RAND_MAX); 00070 00071 nfft_vpr_complex(my_plan.f_hat,my_plan.N_total,"given Fourier coefficients, vector f_hat"); 00072 00074 nnfft_trafo_direct(&my_plan); 00075 nfft_vpr_complex(my_plan.f,my_plan.M_total,"nndft, vector f"); 00076 00078 nnfft_trafo(&my_plan); 00079 nfft_vpr_complex(my_plan.f,my_plan.M_total,"nnfft, vector f"); 00080 00082 nnfft_finalize(&my_plan); 00083 } 00084 00085 static void simple_test_adjoint_nnfft_1d(void) 00086 { 00087 int j; 00088 nnfft_plan my_plan; 00090 int N[1]; 00091 N[0]=12; 00092 00094 nnfft_init(&my_plan, 1, 20, 33, N); 00095 00097 for(j=0;j<my_plan.M_total;j++) 00098 { 00099 my_plan.x[j]=((double)rand())/((double)RAND_MAX)-0.5; 00100 } 00102 for(j=0;j<my_plan.N_total;j++) 00103 { 00104 my_plan.v[j]=((double)rand())/((double)RAND_MAX)-0.5; 00105 } 00106 00108 if(my_plan.nnfft_flags & PRE_PSI) 00109 nnfft_precompute_psi(&my_plan); 00110 00111 if(my_plan.nnfft_flags & PRE_FULL_PSI) 00112 nnfft_precompute_full_psi(&my_plan); 00113 00114 if(my_plan.nnfft_flags & PRE_LIN_PSI) 00115 nnfft_precompute_lin_psi(&my_plan); 00116 00118 if(my_plan.nnfft_flags & PRE_PHI_HUT) 00119 nnfft_precompute_phi_hut(&my_plan); 00120 00122 for(j=0;j<my_plan.M_total;j++) 00123 my_plan.f[j] = ((double)rand())/((double)RAND_MAX) + _Complex_I*((double)rand())/((double)RAND_MAX); 00124 00125 nfft_vpr_complex(my_plan.f,my_plan.M_total,"given Samples, vector f"); 00126 00128 nnfft_adjoint_direct(&my_plan); 00129 nfft_vpr_complex(my_plan.f_hat,my_plan.N_total,"adjoint nndft, vector f_hat"); 00130 00132 nnfft_adjoint(&my_plan); 00133 nfft_vpr_complex(my_plan.f_hat,my_plan.N_total,"adjoint nnfft, vector f_hat"); 00134 00136 nnfft_finalize(&my_plan); 00137 } 00138 00139 static void simple_test_nnfft_2d(void) 00140 { 00141 int j,k; 00142 nnfft_plan my_plan; 00144 int N[2]; 00145 N[0]=12; 00146 N[1]=14; 00147 00149 nnfft_init(&my_plan, 2,12*14,19, N); 00150 00152 for(j=0;j<my_plan.M_total;j++) 00153 { 00154 my_plan.x[2*j]=((double)rand())/((double)RAND_MAX)-0.5; 00155 my_plan.x[2*j+1]=((double)rand())/((double)RAND_MAX)-0.5; 00156 } 00157 00159 for(j=0;j<my_plan.N_total;j++) 00160 { 00161 my_plan.v[2*j]=((double)rand())/((double)RAND_MAX)-0.5; 00162 my_plan.v[2*j+1]=((double)rand())/((double)RAND_MAX)-0.5; 00163 } 00164 00166 if(my_plan.nnfft_flags & PRE_PSI) 00167 nnfft_precompute_psi(&my_plan); 00168 00169 if(my_plan.nnfft_flags & PRE_FULL_PSI) 00170 nnfft_precompute_full_psi(&my_plan); 00171 00172 if(my_plan.nnfft_flags & PRE_LIN_PSI) 00173 nnfft_precompute_lin_psi(&my_plan); 00174 00176 if(my_plan.nnfft_flags & PRE_PHI_HUT) 00177 nnfft_precompute_phi_hut(&my_plan); 00178 00180 for(k=0;k<my_plan.N_total;k++) 00181 my_plan.f_hat[k] = ((double)rand())/((double)RAND_MAX) + _Complex_I*((double)rand())/((double)RAND_MAX); 00182 00183 nfft_vpr_complex(my_plan.f_hat,12, 00184 "given Fourier coefficients, vector f_hat (first 12 entries)"); 00185 00187 nnfft_trafo_direct(&my_plan); 00188 nfft_vpr_complex(my_plan.f,my_plan.M_total,"ndft, vector f"); 00189 00191 nnfft_trafo(&my_plan); 00192 nfft_vpr_complex(my_plan.f,my_plan.M_total,"nfft, vector f"); 00193 00195 nnfft_finalize(&my_plan); 00196 } 00197 00198 static void simple_test_innfft_1d(void) 00199 { 00200 int j,k,l,N=8; 00201 nnfft_plan my_plan; 00202 solver_plan_complex my_iplan; 00205 nnfft_init(&my_plan,1 ,8 ,8 ,&N); 00206 00208 solver_init_advanced_complex(&my_iplan,(nfft_mv_plan_complex*)(&my_plan),CGNR); 00209 00211 for(j=0;j<my_plan.M_total;j++) 00212 my_plan.x[j]=((double)rand())/((double)RAND_MAX)-0.5; 00213 00215 for(k=0;k<my_plan.N_total;k++) 00216 my_plan.v[k]=((double)rand())/((double)RAND_MAX)-0.5; 00217 00219 if(my_plan.nnfft_flags & PRE_PSI) 00220 nnfft_precompute_psi(&my_plan); 00221 00222 if(my_plan.nnfft_flags & PRE_FULL_PSI) 00223 nnfft_precompute_full_psi(&my_plan); 00224 00226 if(my_plan.nnfft_flags & PRE_PHI_HUT) 00227 nnfft_precompute_phi_hut(&my_plan); 00228 00230 for(j=0;j<my_plan.M_total;j++) 00231 my_iplan.y[j] = ((double)rand())/((double)RAND_MAX); 00232 00233 nfft_vpr_complex(my_iplan.y,my_plan.M_total,"given data, vector given_f"); 00234 00236 for(k=0;k<my_plan.N_total;k++) 00237 my_iplan.f_hat_iter[k] = 0.0; 00238 00239 nfft_vpr_complex(my_iplan.f_hat_iter,my_plan.N_total, 00240 "approximate solution, vector f_hat_iter"); 00241 00243 solver_before_loop_complex(&my_iplan); 00244 00245 for(l=0;l<8;l++) 00246 { 00247 printf("iteration l=%d\n",l); 00248 solver_loop_one_step_complex(&my_iplan); 00249 nfft_vpr_complex(my_iplan.f_hat_iter,my_plan.N_total, 00250 "approximate solution, vector f_hat_iter"); 00251 00252 CSWAP(my_iplan.f_hat_iter,my_plan.f_hat); 00253 nnfft_trafo(&my_plan); 00254 nfft_vpr_complex(my_plan.f,my_plan.M_total,"fitting the data, vector f"); 00255 CSWAP(my_iplan.f_hat_iter,my_plan.f_hat); 00256 } 00257 00258 solver_finalize_complex(&my_iplan); 00259 nnfft_finalize(&my_plan); 00260 } 00261 00262 static void measure_time_nnfft_1d(void) 00263 { 00264 int j,k; 00265 nnfft_plan my_plan; 00266 int my_N,N=4; 00267 double t; 00268 double t0, t1; 00269 00270 for(my_N=16; my_N<=16384; my_N*=2) 00271 { 00272 nnfft_init(&my_plan,1,my_N,my_N,&N); 00273 00274 for(j=0;j<my_plan.M_total;j++) 00275 my_plan.x[j]=((double)rand())/((double)RAND_MAX)-0.5; 00276 00277 for(j=0;j<my_plan.N_total;j++) 00278 my_plan.v[j]=((double)rand())/((double)RAND_MAX)-0.5; 00279 00280 if(my_plan.nnfft_flags & PRE_PSI) 00281 nnfft_precompute_psi(&my_plan); 00282 00283 if(my_plan.nnfft_flags & PRE_FULL_PSI) 00284 nnfft_precompute_full_psi(&my_plan); 00285 00286 if(my_plan.nnfft_flags & PRE_PHI_HUT) 00287 nnfft_precompute_phi_hut(&my_plan); 00288 00289 for(k=0;k<my_plan.N_total;k++) 00290 my_plan.f_hat[k] = ((double)rand())/((double)RAND_MAX) + _Complex_I*((double)rand())/((double)RAND_MAX); 00291 00292 t0 = nfft_clock_gettime_seconds(); 00293 nnfft_trafo_direct(&my_plan); 00294 t1 = nfft_clock_gettime_seconds(); 00295 t = t1 - t0; 00296 printf("t_nndft=%e,\t",t); 00297 00298 t0 = nfft_clock_gettime_seconds(); 00299 nnfft_trafo(&my_plan); 00300 t1 = nfft_clock_gettime_seconds(); 00301 t = t1 - t0; 00302 printf("t_nnfft=%e\t",t); 00303 00304 printf("(N=M=%d)\n",my_N); 00305 00306 nnfft_finalize(&my_plan); 00307 } 00308 } 00309 00310 int main(void) 00311 { 00312 system("clear"); 00313 printf("1) computing a one dimensional nndft, nnfft\n\n"); 00314 simple_test_nnfft_1d(); 00315 00316 /*getc(stdin); 00317 00318 system("clear"); 00319 printf("1a) computing a one dimensional adjoint nndft, nnfft\n\n"); 00320 simple_test_adjoint_nnfft_1d(); 00321 00322 getc(stdin); 00323 00324 system("clear"); 00325 printf("2) computing a two dimensional nndft, nnfft\n\n"); 00326 simple_test_nnfft_2d(); 00327 00328 getc(stdin); 00329 00330 system("clear"); 00331 printf("3) computing a one dimensional innfft\n\n"); 00332 simple_test_innfft_1d(); 00333 00334 getc(stdin); 00335 00336 system("clear"); 00337 printf("4) computing times for one dimensional nnfft\n\n"); 00338 measure_time_nnfft_1d(); 00339 */ 00340 return 1; 00341 }