![]() |
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 00024 #include "config.h" 00025 00026 /* Include standard C headers. */ 00027 #include <math.h> 00028 #include <stdlib.h> 00029 #include <stdio.h> 00030 #include <string.h> 00031 #include <time.h> 00032 #ifdef HAVE_COMPLEX_H 00033 #include <complex.h> 00034 #endif 00035 00036 /* Include NFFT 3 utilities headers. */ 00037 00038 /* Include NFFT 3 library header. */ 00039 #include "nfft3.h" 00040 00041 #include "infft.h" 00042 00044 enum boolean {NO = 0, YES = 1}; 00045 00047 enum testtype {ERROR = 0, TIMING = 1}; 00048 00050 enum gridtype {GRID_GAUSS_LEGENDRE = 0, GRID_CLENSHAW_CURTIS = 1, 00051 GRID_HEALPIX = 2, GRID_EQUIDISTRIBUTION = 3, GRID_EQUIDISTRIBUTION_UNIFORM = 4}; 00052 00054 enum functiontype {FUNCTION_RANDOM_BANDLIMITED = 0, FUNCTION_F1 = 1, 00055 FUNCTION_F2 = 2, FUNCTION_F3 = 3, FUNCTION_F4 = 4, FUNCTION_F5 = 5, 00056 FUNCTION_F6 = 6}; 00057 00059 enum modes {USE_GRID = 0, RANDOM = 1}; 00060 00069 int main (int argc, char **argv) 00070 { 00071 int tc; 00072 int tc_max; 00074 int *NQ; 00076 int NQ_max; 00078 int *SQ; 00080 int SQ_max; 00081 int *RQ; 00083 int iNQ; 00084 int iNQ_max; 00085 int testfunction; 00086 int N; 00088 int use_nfsft; 00089 int use_nfft; 00090 int use_fpt; 00091 int cutoff; 00092 double threshold; 00094 int gridtype; 00095 int repetitions; 00096 int mode; 00098 double *w; 00099 double *x_grid; 00100 double *x_compare; 00101 double _Complex *f_grid; 00102 double _Complex *f_compare; 00103 double _Complex *f; 00104 double _Complex *f_hat_gen; 00106 double _Complex *f_hat; 00108 nfsft_plan plan_adjoint; 00109 nfsft_plan plan; 00110 nfsft_plan plan_gen; 00112 double t_avg; 00113 double err_infty_avg; 00114 double err_2_avg; 00116 int i; 00117 int k; 00118 int n; 00119 int d; 00121 int m_theta; 00123 int m_phi; 00125 int m_total; 00126 double *theta; 00128 double *phi; 00130 fftw_plan fplan; 00132 //int nside; /**< The size parameter for the HEALPix grid */ 00133 int d2; 00134 int M; 00135 double theta_s; 00136 double x1,x2,x3,temp; 00137 int m_compare; 00138 nfsft_plan *plan_adjoint_ptr; 00139 nfsft_plan *plan_ptr; 00140 double *w_temp; 00141 int testmode; 00142 ticks t0, t1; 00143 00144 /* Read the number of testcases. */ 00145 fscanf(stdin,"testcases=%d\n",&tc_max); 00146 fprintf(stdout,"%d\n",tc_max); 00147 00148 /* Process each testcase. */ 00149 for (tc = 0; tc < tc_max; tc++) 00150 { 00151 /* Check if the fast transform shall be used. */ 00152 fscanf(stdin,"nfsft=%d\n",&use_nfsft); 00153 fprintf(stdout,"%d\n",use_nfsft); 00154 if (use_nfsft != NO) 00155 { 00156 /* Check if the NFFT shall be used. */ 00157 fscanf(stdin,"nfft=%d\n",&use_nfft); 00158 fprintf(stdout,"%d\n",use_nfsft); 00159 if (use_nfft != NO) 00160 { 00161 /* Read the cut-off parameter. */ 00162 fscanf(stdin,"cutoff=%d\n",&cutoff); 00163 fprintf(stdout,"%d\n",cutoff); 00164 } 00165 else 00166 { 00167 /* TODO remove this */ 00168 /* Initialize unused variable with dummy value. */ 00169 cutoff = 1; 00170 } 00171 /* Check if the fast polynomial transform shall be used. */ 00172 fscanf(stdin,"fpt=%d\n",&use_fpt); 00173 fprintf(stdout,"%d\n",use_fpt); 00174 if (use_fpt != NO) 00175 { 00176 /* Read the NFSFT threshold parameter. */ 00177 fscanf(stdin,"threshold=%lf\n",&threshold); 00178 fprintf(stdout,"%lf\n",threshold); 00179 } 00180 else 00181 { 00182 /* TODO remove this */ 00183 /* Initialize unused variable with dummy value. */ 00184 threshold = 1000.0; 00185 } 00186 } 00187 else 00188 { 00189 /* TODO remove this */ 00190 /* Set dummy values. */ 00191 use_nfft = NO; 00192 use_fpt = NO; 00193 cutoff = 3; 00194 threshold = 1000.0; 00195 } 00196 00197 /* Read the testmode type. */ 00198 fscanf(stdin,"testmode=%d\n",&testmode); 00199 fprintf(stdout,"%d\n",testmode); 00200 00201 if (testmode == ERROR) 00202 { 00203 /* Read the quadrature grid type. */ 00204 fscanf(stdin,"gridtype=%d\n",&gridtype); 00205 fprintf(stdout,"%d\n",gridtype); 00206 00207 /* Read the test function. */ 00208 fscanf(stdin,"testfunction=%d\n",&testfunction); 00209 fprintf(stdout,"%d\n",testfunction); 00210 00211 /* Check if random bandlimited function has been chosen. */ 00212 if (testfunction == FUNCTION_RANDOM_BANDLIMITED) 00213 { 00214 /* Read the bandwidht. */ 00215 fscanf(stdin,"bandlimit=%d\n",&N); 00216 fprintf(stdout,"%d\n",N); 00217 } 00218 else 00219 { 00220 N = 1; 00221 } 00222 00223 /* Read the number of repetitions. */ 00224 fscanf(stdin,"repetitions=%d\n",&repetitions); 00225 fprintf(stdout,"%d\n",repetitions); 00226 00227 fscanf(stdin,"mode=%d\n",&mode); 00228 fprintf(stdout,"%d\n",mode); 00229 00230 if (mode == RANDOM) 00231 { 00232 /* Read the bandwidht. */ 00233 fscanf(stdin,"points=%d\n",&m_compare); 00234 fprintf(stdout,"%d\n",m_compare); 00235 x_compare = (double*) nfft_malloc(2*m_compare*sizeof(double)); 00236 d = 0; 00237 while (d < m_compare) 00238 { 00239 x1 = 2.0*(((double)rand())/RAND_MAX) - 1.0; 00240 x2 = 2.0*(((double)rand())/RAND_MAX) - 1.0; 00241 x3 = 2.0*(((double)rand())/RAND_MAX) - 1.0; 00242 temp = sqrt(x1*x1+x2*x2+x3*x3); 00243 if (temp <= 1) 00244 { 00245 x_compare[2*d+1] = acos(x3); 00246 if (x_compare[2*d+1] == 0 || x_compare[2*d+1] == KPI) 00247 { 00248 x_compare[2*d] = 0.0; 00249 } 00250 else 00251 { 00252 x_compare[2*d] = atan2(x2/sin(x_compare[2*d+1]),x1/sin(x_compare[2*d+1])); 00253 } 00254 x_compare[2*d] *= 1.0/(2.0*KPI); 00255 x_compare[2*d+1] *= 1.0/(2.0*KPI); 00256 d++; 00257 } 00258 } 00259 f_compare = (double _Complex*) nfft_malloc(m_compare*sizeof(double _Complex)); 00260 f = (double _Complex*) nfft_malloc(m_compare*sizeof(double _Complex)); 00261 } 00262 } 00263 00264 /* Initialize maximum cut-off degree and grid size parameter. */ 00265 NQ_max = 0; 00266 SQ_max = 0; 00267 00268 /* Read the number of cut-off degrees. */ 00269 fscanf(stdin,"bandwidths=%d\n",&iNQ_max); 00270 fprintf(stdout,"%d\n",iNQ_max); 00271 00272 /* Allocate memory for the cut-off degrees and grid size parameters. */ 00273 NQ = (int*) nfft_malloc(iNQ_max*sizeof(int)); 00274 SQ = (int*) nfft_malloc(iNQ_max*sizeof(int)); 00275 if (testmode == TIMING) 00276 { 00277 RQ = (int*) nfft_malloc(iNQ_max*sizeof(int)); 00278 } 00279 00280 /* Read the cut-off degrees and grid size parameters. */ 00281 for (iNQ = 0; iNQ < iNQ_max; iNQ++) 00282 { 00283 if (testmode == TIMING) 00284 { 00285 /* Read cut-off degree and grid size parameter. */ 00286 fscanf(stdin,"%d %d %d\n",&NQ[iNQ],&SQ[iNQ],&RQ[iNQ]); 00287 fprintf(stdout,"%d %d %d\n",NQ[iNQ],SQ[iNQ],RQ[iNQ]); 00288 NQ_max = MAX(NQ_max,NQ[iNQ]); 00289 SQ_max = MAX(SQ_max,SQ[iNQ]); 00290 } 00291 else 00292 { 00293 /* Read cut-off degree and grid size parameter. */ 00294 fscanf(stdin,"%d %d\n",&NQ[iNQ],&SQ[iNQ]); 00295 fprintf(stdout,"%d %d\n",NQ[iNQ],SQ[iNQ]); 00296 NQ_max = MAX(NQ_max,NQ[iNQ]); 00297 SQ_max = MAX(SQ_max,SQ[iNQ]); 00298 } 00299 } 00300 00301 /* Do precomputation. */ 00302 //fprintf(stderr,"NFSFT Precomputation\n"); 00303 //fflush(stderr); 00304 nfsft_precompute(NQ_max, threshold, 00305 ((use_nfsft==NO)?(NFSFT_NO_FAST_ALGORITHM):(0U)), 0U); 00306 00307 if (testmode == TIMING) 00308 { 00309 /* Allocate data structures. */ 00310 f_hat = (double _Complex*) nfft_malloc(NFSFT_F_HAT_SIZE(NQ_max)*sizeof(double _Complex)); 00311 f = (double _Complex*) nfft_malloc(SQ_max*sizeof(double _Complex)); 00312 x_grid = (double*) nfft_malloc(2*SQ_max*sizeof(double)); 00313 for (d = 0; d < SQ_max; d++) 00314 { 00315 f[d] = (((double)rand())/RAND_MAX)-0.5 + _Complex_I*((((double)rand())/RAND_MAX)-0.5); 00316 x_grid[2*d] = (((double)rand())/RAND_MAX) - 0.5; 00317 x_grid[2*d+1] = (((double)rand())/RAND_MAX) * 0.5; 00318 } 00319 } 00320 00321 //fprintf(stderr,"Entering loop\n"); 00322 //fflush(stderr); 00323 /* Process all cut-off bandwidths. */ 00324 for (iNQ = 0; iNQ < iNQ_max; iNQ++) 00325 { 00326 if (testmode == TIMING) 00327 { 00328 nfsft_init_guru(&plan,NQ[iNQ],SQ[iNQ], NFSFT_NORMALIZED | 00329 ((use_nfft!=NO)?(0U):(NFSFT_USE_NDFT)) | 00330 ((use_fpt!=NO)?(0U):(NFSFT_USE_DPT)), 00331 PRE_PHI_HUT | PRE_PSI | FFTW_INIT | FFTW_MEASURE | FFT_OUT_OF_PLACE, 00332 cutoff); 00333 00334 plan.f_hat = f_hat; 00335 plan.x = x_grid; 00336 plan.f = f; 00337 00338 nfsft_precompute_x(&plan); 00339 00340 t_avg = 0.0; 00341 00342 for (i = 0; i < RQ[iNQ]; i++) 00343 { 00344 t0 = getticks(); 00345 00346 if (use_nfsft != NO) 00347 { 00348 /* Execute the adjoint NFSFT transformation. */ 00349 nfsft_adjoint(&plan); 00350 } 00351 else 00352 { 00353 /* Execute the adjoint direct NDSFT transformation. */ 00354 nfsft_adjoint_direct(&plan); 00355 } 00356 00357 t1 = getticks(); 00358 t_avg += nfft_elapsed_seconds(t1,t0); 00359 } 00360 00361 t_avg = t_avg/((double)RQ[iNQ]); 00362 00363 nfsft_finalize(&plan); 00364 00365 fprintf(stdout,"%+le\n", t_avg); 00366 fprintf(stderr,"%d: %4d %4d %+le\n", tc, NQ[iNQ], SQ[iNQ], t_avg); 00367 } 00368 else 00369 { 00370 /* Determine the maximum number of nodes. */ 00371 switch (gridtype) 00372 { 00373 case GRID_GAUSS_LEGENDRE: 00374 /* Calculate grid dimensions. */ 00375 m_theta = SQ[iNQ] + 1; 00376 m_phi = 2*SQ[iNQ] + 2; 00377 m_total = m_theta*m_phi; 00378 break; 00379 case GRID_CLENSHAW_CURTIS: 00380 /* Calculate grid dimensions. */ 00381 m_theta = 2*SQ[iNQ] + 1; 00382 m_phi = 2*SQ[iNQ] + 2; 00383 m_total = m_theta*m_phi; 00384 break; 00385 case GRID_HEALPIX: 00386 m_theta = 1; 00387 m_phi = 12*SQ[iNQ]*SQ[iNQ]; 00388 m_total = m_theta * m_phi; 00389 //fprintf("HEALPix: SQ = %d, m_theta = %d, m_phi= %d, m"); 00390 break; 00391 case GRID_EQUIDISTRIBUTION: 00392 case GRID_EQUIDISTRIBUTION_UNIFORM: 00393 m_theta = 2; 00394 //fprintf(stderr,"ed: m_theta = %d\n",m_theta); 00395 for (k = 1; k < SQ[iNQ]; k++) 00396 { 00397 m_theta += (int)floor((2*KPI)/acos((cos(KPI/(double)SQ[iNQ])- 00398 cos(k*KPI/(double)SQ[iNQ])*cos(k*KPI/(double)SQ[iNQ]))/ 00399 (sin(k*KPI/(double)SQ[iNQ])*sin(k*KPI/(double)SQ[iNQ])))); 00400 //fprintf(stderr,"ed: m_theta = %d\n",m_theta); 00401 } 00402 //fprintf(stderr,"ed: m_theta final = %d\n",m_theta); 00403 m_phi = 1; 00404 m_total = m_theta * m_phi; 00405 break; 00406 } 00407 00408 /* Allocate memory for data structures. */ 00409 w = (double*) nfft_malloc(m_theta*sizeof(double)); 00410 x_grid = (double*) nfft_malloc(2*m_total*sizeof(double)); 00411 00412 //fprintf(stderr,"NQ = %d\n",NQ[iNQ]); 00413 //fflush(stderr); 00414 switch (gridtype) 00415 { 00416 case GRID_GAUSS_LEGENDRE: 00417 //fprintf(stderr,"Generating grid for NQ = %d, SQ = %d\n",NQ[iNQ],SQ[iNQ]); 00418 //fflush(stderr); 00419 00420 /* Read quadrature weights. */ 00421 for (k = 0; k < m_theta; k++) 00422 { 00423 fscanf(stdin,"%le\n",&w[k]); 00424 w[k] *= (2.0*KPI)/((double)m_phi); 00425 } 00426 00427 //fprintf(stderr,"Allocating theta and phi\n"); 00428 //fflush(stderr); 00429 /* Allocate memory to store the grid's angles. */ 00430 theta = (double*) nfft_malloc(m_theta*sizeof(double)); 00431 phi = (double*) nfft_malloc(m_phi*sizeof(double)); 00432 00433 //if (theta == NULL || phi == NULL) 00434 //{ 00435 //fprintf(stderr,"Couldn't allocate theta and phi\n"); 00436 //fflush(stderr); 00437 //} 00438 00439 00440 /* Read angles theta. */ 00441 for (k = 0; k < m_theta; k++) 00442 { 00443 fscanf(stdin,"%le\n",&theta[k]); 00444 } 00445 00446 /* Generate the grid angles phi. */ 00447 for (n = 0; n < m_phi; n++) 00448 { 00449 phi[n] = n/((double)m_phi); 00450 phi[n] -= ((phi[n]>=0.5)?(1.0):(0.0)); 00451 } 00452 00453 //fprintf(stderr,"Generating grid nodes\n"); 00454 //fflush(stderr); 00455 00456 /* Generate the grid's nodes. */ 00457 d = 0; 00458 for (k = 0; k < m_theta; k++) 00459 { 00460 for (n = 0; n < m_phi; n++) 00461 { 00462 x_grid[2*d] = phi[n]; 00463 x_grid[2*d+1] = theta[k]; 00464 d++; 00465 } 00466 } 00467 00468 //fprintf(stderr,"Freeing theta and phi\n"); 00469 //fflush(stderr); 00470 /* Free the arrays for the grid's angles. */ 00471 nfft_free(theta); 00472 nfft_free(phi); 00473 00474 break; 00475 00476 case GRID_CLENSHAW_CURTIS: 00477 00478 /* Allocate memory to store the grid's angles. */ 00479 theta = (double*) nfft_malloc(m_theta*sizeof(double)); 00480 phi = (double*) nfft_malloc(m_phi*sizeof(double)); 00481 00482 /* Generate the grid angles theta. */ 00483 for (k = 0; k < m_theta; k++) 00484 { 00485 theta[k] = k/((double)2*(m_theta-1)); 00486 } 00487 00488 /* Generate the grid angles phi. */ 00489 for (n = 0; n < m_phi; n++) 00490 { 00491 phi[n] = n/((double)m_phi); 00492 phi[n] -= ((phi[n]>=0.5)?(1.0):(0.0)); 00493 } 00494 00495 /* Generate quadrature weights. */ 00496 fplan = fftw_plan_r2r_1d(SQ[iNQ]+1, w, w, FFTW_REDFT00, 0U); 00497 for (k = 0; k < SQ[iNQ]+1; k++) 00498 { 00499 w[k] = -2.0/(4*k*k-1); 00500 } 00501 fftw_execute(fplan); 00502 w[0] *= 0.5; 00503 00504 for (k = 0; k < SQ[iNQ]+1; k++) 00505 { 00506 w[k] *= (2.0*KPI)/((double)(m_theta-1)*m_phi); 00507 w[m_theta-1-k] = w[k]; 00508 } 00509 fftw_destroy_plan(fplan); 00510 00511 /* Generate the grid's nodes. */ 00512 d = 0; 00513 for (k = 0; k < m_theta; k++) 00514 { 00515 for (n = 0; n < m_phi; n++) 00516 { 00517 x_grid[2*d] = phi[n]; 00518 x_grid[2*d+1] = theta[k]; 00519 d++; 00520 } 00521 } 00522 00523 /* Free the arrays for the grid's angles. */ 00524 nfft_free(theta); 00525 nfft_free(phi); 00526 00527 break; 00528 00529 case GRID_HEALPIX: 00530 00531 d = 0; 00532 for (k = 1; k <= SQ[iNQ]-1; k++) 00533 { 00534 for (n = 0; n <= 4*k-1; n++) 00535 { 00536 x_grid[2*d+1] = 1 - (k*k)/((double)(3.0*SQ[iNQ]*SQ[iNQ])); 00537 x_grid[2*d] = ((n+0.5)/(4*k)); 00538 x_grid[2*d] -= (x_grid[2*d]>=0.5)?(1.0):(0.0); 00539 d++; 00540 } 00541 } 00542 00543 d2 = d-1; 00544 00545 for (k = SQ[iNQ]; k <= 3*SQ[iNQ]; k++) 00546 { 00547 for (n = 0; n <= 4*SQ[iNQ]-1; n++) 00548 { 00549 x_grid[2*d+1] = 2.0/(3*SQ[iNQ])*(2*SQ[iNQ]-k); 00550 x_grid[2*d] = (n+((k%2==0)?(0.5):(0.0)))/(4*SQ[iNQ]); 00551 x_grid[2*d] -= (x_grid[2*d]>=0.5)?(1.0):(0.0); 00552 d++; 00553 } 00554 } 00555 00556 for (k = 1; k <= SQ[iNQ]-1; k++) 00557 { 00558 for (n = 0; n <= 4*k-1; n++) 00559 { 00560 x_grid[2*d+1] = -x_grid[2*d2+1]; 00561 x_grid[2*d] = x_grid[2*d2]; 00562 d++; 00563 d2--; 00564 } 00565 } 00566 00567 for (d = 0; d < m_total; d++) 00568 { 00569 x_grid[2*d+1] = acos(x_grid[2*d+1])/(2.0*KPI); 00570 } 00571 00572 w[0] = (4.0*KPI)/(m_total); 00573 break; 00574 00575 case GRID_EQUIDISTRIBUTION: 00576 case GRID_EQUIDISTRIBUTION_UNIFORM: 00577 /* TODO Compute the weights. */ 00578 00579 if (gridtype == GRID_EQUIDISTRIBUTION) 00580 { 00581 w_temp = (double*) nfft_malloc((SQ[iNQ]+1)*sizeof(double)); 00582 fplan = fftw_plan_r2r_1d(SQ[iNQ]/2+1, w_temp, w_temp, FFTW_REDFT00, 0U); 00583 for (k = 0; k < SQ[iNQ]/2+1; k++) 00584 { 00585 w_temp[k] = -2.0/(4*k*k-1); 00586 } 00587 fftw_execute(fplan); 00588 w_temp[0] *= 0.5; 00589 00590 for (k = 0; k < SQ[iNQ]/2+1; k++) 00591 { 00592 w_temp[k] *= (2.0*KPI)/((double)(SQ[iNQ])); 00593 w_temp[SQ[iNQ]-k] = w_temp[k]; 00594 } 00595 fftw_destroy_plan(fplan); 00596 } 00597 00598 d = 0; 00599 x_grid[2*d] = -0.5; 00600 x_grid[2*d+1] = 0.0; 00601 if (gridtype == GRID_EQUIDISTRIBUTION) 00602 { 00603 w[d] = w_temp[0]; 00604 } 00605 else 00606 { 00607 w[d] = (4.0*KPI)/(m_total); 00608 } 00609 d = 1; 00610 x_grid[2*d] = -0.5; 00611 x_grid[2*d+1] = 0.5; 00612 if (gridtype == GRID_EQUIDISTRIBUTION) 00613 { 00614 w[d] = w_temp[SQ[iNQ]]; 00615 } 00616 else 00617 { 00618 w[d] = (4.0*KPI)/(m_total); 00619 } 00620 d = 2; 00621 00622 for (k = 1; k < SQ[iNQ]; k++) 00623 { 00624 theta_s = (double)k*KPI/(double)SQ[iNQ]; 00625 M = (int)floor((2.0*KPI)/acos((cos(KPI/(double)SQ[iNQ])- 00626 cos(theta_s)*cos(theta_s))/(sin(theta_s)*sin(theta_s)))); 00627 00628 for (n = 0; n < M; n++) 00629 { 00630 x_grid[2*d] = (n + 0.5)/M; 00631 x_grid[2*d] -= (x_grid[2*d]>=0.5)?(1.0):(0.0); 00632 x_grid[2*d+1] = theta_s/(2.0*KPI); 00633 if (gridtype == GRID_EQUIDISTRIBUTION) 00634 { 00635 w[d] = w_temp[k]/((double)(M)); 00636 } 00637 else 00638 { 00639 w[d] = (4.0*KPI)/(m_total); 00640 } 00641 d++; 00642 } 00643 } 00644 00645 if (gridtype == GRID_EQUIDISTRIBUTION) 00646 { 00647 nfft_free(w_temp); 00648 } 00649 break; 00650 00651 default: 00652 break; 00653 } 00654 00655 /* Allocate memory for grid values. */ 00656 f_grid = (double _Complex*) nfft_malloc(m_total*sizeof(double _Complex)); 00657 00658 if (mode == RANDOM) 00659 { 00660 } 00661 else 00662 { 00663 m_compare = m_total; 00664 f_compare = (double _Complex*) nfft_malloc(m_compare*sizeof(double _Complex)); 00665 x_compare = x_grid; 00666 f = f_grid; 00667 } 00668 00669 //fprintf(stderr,"Generating test function\n"); 00670 //fflush(stderr); 00671 switch (testfunction) 00672 { 00673 case FUNCTION_RANDOM_BANDLIMITED: 00674 f_hat_gen = (double _Complex*) nfft_malloc(NFSFT_F_HAT_SIZE(N)*sizeof(double _Complex)); 00675 //fprintf(stderr,"Generating random test function\n"); 00676 //fflush(stderr); 00677 /* Generate random function samples by sampling a bandlimited 00678 * function. */ 00679 nfsft_init_guru(&plan_gen,N,m_total, NFSFT_NORMALIZED | 00680 ((use_nfft!=NO)?(0U):(NFSFT_USE_NDFT)) | 00681 ((use_fpt!=NO)?(0U):(NFSFT_USE_DPT)), 00682 ((N>512)?(0U):(PRE_PHI_HUT | PRE_PSI)) | FFTW_INIT | 00683 FFT_OUT_OF_PLACE, cutoff); 00684 00685 plan_gen.f_hat = f_hat_gen; 00686 plan_gen.x = x_grid; 00687 plan_gen.f = f_grid; 00688 00689 nfsft_precompute_x(&plan_gen); 00690 00691 for (k = 0; k < plan_gen.N_total; k++) 00692 { 00693 f_hat_gen[k] = 0.0; 00694 } 00695 00696 for (k = 0; k <= N; k++) 00697 { 00698 for (n = -k; n <= k; n++) 00699 { 00700 f_hat_gen[NFSFT_INDEX(k,n,&plan_gen)] = 00701 (((double)rand())/RAND_MAX)-0.5 + _Complex_I*((((double)rand())/RAND_MAX)-0.5); 00702 } 00703 } 00704 00705 if (use_nfsft != NO) 00706 { 00707 /* Execute the NFSFT transformation. */ 00708 nfsft_trafo(&plan_gen); 00709 } 00710 else 00711 { 00712 /* Execute the direct NDSFT transformation. */ 00713 nfsft_trafo_direct(&plan_gen); 00714 } 00715 00716 nfsft_finalize(&plan_gen); 00717 00718 if (mode == RANDOM) 00719 { 00720 nfsft_init_guru(&plan_gen,N,m_compare, NFSFT_NORMALIZED | 00721 ((use_nfft!=NO)?(0U):(NFSFT_USE_NDFT)) | 00722 ((use_fpt!=NO)?(0U):(NFSFT_USE_DPT)), 00723 ((N>512)?(0U):(PRE_PHI_HUT | PRE_PSI)) | FFTW_INIT | 00724 FFT_OUT_OF_PLACE, cutoff); 00725 00726 plan_gen.f_hat = f_hat_gen; 00727 plan_gen.x = x_compare; 00728 plan_gen.f = f_compare; 00729 00730 nfsft_precompute_x(&plan_gen); 00731 00732 if (use_nfsft != NO) 00733 { 00734 /* Execute the NFSFT transformation. */ 00735 nfsft_trafo(&plan_gen); 00736 } 00737 else 00738 { 00739 /* Execute the direct NDSFT transformation. */ 00740 nfsft_trafo_direct(&plan_gen); 00741 } 00742 00743 nfsft_finalize(&plan_gen); 00744 } 00745 else 00746 { 00747 memcpy(f_compare,f_grid,m_total*sizeof(double _Complex)); 00748 } 00749 00750 nfft_free(f_hat_gen); 00751 00752 break; 00753 00754 case FUNCTION_F1: 00755 for (d = 0; d < m_total; d++) 00756 { 00757 x1 = sin(x_grid[2*d+1]*2.0*KPI)*cos(x_grid[2*d]*2.0*KPI); 00758 x2 = sin(x_grid[2*d+1]*2.0*KPI)*sin(x_grid[2*d]*2.0*KPI); 00759 x3 = cos(x_grid[2*d+1]*2.0*KPI); 00760 f_grid[d] = x1*x2*x3; 00761 } 00762 if (mode == RANDOM) 00763 { 00764 for (d = 0; d < m_compare; d++) 00765 { 00766 x1 = sin(x_compare[2*d+1]*2.0*KPI)*cos(x_compare[2*d]*2.0*KPI); 00767 x2 = sin(x_compare[2*d+1]*2.0*KPI)*sin(x_compare[2*d]*2.0*KPI); 00768 x3 = cos(x_compare[2*d+1]*2.0*KPI); 00769 f_compare[d] = x1*x2*x3; 00770 } 00771 } 00772 else 00773 { 00774 memcpy(f_compare,f_grid,m_total*sizeof(double _Complex)); 00775 } 00776 break; 00777 case FUNCTION_F2: 00778 for (d = 0; d < m_total; d++) 00779 { 00780 x1 = sin(x_grid[2*d+1]*2.0*KPI)*cos(x_grid[2*d]*2.0*KPI); 00781 x2 = sin(x_grid[2*d+1]*2.0*KPI)*sin(x_grid[2*d]*2.0*KPI); 00782 x3 = cos(x_grid[2*d+1]*2.0*KPI); 00783 f_grid[d] = 0.1*exp(x1+x2+x3); 00784 } 00785 if (mode == RANDOM) 00786 { 00787 for (d = 0; d < m_compare; d++) 00788 { 00789 x1 = sin(x_compare[2*d+1]*2.0*KPI)*cos(x_compare[2*d]*2.0*KPI); 00790 x2 = sin(x_compare[2*d+1]*2.0*KPI)*sin(x_compare[2*d]*2.0*KPI); 00791 x3 = cos(x_compare[2*d+1]*2.0*KPI); 00792 f_compare[d] = 0.1*exp(x1+x2+x3); 00793 } 00794 } 00795 else 00796 { 00797 memcpy(f_compare,f_grid,m_total*sizeof(double _Complex)); 00798 } 00799 break; 00800 case FUNCTION_F3: 00801 for (d = 0; d < m_total; d++) 00802 { 00803 x1 = sin(x_grid[2*d+1]*2.0*KPI)*cos(x_grid[2*d]*2.0*KPI); 00804 x2 = sin(x_grid[2*d+1]*2.0*KPI)*sin(x_grid[2*d]*2.0*KPI); 00805 x3 = cos(x_grid[2*d+1]*2.0*KPI); 00806 temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3); 00807 f_grid[d] = 0.1*temp; 00808 } 00809 if (mode == RANDOM) 00810 { 00811 for (d = 0; d < m_compare; d++) 00812 { 00813 x1 = sin(x_compare[2*d+1]*2.0*KPI)*cos(x_compare[2*d]*2.0*KPI); 00814 x2 = sin(x_compare[2*d+1]*2.0*KPI)*sin(x_compare[2*d]*2.0*KPI); 00815 x3 = cos(x_compare[2*d+1]*2.0*KPI); 00816 temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3); 00817 f_compare[d] = 0.1*temp; 00818 } 00819 } 00820 else 00821 { 00822 memcpy(f_compare,f_grid,m_total*sizeof(double _Complex)); 00823 } 00824 break; 00825 case FUNCTION_F4: 00826 for (d = 0; d < m_total; d++) 00827 { 00828 x1 = sin(x_grid[2*d+1]*2.0*KPI)*cos(x_grid[2*d]*2.0*KPI); 00829 x2 = sin(x_grid[2*d+1]*2.0*KPI)*sin(x_grid[2*d]*2.0*KPI); 00830 x3 = cos(x_grid[2*d+1]*2.0*KPI); 00831 temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3); 00832 f_grid[d] = 1.0/(temp); 00833 } 00834 if (mode == RANDOM) 00835 { 00836 for (d = 0; d < m_compare; d++) 00837 { 00838 x1 = sin(x_compare[2*d+1]*2.0*KPI)*cos(x_compare[2*d]*2.0*KPI); 00839 x2 = sin(x_compare[2*d+1]*2.0*KPI)*sin(x_compare[2*d]*2.0*KPI); 00840 x3 = cos(x_compare[2*d+1]*2.0*KPI); 00841 temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3); 00842 f_compare[d] = 1.0/(temp); 00843 } 00844 } 00845 else 00846 { 00847 memcpy(f_compare,f_grid,m_total*sizeof(double _Complex)); 00848 } 00849 break; 00850 case FUNCTION_F5: 00851 for (d = 0; d < m_total; d++) 00852 { 00853 x1 = sin(x_grid[2*d+1]*2.0*KPI)*cos(x_grid[2*d]*2.0*KPI); 00854 x2 = sin(x_grid[2*d+1]*2.0*KPI)*sin(x_grid[2*d]*2.0*KPI); 00855 x3 = cos(x_grid[2*d+1]*2.0*KPI); 00856 temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3); 00857 f_grid[d] = 0.1*sin(1+temp)*sin(1+temp); 00858 } 00859 if (mode == RANDOM) 00860 { 00861 for (d = 0; d < m_compare; d++) 00862 { 00863 x1 = sin(x_compare[2*d+1]*2.0*KPI)*cos(x_compare[2*d]*2.0*KPI); 00864 x2 = sin(x_compare[2*d+1]*2.0*KPI)*sin(x_compare[2*d]*2.0*KPI); 00865 x3 = cos(x_compare[2*d+1]*2.0*KPI); 00866 temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3); 00867 f_compare[d] = 0.1*sin(1+temp)*sin(1+temp); 00868 } 00869 } 00870 else 00871 { 00872 memcpy(f_compare,f_grid,m_total*sizeof(double _Complex)); 00873 } 00874 break; 00875 case FUNCTION_F6: 00876 for (d = 0; d < m_total; d++) 00877 { 00878 if (x_grid[2*d+1] <= 0.25) 00879 { 00880 f_grid[d] = 1.0; 00881 } 00882 else 00883 { 00884 f_grid[d] = 1.0/(sqrt(1+3*cos(2.0*KPI*x_grid[2*d+1])*cos(2.0*KPI*x_grid[2*d+1]))); 00885 } 00886 } 00887 if (mode == RANDOM) 00888 { 00889 for (d = 0; d < m_compare; d++) 00890 { 00891 if (x_compare[2*d+1] <= 0.25) 00892 { 00893 f_compare[d] = 1.0; 00894 } 00895 else 00896 { 00897 f_compare[d] = 1.0/(sqrt(1+3*cos(2.0*KPI*x_compare[2*d+1])*cos(2.0*KPI*x_compare[2*d+1]))); 00898 } 00899 } 00900 } 00901 else 00902 { 00903 memcpy(f_compare,f_grid,m_total*sizeof(double _Complex)); 00904 } 00905 break; 00906 default: 00907 //fprintf(stderr,"Generating one function\n"); 00908 //fflush(stderr); 00909 for (d = 0; d < m_total; d++) 00910 { 00911 f_grid[d] = 1.0; 00912 } 00913 if (mode == RANDOM) 00914 { 00915 for (d = 0; d < m_compare; d++) 00916 { 00917 f_compare[d] = 1.0; 00918 } 00919 } 00920 else 00921 { 00922 memcpy(f_compare,f_grid,m_total*sizeof(double _Complex)); 00923 } 00924 break; 00925 } 00926 00927 //fprintf(stderr,"Initializing trafo\n"); 00928 //fflush(stderr); 00929 /* Init transform plan. */ 00930 nfsft_init_guru(&plan_adjoint,NQ[iNQ],m_total, NFSFT_NORMALIZED | 00931 ((use_nfft!=NO)?(0U):(NFSFT_USE_NDFT)) | 00932 ((use_fpt!=NO)?(0U):(NFSFT_USE_DPT)), 00933 ((NQ[iNQ]>512)?(0U):(PRE_PHI_HUT | PRE_PSI)) | FFTW_INIT | 00934 FFT_OUT_OF_PLACE, cutoff); 00935 00936 plan_adjoint_ptr = &plan_adjoint; 00937 00938 if (mode == RANDOM) 00939 { 00940 nfsft_init_guru(&plan,NQ[iNQ],m_compare, NFSFT_NORMALIZED | 00941 ((use_nfft!=NO)?(0U):(NFSFT_USE_NDFT)) | 00942 ((use_fpt!=NO)?(0U):(NFSFT_USE_DPT)), 00943 ((NQ[iNQ]>512)?(0U):(PRE_PHI_HUT | PRE_PSI)) | FFTW_INIT | 00944 FFT_OUT_OF_PLACE, cutoff); 00945 plan_ptr = &plan; 00946 } 00947 else 00948 { 00949 plan_ptr = &plan_adjoint; 00950 } 00951 00952 f_hat = (double _Complex*) nfft_malloc(NFSFT_F_HAT_SIZE(NQ[iNQ])*sizeof(double _Complex)); 00953 00954 plan_adjoint_ptr->f_hat = f_hat; 00955 plan_adjoint_ptr->x = x_grid; 00956 plan_adjoint_ptr->f = f_grid; 00957 00958 plan_ptr->f_hat = f_hat; 00959 plan_ptr->x = x_compare; 00960 plan_ptr->f = f; 00961 00962 //fprintf(stderr,"Precomputing for x\n"); 00963 //fflush(stderr); 00964 nfsft_precompute_x(plan_adjoint_ptr); 00965 if (plan_adjoint_ptr != plan_ptr) 00966 { 00967 nfsft_precompute_x(plan_ptr); 00968 } 00969 00970 /* Initialize cumulative time variable. */ 00971 t_avg = 0.0; 00972 err_infty_avg = 0.0; 00973 err_2_avg = 0.0; 00974 00975 /* Cycle through all runs. */ 00976 for (i = 0; i < 1/*repetitions*/; i++) 00977 { 00978 //fprintf(stderr,"Copying original values\n"); 00979 //fflush(stderr); 00980 /* Copy exact funtion values to working array. */ 00981 //memcpy(f,f_grid,m_total*sizeof(double _Complex)); 00982 00983 /* Initialize time measurement. */ 00984 t0 = getticks(); 00985 00986 //fprintf(stderr,"Multiplying with quadrature weights\n"); 00987 //fflush(stderr); 00988 /* Multiplication with the quadrature weights. */ 00989 /*fprintf(stderr,"\n");*/ 00990 d = 0; 00991 for (k = 0; k < m_theta; k++) 00992 { 00993 for (n = 0; n < m_phi; n++) 00994 { 00995 /*fprintf(stderr,"f_ref[%d] = %le + I*%le,\t f[%d] = %le + I*%le, \t w[%d] = %le\n", 00996 d,creal(f_ref[d]),cimag(f_ref[d]),d,creal(f[d]),cimag(f[d]),k, 00997 w[k]);*/ 00998 f_grid[d] *= w[k]; 00999 d++; 01000 } 01001 } 01002 01003 t1 = getticks(); 01004 t_avg += nfft_elapsed_seconds(t1,t0); 01005 01006 nfft_free(w); 01007 01008 t0 = getticks(); 01009 01010 /*fprintf(stderr,"\n"); 01011 d = 0; 01012 for (d = 0; d < grid_total; d++) 01013 { 01014 fprintf(stderr,"f[%d] = %le + I*%le, theta[%d] = %le, phi[%d] = %le\n", 01015 d,creal(f[d]),cimag(f[d]),d,x[2*d+1],d,x[2*d]); 01016 }*/ 01017 01018 //fprintf(stderr,"Executing adjoint\n"); 01019 //fflush(stderr); 01020 /* Check if the fast NFSFT algorithm shall be tested. */ 01021 if (use_nfsft != NO) 01022 { 01023 /* Execute the adjoint NFSFT transformation. */ 01024 nfsft_adjoint(plan_adjoint_ptr); 01025 } 01026 else 01027 { 01028 /* Execute the adjoint direct NDSFT transformation. */ 01029 nfsft_adjoint_direct(plan_adjoint_ptr); 01030 } 01031 01032 /* Multiplication with the Fourier-Legendre coefficients. */ 01033 /*for (k = 0; k <= m[im]; k++) 01034 { 01035 for (n = -k; n <= k; n++) 01036 { 01037 fprintf(stderr,"f_hat[%d,%d] = %le\t + I*%le\n",k,n, 01038 creal(f_hat[NFSFT_INDEX(k,n,&plan_adjoint)]), 01039 cimag(f_hat[NFSFT_INDEX(k,n,&plan_adjoint)])); 01040 } 01041 }*/ 01042 01043 //fprintf(stderr,"Executing trafo\n"); 01044 //fflush(stderr); 01045 if (use_nfsft != NO) 01046 { 01047 /* Execute the NFSFT transformation. */ 01048 nfsft_trafo(plan_ptr); 01049 } 01050 else 01051 { 01052 /* Execute the direct NDSFT transformation. */ 01053 nfsft_trafo_direct(plan_ptr); 01054 } 01055 01056 t1 = getticks(); 01057 t_avg += nfft_elapsed_seconds(t1,t0); 01058 01059 //fprintf(stderr,"Finalizing\n"); 01060 //fflush(stderr); 01061 /* Finalize the NFSFT plans */ 01062 nfsft_finalize(plan_adjoint_ptr); 01063 if (plan_ptr != plan_adjoint_ptr) 01064 { 01065 nfsft_finalize(plan_ptr); 01066 } 01067 01068 /* Free data arrays. */ 01069 nfft_free(f_hat); 01070 nfft_free(x_grid); 01071 01072 err_infty_avg += X(error_l_infty_complex)(f, f_compare, m_compare); 01073 err_2_avg += X(error_l_2_complex)(f, f_compare, m_compare); 01074 01075 nfft_free(f_grid); 01076 01077 if (mode == RANDOM) 01078 { 01079 } 01080 else 01081 { 01082 nfft_free(f_compare); 01083 } 01084 01085 /*for (d = 0; d < m_total; d++) 01086 { 01087 fprintf(stderr,"f_ref[%d] = %le + I*%le,\t f[%d] = %le + I*%le\n", 01088 d,creal(f_ref[d]),cimag(f_ref[d]),d,creal(f[d]),cimag(f[d])); 01089 }*/ 01090 } 01091 01092 //fprintf(stderr,"Calculating the error\n"); 01093 //fflush(stderr); 01094 /* Calculate average time needed. */ 01095 t_avg = t_avg/((double)repetitions); 01096 01097 /* Calculate the average error. */ 01098 err_infty_avg = err_infty_avg/((double)repetitions); 01099 01100 /* Calculate the average error. */ 01101 err_2_avg = err_2_avg/((double)repetitions); 01102 01103 /* Print out the error measurements. */ 01104 fprintf(stdout,"%+le %+le %+le\n", t_avg, err_infty_avg, err_2_avg); 01105 fprintf(stderr,"%d: %4d %4d %+le %+le %+le\n", tc, NQ[iNQ], SQ[iNQ], 01106 t_avg, err_infty_avg, err_2_avg); 01107 } 01108 } /* for (im = 0; im < im_max; im++) - Process all cut-off 01109 * bandwidths.*/ 01110 fprintf(stderr,"\n"); 01111 01112 /* Delete precomputed data. */ 01113 nfsft_forget(); 01114 01115 /* Free memory for cut-off bandwidths and grid size parameters. */ 01116 nfft_free(NQ); 01117 nfft_free(SQ); 01118 if (testmode == TIMING) 01119 { 01120 nfft_free(RQ); 01121 } 01122 01123 if (mode == RANDOM) 01124 { 01125 nfft_free(x_compare); 01126 nfft_free(f_compare); 01127 nfft_free(f); 01128 } 01129 01130 if (testmode == TIMING) 01131 { 01132 /* Allocate data structures. */ 01133 nfft_free(f_hat); 01134 nfft_free(f); 01135 nfft_free(x_grid); 01136 } 01137 01138 } /* for (tc = 0; tc < tc_max; tc++) - Process each testcase. */ 01139 01140 /* Return exit code for successful run. */ 01141 return EXIT_SUCCESS; 01142 } 01143 /* \} */