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