![]() |
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 #include "config.h" 00019 00020 #include <stdio.h> 00021 #include <math.h> 00022 #include <string.h> 00023 #include <stdlib.h> 00024 #ifdef HAVE_COMPLEX_H 00025 #include <complex.h> 00026 #endif 00027 00028 #include "nfft3.h" 00029 #include "infft.h" 00030 00031 int global_n; 00032 int global_d; 00033 00034 static int comp1(const void *x, const void *y) 00035 { 00036 return ((*(const R*) x) < (*(const R*) y) ? -1 : 1); 00037 } 00038 00039 static int comp2(const void *x, const void *y) 00040 { 00041 int nx0, nx1, ny0, ny1; 00042 nx0 = global_n * (int)LRINT(*((const R*) x + 0)); 00043 nx1 = global_n * (int)LRINT(*((const R*) x + 1)); 00044 ny0 = global_n * (int)LRINT(*((const R*) y + 0)); 00045 ny1 = global_n * (int)LRINT(*((const R*) y + 1)); 00046 00047 if (nx0 < ny0) 00048 return -1; 00049 else if (nx0 == ny0) 00050 if (nx1 < ny1) 00051 return -1; 00052 else 00053 return 1; 00054 else 00055 return 1; 00056 } 00057 00058 static int comp3(const void *x, const void *y) 00059 { 00060 int nx0, nx1, nx2, ny0, ny1, ny2; 00061 nx0 = global_n * (int)LRINT(*((const R*) x + 0)); 00062 nx1 = global_n * (int)LRINT(*((const R*) x + 1)); 00063 nx2 = global_n * (int)LRINT(*((const R*) x + 2)); 00064 ny0 = global_n * (int)LRINT(*((const R*) y + 0)); 00065 ny1 = global_n * (int)LRINT(*((const R*) y + 1)); 00066 ny2 = global_n * (int)LRINT(*((const R*) y + 2)); 00067 00068 if (nx0 < ny0) 00069 return -1; 00070 else if (nx0 == ny0) 00071 if (nx1 < ny1) 00072 return -1; 00073 else if (nx1 == ny1) 00074 if (nx2 < ny2) 00075 return -1; 00076 else 00077 return 1; 00078 else 00079 return 1; 00080 else 00081 return 1; 00082 } 00083 00084 static void measure_time_nfft(int d, int N, unsigned test_ndft) 00085 { 00086 int r, M, NN[d], nn[d]; 00087 R t, t_fft, t_ndft, t_nfft; 00088 ticks t0, t1; 00089 00090 NFFT(plan) p; 00091 FFTW(plan) p_fft; 00092 00093 printf("\\verb+%ld+&\t", LRINT(LOG((R)(N)) / LOG((R)(2)) * (R)(d) + K(0.5))); 00094 00095 for (r = 0, M = 1; r < d; r++) 00096 { 00097 M = N * M; 00098 NN[r] = N; 00099 nn[r] = 2 * N; 00100 } 00101 00102 NFFT(init_guru)(&p, d, NN, M, nn, 2, 00103 PRE_PHI_HUT | 00104 PRE_FULL_PSI | MALLOC_F_HAT | MALLOC_X | MALLOC_F | 00105 FFTW_INIT | FFT_OUT_OF_PLACE, 00106 FFTW_MEASURE | FFTW_DESTROY_INPUT); 00107 00108 p_fft = FFTW(plan_dft)(d, NN, p.f_hat, p.f, FFTW_FORWARD, FFTW_MEASURE); 00109 00111 NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total); 00112 00113 global_n = nn[0]; 00114 global_d = d; 00115 00116 switch (d) 00117 { 00118 case 1: 00119 { 00120 qsort(p.x, (size_t)(p.M_total), (size_t)(d) * sizeof(R), comp1); 00121 break; 00122 } 00123 case 2: 00124 { 00125 qsort(p.x, (size_t)(p.M_total), (size_t)(d) * sizeof(R), comp2); 00126 break; 00127 } 00128 case 3: 00129 { 00130 qsort(p.x, (size_t)(p.M_total), (size_t)(d) * sizeof(R), comp3); 00131 break; 00132 } 00133 } 00134 00135 NFFT(precompute_one_psi)(&p); 00136 00137 /* init pseudo random Fourier coefficients */ 00138 NFFT(vrand_unit_complex)(p.f_hat, p.N_total); 00139 00141 t_fft = K(0.0); 00142 r = 0; 00143 00144 while (t_fft < K(1.0)) 00145 { 00146 r++; 00147 t0 = getticks(); 00148 FFTW(execute)(p_fft); 00149 t1 = getticks(); 00150 t = NFFT(elapsed_seconds)(t1, t0); 00151 t_fft += t; 00152 } 00153 t_fft /= (R)(r); 00154 00155 // printf("\\verb+%.1" __FES__ "+ & \\verb+(%.1f)+ &\t",t_fft,t_fft/(p.N_total*(log(N)/log(2)*d))*auxC); 00156 printf("\\verb+%.1" __FES__ "+ &\t", t_fft); 00157 00159 if (test_ndft) 00160 { 00161 t_ndft = K(0.0); 00162 r = 0; 00163 while (t_ndft < K(1.0)) 00164 { 00165 r++; 00166 t0 = getticks(); 00167 NFFT(trafo_direct)(&p); 00168 t1 = getticks(); 00169 t = NFFT(elapsed_seconds)(t1, t0); 00170 t_ndft += t; 00171 } 00172 t_ndft /= (R)(r); 00173 //printf("\\verb+%.1" __FES__ "+ & \\verb+(%d)+&\t",t_ndft,(int)round(t_ndft/(p.N_total*p.N_total)*auxC)); 00174 printf("\\verb+%.1" __FES__ "+ &\t", t_ndft); 00175 } 00176 else 00177 // printf("\\verb+*+\t&\t&\t"); 00178 printf("\\verb+*+\t&\t"); 00179 00181 t_nfft = K(0.0); 00182 r = 0; 00183 while (t_nfft < K(1.0)) 00184 { 00185 r++; 00186 t0 = getticks(); 00187 switch (d) 00188 { 00189 case 1: 00190 { 00191 NFFT(trafo_1d)(&p); 00192 break; 00193 } 00194 case 2: 00195 { 00196 NFFT(trafo_2d)(&p); 00197 break; 00198 } 00199 case 3: 00200 { 00201 NFFT(trafo_3d)(&p); 00202 break; 00203 } 00204 default: 00205 NFFT(trafo)(&p); 00206 } 00207 t1 = getticks(); 00208 t = NFFT(elapsed_seconds)(t1, t0); 00209 t_nfft += t; 00210 } 00211 t_nfft /= (R)(r); 00212 00213 // printf("\\verb+%.1" __FES__ "+ & \\verb+(%d)+ & \\verb+(%.1" __FES__ ")+\\\\\n",t_nfft,(int)round(t_nfft/(p.N_total*(log(N)/log(2)*d))*auxC),t_nfft/t_fft); 00214 printf("\\verb+%.1" __FES__ "+ & \\verb+(%3.1" __FIS__ ")+\\\\\n", t_nfft, t_nfft / t_fft); 00215 00216 FFTW(destroy_plan)(p_fft); 00217 NFFT(finalize)(&p); 00218 } 00219 00220 static void measure_time_nfft_XXX2(int d, int N, unsigned test_ndft) 00221 { 00222 int r, M, NN[d], nn[d]; 00223 R t, t_fft, t_ndft, t_nfft; 00224 ticks t0, t1; 00225 00226 NFFT(plan) p; 00227 FFTW(plan) p_fft; 00228 00229 printf("%ld\t", LRINT(LOG((R)(N)) / LOG((R)(2)) * (R)(d) + K(0.5))); 00230 fflush(stdout); 00231 00232 for (r = 0, M = 1; r < d; r++) 00233 { 00234 M = N * M; 00235 NN[r] = N; 00236 nn[r] = 2 * N; 00237 } 00238 00239 NFFT(init_guru)(&p, d, NN, M, nn, 2, 00240 PRE_PHI_HUT | 00241 PRE_PSI | 00242 MALLOC_F_HAT | MALLOC_X | MALLOC_F | 00243 FFTW_INIT | FFT_OUT_OF_PLACE, 00244 FFTW_MEASURE | FFTW_DESTROY_INPUT); 00245 00246 p_fft = FFTW(plan_dft)(d, NN, p.f_hat, p.f, FFTW_FORWARD, FFTW_MEASURE); 00247 00248 C *swapndft = (C*) NFFT(malloc)((size_t)(p.M_total) * sizeof(C)); 00249 00251 NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total); 00252 00253 qsort(p.x, (size_t)(p.M_total), (size_t)(d) * sizeof(R), comp1); 00254 //nfft_vpr_double(p.x,p.M_total,"nodes x"); 00255 00256 NFFT(precompute_one_psi)(&p); 00257 00259 NFFT(vrand_unit_complex)(p.f_hat, p.N_total); 00260 00262 t_fft = K(0.0); 00263 r = 0; 00264 while (t_fft < K(0.1)) 00265 { 00266 r++; 00267 t0 = getticks(); 00268 FFTW(execute)(p_fft); 00269 t1 = getticks(); 00270 t = NFFT(elapsed_seconds)(t1, t0); 00271 t_fft += t; 00272 } 00273 t_fft /= (R)(r); 00274 00275 printf("%.1" __FES__ "\t", t_fft); 00276 00278 if (test_ndft) 00279 { 00280 CSWAP(p.f, swapndft); 00281 t_ndft = K(0.0); 00282 r = 0; 00283 while (t_ndft < K(0.1)) 00284 { 00285 r++; 00286 t0 = getticks(); 00287 NFFT(trafo_direct)(&p); 00288 t1 = getticks(); 00289 t = NFFT(elapsed_seconds)(t1, t0); 00290 t_ndft += t; 00291 } 00292 t_ndft /= (R)(r); 00293 printf("%.1" __FES__ "\t", t_ndft); 00294 CSWAP(p.f, swapndft); 00295 } 00296 else 00297 printf("\t"); 00298 00300 t_nfft = K(0.0); 00301 r = 0; 00302 while (t_nfft < K(0.1)) 00303 { 00304 r++; 00305 t0 = getticks(); 00306 NFFT(trafo)(&p); 00307 t1 = getticks(); 00308 t = NFFT(elapsed_seconds)(t1, t0); 00309 t_nfft += t; 00310 } 00311 t_nfft /= (R)(r); 00312 printf("%.1" __FES__ "\t", t_nfft); 00313 if (test_ndft) 00314 printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f, p.M_total)); 00315 00317 t_nfft = K(0.0); 00318 r = 0; 00319 while (t_nfft < K(0.1)) 00320 { 00321 r++; 00322 t0 = getticks(); 00323 NFFT(trafo_1d)(&p); 00324 t1 = getticks(); 00325 t = NFFT(elapsed_seconds)(t1, t0); 00326 t_nfft += t; 00327 } 00328 t_nfft /= (R)(r); 00329 printf("%.1" __FES__ "\t", t_nfft); 00330 if (test_ndft) 00331 printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f, p.M_total)); 00332 00333 printf("\n"); 00334 00335 NFFT(free)(swapndft); 00336 FFTW(destroy_plan)(p_fft); 00337 NFFT(finalize)(&p); 00338 } 00339 00340 static void measure_time_nfft_XXX3(int d, int N, unsigned test_ndft) 00341 { 00342 int r, M, NN[d], nn[d]; 00343 R t, t_fft, t_ndft, t_nfft; 00344 ticks t0, t1; 00345 00346 NFFT(plan) p; 00347 FFTW(plan) p_fft; 00348 00349 printf("%ld\t", LRINT(LOG((R)(N)) / LOG((R)(2)) * (R)(d) + K(0.5))); 00350 fflush(stdout); 00351 00352 for (r = 0, M = 1; r < d; r++) 00353 { 00354 M = N * M; 00355 NN[r] = N; 00356 nn[r] = 2 * N; 00357 } 00358 00359 NFFT(init_guru)(&p, d, NN, M, nn, 2, 00360 PRE_PHI_HUT | 00361 PRE_PSI | 00362 MALLOC_F_HAT | MALLOC_X | MALLOC_F | 00363 FFTW_INIT | FFT_OUT_OF_PLACE, 00364 FFTW_MEASURE | FFTW_DESTROY_INPUT); 00365 00366 p_fft = FFTW(plan_dft)(d, NN, p.f, p.f_hat, FFTW_BACKWARD, FFTW_MEASURE); 00367 00368 C *swapndft = (C*) NFFT(malloc)((size_t)(p.N_total) * sizeof(C)); 00369 00371 NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total); 00372 00373 qsort(p.x, (size_t)(p.M_total), (size_t)(d) * sizeof(R), comp1); 00374 //nfft_vpr_double(p.x,p.M_total,"nodes x"); 00375 00376 NFFT(precompute_one_psi)(&p); 00377 00379 NFFT(vrand_unit_complex)(p.f, p.N_total); 00380 00382 t_fft = K(0.0); 00383 r = 0; 00384 while (t_fft < K(0.1)) 00385 { 00386 r++; 00387 t0 = getticks(); 00388 FFTW(execute)(p_fft); 00389 t1 = getticks(); 00390 t = NFFT(elapsed_seconds)(t1, t0); 00391 t_fft += t; 00392 } 00393 t_fft /= (R)(r); 00394 00395 printf("%.1" __FES__ "\t", t_fft); 00396 00398 if (test_ndft) 00399 { 00400 CSWAP(p.f_hat, swapndft); 00401 t_ndft = K(0.0); 00402 r = 0; 00403 while (t_ndft < K(0.1)) 00404 { 00405 r++; 00406 t0 = getticks(); 00407 NFFT(adjoint_direct)(&p); 00408 t1 = getticks(); 00409 t = NFFT(elapsed_seconds)(t1, t0); 00410 t_ndft += t; 00411 } 00412 t_ndft /= (R)(r); 00413 printf("%.1" __FES__ "\t", t_ndft); 00414 CSWAP(p.f_hat, swapndft); 00415 } 00416 else 00417 printf("\t"); 00418 00420 t_nfft = K(0.0); 00421 r = 0; 00422 while (t_nfft < K(0.1)) 00423 { 00424 r++; 00425 t0 = getticks(); 00426 NFFT(adjoint)(&p); 00427 t1 = getticks(); 00428 t = NFFT(elapsed_seconds)(t1, t0); 00429 t_nfft += t; 00430 } 00431 t_nfft /= (R)(r); 00432 printf("%.1" __FES__ "\t", t_nfft); 00433 if (test_ndft) 00434 printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f_hat, p.N_total)); 00435 00437 t_nfft = K(0.0); 00438 r = 0; 00439 while (t_nfft < K(0.1)) 00440 { 00441 r++; 00442 t0 = getticks(); 00443 NFFT(adjoint_1d)(&p); 00444 t1 = getticks(); 00445 t = NFFT(elapsed_seconds)(t1, t0); 00446 t_nfft += t; 00447 } 00448 t_nfft /= (R)(r); 00449 printf("%.1" __FES__ "\t", t_nfft); 00450 if (test_ndft) 00451 printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f_hat, p.N_total)); 00452 00453 printf("\n"); 00454 00455 NFFT(free)(swapndft); 00456 FFTW(destroy_plan)(p_fft); 00457 NFFT(finalize)(&p); 00458 } 00459 00460 static void measure_time_nfft_XXX4(int d, int N, unsigned test_ndft) 00461 { 00462 int r, M, NN[d], nn[d]; 00463 R t, t_fft, t_ndft, t_nfft; 00464 ticks t0, t1; 00465 00466 NFFT(plan) p; 00467 FFTW(plan) p_fft; 00468 00469 printf("%ld\t", LRINT(LOG((R)(N)) / LOG((R)(2)) * (R)(d) + K(0.5))); 00470 fflush(stdout); 00471 00472 for (r = 0, M = 1; r < d; r++) 00473 { 00474 M = N * M; 00475 NN[r] = N; 00476 nn[r] = 2 * N; 00477 } 00478 00479 NFFT(init_guru)(&p, d, NN, M, nn, 4, 00480 PRE_PHI_HUT | 00481 PRE_PSI | 00482 MALLOC_F_HAT | MALLOC_X | MALLOC_F | 00483 FFTW_INIT | FFT_OUT_OF_PLACE, 00484 FFTW_MEASURE | FFTW_DESTROY_INPUT); 00485 00486 p_fft = FFTW(plan_dft)(d, NN, p.f_hat, p.f, FFTW_FORWARD, FFTW_MEASURE); 00487 00488 C *swapndft = (C*) NFFT(malloc)((size_t)(p.M_total) * sizeof(C)); 00489 00491 NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total); 00492 00493 //for(j=0;j<2*M;j+=2) 00494 // p.x[j]=0.5*p.x[j]+0.25*(p.x[j]>=0?1:-1); 00495 00496 //qsort(p.x,p.M_total,d*sizeof(double),comp1); 00497 //nfft_vpr_double(p.x,p.M_total,"nodes x"); 00498 00499 NFFT(precompute_one_psi)(&p); 00500 00502 NFFT(vrand_unit_complex)(p.f_hat, p.N_total); 00503 00505 t_fft = K(0.0); 00506 r = 0; 00507 while (t_fft < K(0.1)) 00508 { 00509 r++; 00510 t0 = getticks(); 00511 FFTW(execute)(p_fft); 00512 t1 = getticks(); 00513 t = NFFT(elapsed_seconds)(t1, t0); 00514 t_fft += t; 00515 } 00516 t_fft /= (R)(r); 00517 00518 printf("%.1" __FES__ "\t", t_fft); 00519 00521 NFFT(vrand_unit_complex)(p.f_hat, p.N_total); 00522 00524 if (test_ndft) 00525 { 00526 CSWAP(p.f, swapndft); 00527 t_ndft = K(0.0); 00528 r = 0; 00529 while (t_ndft < K(0.1)) 00530 { 00531 r++; 00532 t0 = getticks(); 00533 NFFT(trafo_direct)(&p); 00534 t1 = getticks(); 00535 t = NFFT(elapsed_seconds)(t1, t0); 00536 t_ndft += t; 00537 } 00538 t_ndft /= (R)(r); 00539 printf("%.1" __FES__ "\t", t_ndft); 00540 00541 //printf("f=%e+i%e\t",creal(p.f[0]),cimag(p.f[0])); 00542 00543 CSWAP(p.f, swapndft); 00544 } 00545 else 00546 printf("\t"); 00547 00549 t_nfft = K(0.0); 00550 r = 0; 00551 while (t_nfft < K(0.1)) 00552 { 00553 r++; 00554 t0 = getticks(); 00555 NFFT(trafo)(&p); 00556 t1 = getticks(); 00557 t = NFFT(elapsed_seconds)(t1, t0); 00558 t_nfft += t; 00559 } 00560 t_nfft /= (R)(r); 00561 printf("%.1" __FES__ "\t", t_nfft); 00562 if (test_ndft) 00563 printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f, p.M_total)); 00564 00565 //printf("f=%e+i%e\t",creal(p.f[0]),cimag(p.f[0])); 00566 00568 t_nfft = K(0.0); 00569 r = 0; 00570 while (t_nfft < K(0.1)) 00571 { 00572 r++; 00573 t0 = getticks(); 00574 NFFT(trafo_2d)(&p); 00575 t1 = getticks(); 00576 t = NFFT(elapsed_seconds)(t1, t0); 00577 t_nfft += t; 00578 } 00579 t_nfft /= (R)(r); 00580 printf("%.1" __FES__ "\t", t_nfft); 00581 if (test_ndft) 00582 printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f, p.M_total)); 00583 00584 //printf("f=%e+i%e\t",creal(p.f[0]),cimag(p.f[0])); 00585 00586 printf("\n"); 00587 00588 NFFT(free)(swapndft); 00589 FFTW(destroy_plan)(p_fft); 00590 NFFT(finalize)(&p); 00591 } 00592 00593 static void measure_time_nfft_XXX5(int d, int N, unsigned test_ndft) 00594 { 00595 int r, M, NN[d], nn[d]; 00596 R t, t_fft, t_ndft, t_nfft; 00597 ticks t0, t1; 00598 00599 NFFT(plan) p; 00600 FFTW(plan) p_fft; 00601 00602 printf("%ld\t", LRINT(LOG((R)(N)) / LOG((R)(2)) * (R)(d) + K(0.5))); 00603 fflush(stdout); 00604 00605 for (r = 0, M = 1; r < d; r++) 00606 { 00607 M = N * M; 00608 NN[r] = N; 00609 nn[r] = 2 * N; 00610 } 00611 00612 NFFT(init_guru)(&p, d, NN, M, nn, 4, 00613 PRE_PHI_HUT | 00614 PRE_PSI | 00615 MALLOC_F_HAT | MALLOC_X | MALLOC_F | 00616 FFTW_INIT | FFT_OUT_OF_PLACE, 00617 FFTW_MEASURE | FFTW_DESTROY_INPUT); 00618 00619 p_fft = FFTW(plan_dft)(d, NN, p.f, p.f_hat, FFTW_FORWARD, FFTW_MEASURE); 00620 00621 C *swapndft = (C*) NFFT(malloc)((size_t)(p.N_total) * sizeof(C)); 00622 00624 NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total); 00625 00626 //sort_nodes(p.x,p.d,p.M_total, 00627 00628 NFFT(precompute_one_psi)(&p); 00629 00631 NFFT(vrand_unit_complex)(p.f, p.M_total); 00632 00634 t_fft = K(0.0); 00635 r = 0; 00636 while (t_fft < K(0.1)) 00637 { 00638 r++; 00639 t0 = getticks(); 00640 FFTW(execute)(p_fft); 00641 t1 = getticks(); 00642 t = NFFT(elapsed_seconds)(t1, t0); 00643 t_fft += t; 00644 } 00645 t_fft /= (R)(r); 00646 00647 printf("%.1" __FES__ "\t", t_fft); 00648 00650 NFFT(vrand_unit_complex)(p.f, p.M_total); 00651 00653 if (test_ndft) 00654 { 00655 CSWAP(p.f_hat, swapndft); 00656 t_ndft = K(0.0); 00657 r = 0; 00658 while (t_ndft < K(0.1)) 00659 { 00660 r++; 00661 t0 = getticks(); 00662 NFFT(adjoint_direct)(&p); 00663 t1 = getticks(); 00664 t = NFFT(elapsed_seconds)(t1, t0); 00665 t_ndft += t; 00666 } 00667 t_ndft /= (R)(r); 00668 printf("%.1" __FES__ "\t", t_ndft); 00669 00670 //printf("\nf_hat=%e+i%e\t",creal(p.f_hat[0]),cimag(p.f_hat[0])); 00671 00672 CSWAP(p.f_hat, swapndft); 00673 } 00674 else 00675 printf("\t"); 00676 00678 t_nfft = K(0.0); 00679 r = 0; 00680 while (t_nfft < K(0.1)) 00681 { 00682 r++; 00683 t0 = getticks(); 00684 NFFT(adjoint)(&p); 00685 t1 = getticks(); 00686 t = NFFT(elapsed_seconds)(t1, t0); 00687 t_nfft += t; 00688 } 00689 t_nfft /= (R)(r); 00690 printf("%.1" __FES__ "\t", t_nfft); 00691 if (test_ndft) 00692 printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f_hat, p.N_total)); 00693 00694 //printf("\nf_hat=%e+i%e\t",creal(p.f_hat[0]),cimag(p.f_hat[0])); 00695 00697 t_nfft = K(0.0); 00698 r = 0; 00699 while (t_nfft < K(0.1)) 00700 { 00701 r++; 00702 t0 = getticks(); 00703 NFFT(adjoint_2d)(&p); 00704 t1 = getticks(); 00705 t = NFFT(elapsed_seconds)(t1, t0); 00706 t_nfft += t; 00707 } 00708 t_nfft /= (R)(r); 00709 printf("%.1" __FES__ "\t", t_nfft); 00710 if (test_ndft) 00711 printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f_hat, p.N_total)); 00712 00713 //printf("\nf_hat=%e+i%e\t",creal(p.f_hat[0]),cimag(p.f_hat[0])); 00714 00715 printf("\n"); 00716 00717 NFFT(free)(swapndft); 00718 FFTW(destroy_plan)(p_fft); 00719 NFFT(finalize)(&p); 00720 } 00721 00722 static void measure_time_nfft_XXX6(int d, int N, unsigned test_ndft) 00723 { 00724 int r, M, NN[d], nn[d]; 00725 R t, t_fft, t_ndft, t_nfft; 00726 ticks t0, t1; 00727 00728 NFFT(plan) p; 00729 FFTW(plan) p_fft; 00730 00731 printf("%ld\t", LRINT(LOG((R)(N)) / LOG((R)(2)) * (R)(d) + K(0.5))); 00732 fflush(stdout); 00733 00734 for (r = 0, M = 1; r < d; r++) 00735 { 00736 M = N * M; 00737 NN[r] = N; 00738 nn[r] = 2 * N; 00739 } 00740 00741 NFFT(init_guru)(&p, d, NN, M, nn, 2, 00742 PRE_PHI_HUT | 00743 FG_PSI | 00744 MALLOC_F_HAT | MALLOC_X | MALLOC_F | 00745 FFTW_INIT | FFT_OUT_OF_PLACE, 00746 FFTW_MEASURE | FFTW_DESTROY_INPUT); 00747 00748 p_fft = FFTW(plan_dft)(d, NN, p.f_hat, p.f, FFTW_FORWARD, FFTW_MEASURE); 00749 00750 C *swapndft = (C*) NFFT(malloc)((size_t)(p.M_total) * sizeof(C)); 00751 00753 NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total); 00754 00755 //qsort(p.x,p.M_total,d*sizeof(double),comp1); 00756 //nfft_vpr_double(p.x,p.M_total,"nodes x"); 00757 00758 NFFT(precompute_one_psi)(&p); 00759 00761 NFFT(vrand_unit_complex)(p.f_hat, p.N_total); 00762 00764 t_fft = K(0.0); 00765 r = 0; 00766 while (t_fft < K(0.1)) 00767 { 00768 r++; 00769 t0 = getticks(); 00770 FFTW(execute)(p_fft); 00771 t1 = getticks(); 00772 t = NFFT(elapsed_seconds)(t1, t0); 00773 t_fft += t; 00774 } 00775 t_fft /= (R)(r); 00776 00777 printf("%.1" __FES__ "\t", t_fft); 00778 00780 NFFT(vrand_unit_complex)(p.f_hat, p.N_total); 00781 00783 if (test_ndft) 00784 { 00785 CSWAP(p.f, swapndft); 00786 t_ndft = K(0.0); 00787 r = 0; 00788 while (t_ndft < K(0.1)) 00789 { 00790 r++; 00791 t0 = getticks(); 00792 NFFT(trafo_direct)(&p); 00793 t1 = getticks(); 00794 t = NFFT(elapsed_seconds)(t1, t0); 00795 t_ndft += t; 00796 } 00797 t_ndft /= (R)(r); 00798 printf("%.1" __FES__ "\t", t_ndft); 00799 00800 //printf("f=%e+i%e\t",creal(p.f[0]),cimag(p.f[0])); 00801 00802 CSWAP(p.f, swapndft); 00803 } 00804 else 00805 printf("\t"); 00806 00808 t_nfft = K(0.0); 00809 r = 0; 00810 while (t_nfft < K(0.1)) 00811 { 00812 r++; 00813 t0 = getticks(); 00814 NFFT(trafo)(&p); 00815 t1 = getticks(); 00816 t = NFFT(elapsed_seconds)(t1, t0); 00817 t_nfft += t; 00818 } 00819 t_nfft /= (R)(r); 00820 printf("%.1" __FES__ "\t", t_nfft); 00821 if (test_ndft) 00822 printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f, p.M_total)); 00823 00824 //printf("f=%e+i%e\t",creal(p.f[0]),cimag(p.f[0])); 00825 00827 t_nfft = K(0.0); 00828 r = 0; 00829 while (t_nfft < K(0.1)) 00830 { 00831 r++; 00832 t0 = getticks(); 00833 NFFT(trafo_3d)(&p); 00834 t1 = getticks(); 00835 t = NFFT(elapsed_seconds)(t1, t0); 00836 t_nfft += t; 00837 } 00838 t_nfft /= (R)(r); 00839 printf("%.1" __FES__ "\t", t_nfft); 00840 if (test_ndft) 00841 printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f, p.M_total)); 00842 00843 //printf("f=%e+i%e\t",creal(p.f[0]),cimag(p.f[0])); 00844 00845 printf("\n"); 00846 00847 NFFT(free)(swapndft); 00848 FFTW(destroy_plan)(p_fft); 00849 NFFT(finalize)(&p); 00850 } 00851 00852 static void measure_time_nfft_XXX7(int d, int N, unsigned test_ndft) 00853 { 00854 int r, M, NN[d], nn[d]; 00855 R t, t_fft, t_ndft, t_nfft; 00856 ticks t0, t1; 00857 00858 NFFT(plan) p; 00859 FFTW(plan) p_fft; 00860 00861 printf("%ld\t", LRINT(LOG((R)(N)) / LOG((R)(2)) * (R)(d) + K(0.5))); 00862 fflush(stdout); 00863 00864 for (r = 0, M = 1; r < d; r++) 00865 { 00866 M = N * M; 00867 NN[r] = N; 00868 nn[r] = 2 * N; 00869 } 00870 00871 NFFT(init_guru)(&p, d, NN, M, nn, 2, 00872 PRE_PHI_HUT | 00873 FG_PSI | 00874 MALLOC_F_HAT | MALLOC_X | MALLOC_F | 00875 FFTW_INIT | FFT_OUT_OF_PLACE, 00876 FFTW_MEASURE | FFTW_DESTROY_INPUT); 00877 00878 p_fft = FFTW(plan_dft)(d, NN, p.f, p.f_hat, FFTW_FORWARD, FFTW_MEASURE); 00879 00880 C *swapndft = (C*) NFFT(malloc)((size_t)(p.N_total) * sizeof(C)); 00881 00883 NFFT(vrand_shifted_unit_double)(p.x, p.d * p.M_total); 00884 00885 //sort_nodes(p.x,p.d,p.M_total, 00886 00887 NFFT(precompute_one_psi)(&p); 00888 00890 NFFT(vrand_unit_complex)(p.f, p.M_total); 00891 00893 t_fft = K(0.0); 00894 r = 0; 00895 while (t_fft < K(0.1)) 00896 { 00897 r++; 00898 t0 = getticks(); 00899 FFTW(execute)(p_fft); 00900 t1 = getticks(); 00901 t = NFFT(elapsed_seconds)(t1, t0); 00902 t_fft += t; 00903 } 00904 t_fft /= (R)(r); 00905 00906 printf("%.1" __FES__ "\t", t_fft); 00907 00909 NFFT(vrand_unit_complex)(p.f, p.M_total); 00910 00912 if (test_ndft) 00913 { 00914 CSWAP(p.f_hat, swapndft); 00915 t_ndft = K(0.0); 00916 r = 0; 00917 while (t_ndft < K(0.1)) 00918 { 00919 r++; 00920 t0 = getticks(); 00921 NFFT(adjoint_direct)(&p); 00922 t1 = getticks(); 00923 t = NFFT(elapsed_seconds)(t1, t0); 00924 t_ndft += t; 00925 } 00926 t_ndft /= (R)(r); 00927 printf("%.1" __FES__ "\t", t_ndft); 00928 00929 //printf("\nf_hat=%e+i%e\t",creal(p.f_hat[0]),cimag(p.f_hat[0])); 00930 00931 CSWAP(p.f_hat, swapndft); 00932 } 00933 else 00934 printf("\t"); 00935 00937 t_nfft = K(0.0); 00938 r = 0; 00939 while (t_nfft < K(0.1)) 00940 { 00941 r++; 00942 t0 = getticks(); 00943 NFFT(adjoint)(&p); 00944 t1 = getticks(); 00945 t = NFFT(elapsed_seconds)(t1, t0); 00946 t_nfft += t; 00947 } 00948 t_nfft /= (R)(r); 00949 printf("%.1" __FES__ "\t", t_nfft); 00950 if (test_ndft) 00951 printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f_hat, p.N_total)); 00952 00953 //printf("\nf_hat=%e+i%e\t",creal(p.f_hat[0]),cimag(p.f_hat[0])); 00954 00956 t_nfft = K(0.0); 00957 r = 0; 00958 while (t_nfft < K(0.1)) 00959 { 00960 r++; 00961 t0 = getticks(); 00962 NFFT(adjoint_3d)(&p); 00963 t1 = getticks(); 00964 t = NFFT(elapsed_seconds)(t1, t0); 00965 t_nfft += t; 00966 } 00967 t_nfft /= (R)(r); 00968 printf("%.1" __FES__ "\t", t_nfft); 00969 if (test_ndft) 00970 printf("(%.1" __FES__ ")\t", NFFT(error_l_2_complex)(swapndft, p.f_hat, p.N_total)); 00971 00972 //printf("\nf_hat=%e+i%e\t",creal(p.f_hat[0]),cimag(p.f_hat[0])); 00973 00974 printf("\n"); 00975 00976 NFFT(free)(swapndft); 00977 FFTW(destroy_plan)(p_fft); 00978 NFFT(finalize)(&p); 00979 } 00980 00981 //static int main(void) 00982 //{ 00983 // int l, d, logIN; 00984 // 00985 // for (l = 3; l <= 6; l++) 00986 // { 00987 // d = 3; 00988 // logIN = d * l; 00989 // int N = (int)(1U << (logIN / d)); 00990 // measure_time_nfft_XXX6(d, N, logIN <= 15 ? 1 : 0); 00991 // measure_time_nfft_XXX7(d, N, logIN <= 15 ? 1 : 0); 00992 // } 00993 // 00994 // for (l = 7; l <= 12; l++) 00995 // { 00996 // d = 2; 00997 // logIN = d * l; 00998 // int N = (int)(1U << (logIN / d)); 00999 // measure_time_nfft_XXX4(d, N, logIN <= 15 ? 1 : 0); 01000 // measure_time_nfft_XXX5(d, N, logIN <= 15 ? 1 : 0); 01001 // } 01002 // 01003 // for (l = 3; l <= 12; l++) 01004 // { 01005 // logIN = l; 01006 // int N = (int)(1U << (logIN)); 01007 // measure_time_nfft_XXX2(1, N, logIN <= 15 ? 1 : 0); 01008 // measure_time_nfft_XXX3(1, N, logIN <= 15 ? 1 : 0); 01009 // } 01010 // 01011 // return EXIT_SUCCESS; 01012 //} 01013 01014 int main(void) 01015 { 01016 int l, d, logIN; 01017 01018 UNUSED(measure_time_nfft_XXX2); 01019 UNUSED(measure_time_nfft_XXX3); 01020 UNUSED(measure_time_nfft_XXX4); 01021 UNUSED(measure_time_nfft_XXX5); 01022 UNUSED(measure_time_nfft_XXX6); 01023 UNUSED(measure_time_nfft_XXX7); 01024 01025 printf("\\hline $l_N$ & FFT & NDFT & NFFT & NFFT/FFT\\\\\n"); 01026 printf("\\hline \\hline \\multicolumn{5}{|c|}{$d=1$} \\\\ \\hline\n"); 01027 for (l = 3; l <= 22; l++) 01028 { 01029 d = 1; 01030 logIN = l; 01031 int N = (int)(1U << (logIN / d)); 01032 measure_time_nfft(d, N, logIN <= 15 ? 1 : 0); 01033 01034 fflush(stdout); 01035 } 01036 01037 printf("\\hline $l_N$ & FFT & NDFT & NFFT & NFFT/FFT\\\\\n"); 01038 printf("\\hline \\hline \\multicolumn{5}{|c|}{$d=2$} \\\\ \\hline\n"); 01039 for (l = 3; l <= 11; l++) 01040 { 01041 d = 2; 01042 logIN = d * l; 01043 int N = (int)(1U << (logIN / d)); 01044 measure_time_nfft(d, N, logIN <= 15 ? 1 : 0); 01045 01046 fflush(stdout); 01047 } 01048 01049 printf("\\hline \\hline \\multicolumn{5}{|c|}{$d=3$} \\\\ \\hline\n"); 01050 for (l = 3; l <= 7; l++) 01051 { 01052 d = 3; 01053 logIN = d * l; 01054 int N = (int)(1U << (logIN / d)); 01055 measure_time_nfft(d, N, logIN <= 15 ? 1 : 0); 01056 01057 fflush(stdout); 01058 } 01059 01060 return 1; 01061 }