00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00027
00028 #include <math.h>
00029 #include <stdlib.h>
00030
00031 #include <stdio.h>
00032 #include <string.h>
00033 #include <time.h>
00034
00035 #include <complex.h>
00036
00037
00038 #include "nfft3util.h"
00039
00040
00041 #include "nfft3.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;
00114 double t_avg;
00115 double err_infty_avg;
00116 double err_2_avg;
00118 int i;
00119 int k;
00120 int n;
00121 int d;
00123 int m_theta;
00125 int m_phi;
00127 int m_total;
00128 double *theta;
00130 double *phi;
00132 fftw_plan fplan;
00134
00135 int d2;
00136 int M;
00137 double theta_s;
00138 double x1,x2,x3,temp;
00139 int m_compare;
00140 nfsft_plan *plan_adjoint_ptr;
00141 nfsft_plan *plan_ptr;
00142 double *w_temp;
00143 int testmode;
00144
00145
00146 fscanf(stdin,"testcases=%d\n",&tc_max);
00147 fprintf(stdout,"%d\n",tc_max);
00148
00149
00150 for (tc = 0; tc < tc_max; tc++)
00151 {
00152
00153 fscanf(stdin,"nfsft=%d\n",&use_nfsft);
00154 fprintf(stdout,"%d\n",use_nfsft);
00155 if (use_nfsft != NO)
00156 {
00157
00158 fscanf(stdin,"nfft=%d\n",&use_nfft);
00159 fprintf(stdout,"%d\n",use_nfsft);
00160 if (use_nfft != NO)
00161 {
00162
00163 fscanf(stdin,"cutoff=%d\n",&cutoff);
00164 fprintf(stdout,"%d\n",cutoff);
00165 }
00166 else
00167 {
00168
00169
00170 cutoff = 1;
00171 }
00172
00173 fscanf(stdin,"fpt=%d\n",&use_fpt);
00174 fprintf(stdout,"%d\n",use_fpt);
00175 if (use_fpt != NO)
00176 {
00177
00178 fscanf(stdin,"threshold=%lf\n",&threshold);
00179 fprintf(stdout,"%lf\n",threshold);
00180 }
00181 else
00182 {
00183
00184
00185 threshold = 1000.0;
00186 }
00187 }
00188 else
00189 {
00190
00191
00192 use_nfft = NO;
00193 use_fpt = NO;
00194 cutoff = 3;
00195 threshold = 1000.0;
00196 }
00197
00198
00199 fscanf(stdin,"testmode=%d\n",&testmode);
00200 fprintf(stdout,"%d\n",testmode);
00201
00202 if (testmode == ERROR)
00203 {
00204
00205 fscanf(stdin,"gridtype=%d\n",&gridtype);
00206 fprintf(stdout,"%d\n",gridtype);
00207
00208
00209 fscanf(stdin,"testfunction=%d\n",&testfunction);
00210 fprintf(stdout,"%d\n",testfunction);
00211
00212
00213 if (testfunction == FUNCTION_RANDOM_BANDLIMITED)
00214 {
00215
00216 fscanf(stdin,"bandlimit=%d\n",&N);
00217 fprintf(stdout,"%d\n",N);
00218 }
00219 else
00220 {
00221 N = 1;
00222 }
00223
00224
00225 fscanf(stdin,"repetitions=%d\n",&repetitions);
00226 fprintf(stdout,"%d\n",repetitions);
00227
00228 fscanf(stdin,"mode=%d\n",&mode);
00229 fprintf(stdout,"%d\n",mode);
00230
00231 if (mode == RANDOM)
00232 {
00233
00234 fscanf(stdin,"points=%d\n",&m_compare);
00235 fprintf(stdout,"%d\n",m_compare);
00236 x_compare = (double*) nfft_malloc(2*m_compare*sizeof(double));
00237 d = 0;
00238 while (d < m_compare)
00239 {
00240 x1 = 2.0*(((double)rand())/RAND_MAX) - 1.0;
00241 x2 = 2.0*(((double)rand())/RAND_MAX) - 1.0;
00242 x3 = 2.0*(((double)rand())/RAND_MAX) - 1.0;
00243 temp = sqrt(x1*x1+x2*x2+x3*x3);
00244 if (temp <= 1)
00245 {
00246 x_compare[2*d+1] = acos(x3);
00247 if (x_compare[2*d+1] == 0 || x_compare[2*d+1] == PI)
00248 {
00249 x_compare[2*d] = 0.0;
00250 }
00251 else
00252 {
00253 x_compare[2*d] = atan2(x2/sin(x_compare[2*d+1]),x1/sin(x_compare[2*d+1]));
00254 }
00255 x_compare[2*d] *= 1.0/(2.0*PI);
00256 x_compare[2*d+1] *= 1.0/(2.0*PI);
00257 d++;
00258 }
00259 }
00260 f_compare = (double _Complex*) nfft_malloc(m_compare*sizeof(double _Complex));
00261 f = (double _Complex*) nfft_malloc(m_compare*sizeof(double _Complex));
00262 }
00263 }
00264
00265
00266 NQ_max = 0;
00267 SQ_max = 0;
00268
00269
00270 fscanf(stdin,"bandwidths=%d\n",&iNQ_max);
00271 fprintf(stdout,"%d\n",iNQ_max);
00272
00273
00274 NQ = (int*) nfft_malloc(iNQ_max*sizeof(int));
00275 SQ = (int*) nfft_malloc(iNQ_max*sizeof(int));
00276 if (testmode == TIMING)
00277 {
00278 RQ = (int*) nfft_malloc(iNQ_max*sizeof(int));
00279 }
00280
00281
00282 for (iNQ = 0; iNQ < iNQ_max; iNQ++)
00283 {
00284 if (testmode == TIMING)
00285 {
00286
00287 fscanf(stdin,"%d %d %d\n",&NQ[iNQ],&SQ[iNQ],&RQ[iNQ]);
00288 fprintf(stdout,"%d %d %d\n",NQ[iNQ],SQ[iNQ],RQ[iNQ]);
00289 NQ_max = NFFT_MAX(NQ_max,NQ[iNQ]);
00290 SQ_max = NFFT_MAX(SQ_max,SQ[iNQ]);
00291 }
00292 else
00293 {
00294
00295 fscanf(stdin,"%d %d\n",&NQ[iNQ],&SQ[iNQ]);
00296 fprintf(stdout,"%d %d\n",NQ[iNQ],SQ[iNQ]);
00297 NQ_max = NFFT_MAX(NQ_max,NQ[iNQ]);
00298 SQ_max = NFFT_MAX(SQ_max,SQ[iNQ]);
00299 }
00300 }
00301
00302
00303
00304
00305 nfsft_precompute(NQ_max, threshold,
00306 ((use_nfsft==NO)?(NFSFT_NO_FAST_ALGORITHM):(0U)), 0U);
00307
00308 if (testmode == TIMING)
00309 {
00310
00311 f_hat = (double _Complex*) nfft_malloc(NFSFT_F_HAT_SIZE(NQ_max)*sizeof(double _Complex));
00312 f = (double _Complex*) nfft_malloc(SQ_max*sizeof(double _Complex));
00313 x_grid = (double*) nfft_malloc(2*SQ_max*sizeof(double));
00314 for (d = 0; d < SQ_max; d++)
00315 {
00316 f[d] = (((double)rand())/RAND_MAX)-0.5 + _Complex_I*((((double)rand())/RAND_MAX)-0.5);
00317 x_grid[2*d] = (((double)rand())/RAND_MAX) - 0.5;
00318 x_grid[2*d+1] = (((double)rand())/RAND_MAX) * 0.5;
00319 }
00320 }
00321
00322
00323
00324
00325 for (iNQ = 0; iNQ < iNQ_max; iNQ++)
00326 {
00327 if (testmode == TIMING)
00328 {
00329 nfsft_init_guru(&plan,NQ[iNQ],SQ[iNQ], NFSFT_NORMALIZED |
00330 ((use_nfft!=NO)?(0U):(NFSFT_USE_NDFT)) |
00331 ((use_fpt!=NO)?(0U):(NFSFT_USE_DPT)),
00332 PRE_PHI_HUT | PRE_PSI | FFTW_INIT | FFTW_MEASURE | FFT_OUT_OF_PLACE,
00333 cutoff);
00334
00335 plan.f_hat = f_hat;
00336 plan.x = x_grid;
00337 plan.f = f;
00338
00339 nfsft_precompute_x(&plan);
00340
00341 t_avg = 0.0;
00342
00343 for (i = 0; i < RQ[iNQ]; i++)
00344 {
00345 t = nfft_second();
00346
00347 if (use_nfsft != NO)
00348 {
00349
00350 nfsft_adjoint(&plan);
00351 }
00352 else
00353 {
00354
00355 ndsft_adjoint(&plan);
00356 }
00357
00358 t_avg += nfft_second() - t;
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
00371 switch (gridtype)
00372 {
00373 case GRID_GAUSS_LEGENDRE:
00374
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
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
00390 break;
00391 case GRID_EQUIDISTRIBUTION:
00392 case GRID_EQUIDISTRIBUTION_UNIFORM:
00393 m_theta = 2;
00394
00395 for (k = 1; k < SQ[iNQ]; k++)
00396 {
00397 m_theta += (int)floor((2*PI)/acos((cos(PI/(double)SQ[iNQ])-
00398 cos(k*PI/(double)SQ[iNQ])*cos(k*PI/(double)SQ[iNQ]))/
00399 (sin(k*PI/(double)SQ[iNQ])*sin(k*PI/(double)SQ[iNQ]))));
00400
00401 }
00402
00403 m_phi = 1;
00404 m_total = m_theta * m_phi;
00405 break;
00406 }
00407
00408
00409 w = (double*) nfft_malloc(m_theta*sizeof(double));
00410 x_grid = (double*) nfft_malloc(2*m_total*sizeof(double));
00411
00412
00413
00414 switch (gridtype)
00415 {
00416 case GRID_GAUSS_LEGENDRE:
00417
00418
00419
00420
00421 for (k = 0; k < m_theta; k++)
00422 {
00423 fscanf(stdin,"%le\n",&w[k]);
00424 w[k] *= (2.0*PI)/((double)m_phi);
00425 }
00426
00427
00428
00429
00430 theta = (double*) nfft_malloc(m_theta*sizeof(double));
00431 phi = (double*) nfft_malloc(m_phi*sizeof(double));
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441 for (k = 0; k < m_theta; k++)
00442 {
00443 fscanf(stdin,"%le\n",&theta[k]);
00444 }
00445
00446
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
00454
00455
00456
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
00469
00470
00471 nfft_free(theta);
00472 nfft_free(phi);
00473
00474 break;
00475
00476 case GRID_CLENSHAW_CURTIS:
00477
00478
00479 theta = (double*) nfft_malloc(m_theta*sizeof(double));
00480 phi = (double*) nfft_malloc(m_phi*sizeof(double));
00481
00482
00483 for (k = 0; k < m_theta; k++)
00484 {
00485 theta[k] = k/((double)2*(m_theta-1));
00486 }
00487
00488
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
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*PI)/((double)(m_theta-1)*m_phi);
00507 w[m_theta-1-k] = w[k];
00508 }
00509 fftw_destroy_plan(fplan);
00510
00511
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
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*PI);
00570 }
00571
00572 w[0] = (4.0*PI)/(m_total);
00573 break;
00574
00575 case GRID_EQUIDISTRIBUTION:
00576 case GRID_EQUIDISTRIBUTION_UNIFORM:
00577
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*PI)/((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*PI)/(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*PI)/(m_total);
00619 }
00620 d = 2;
00621
00622 for (k = 1; k < SQ[iNQ]; k++)
00623 {
00624 theta_s = (double)k*PI/(double)SQ[iNQ];
00625 M = (int)floor((2.0*PI)/acos((cos(PI/(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*PI);
00633 if (gridtype == GRID_EQUIDISTRIBUTION)
00634 {
00635 w[d] = w_temp[k]/((double)(M));
00636 }
00637 else
00638 {
00639 w[d] = (4.0*PI)/(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
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
00670
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
00676
00677
00678
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
00708 nfsft_trafo(&plan_gen);
00709 }
00710 else
00711 {
00712
00713 ndsft_trafo(&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
00735 nfsft_trafo(&plan_gen);
00736 }
00737 else
00738 {
00739
00740 ndsft_trafo(&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*PI)*cos(x_grid[2*d]*2.0*PI);
00758 x2 = sin(x_grid[2*d+1]*2.0*PI)*sin(x_grid[2*d]*2.0*PI);
00759 x3 = cos(x_grid[2*d+1]*2.0*PI);
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*PI)*cos(x_compare[2*d]*2.0*PI);
00767 x2 = sin(x_compare[2*d+1]*2.0*PI)*sin(x_compare[2*d]*2.0*PI);
00768 x3 = cos(x_compare[2*d+1]*2.0*PI);
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*PI)*cos(x_grid[2*d]*2.0*PI);
00781 x2 = sin(x_grid[2*d+1]*2.0*PI)*sin(x_grid[2*d]*2.0*PI);
00782 x3 = cos(x_grid[2*d+1]*2.0*PI);
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*PI)*cos(x_compare[2*d]*2.0*PI);
00790 x2 = sin(x_compare[2*d+1]*2.0*PI)*sin(x_compare[2*d]*2.0*PI);
00791 x3 = cos(x_compare[2*d+1]*2.0*PI);
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*PI)*cos(x_grid[2*d]*2.0*PI);
00804 x2 = sin(x_grid[2*d+1]*2.0*PI)*sin(x_grid[2*d]*2.0*PI);
00805 x3 = cos(x_grid[2*d+1]*2.0*PI);
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*PI)*cos(x_compare[2*d]*2.0*PI);
00814 x2 = sin(x_compare[2*d+1]*2.0*PI)*sin(x_compare[2*d]*2.0*PI);
00815 x3 = cos(x_compare[2*d+1]*2.0*PI);
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*PI)*cos(x_grid[2*d]*2.0*PI);
00829 x2 = sin(x_grid[2*d+1]*2.0*PI)*sin(x_grid[2*d]*2.0*PI);
00830 x3 = cos(x_grid[2*d+1]*2.0*PI);
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*PI)*cos(x_compare[2*d]*2.0*PI);
00839 x2 = sin(x_compare[2*d+1]*2.0*PI)*sin(x_compare[2*d]*2.0*PI);
00840 x3 = cos(x_compare[2*d+1]*2.0*PI);
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*PI)*cos(x_grid[2*d]*2.0*PI);
00854 x2 = sin(x_grid[2*d+1]*2.0*PI)*sin(x_grid[2*d]*2.0*PI);
00855 x3 = cos(x_grid[2*d+1]*2.0*PI);
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*PI)*cos(x_compare[2*d]*2.0*PI);
00864 x2 = sin(x_compare[2*d+1]*2.0*PI)*sin(x_compare[2*d]*2.0*PI);
00865 x3 = cos(x_compare[2*d+1]*2.0*PI);
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*PI*x_grid[2*d+1])*cos(2.0*PI*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*PI*x_compare[2*d+1])*cos(2.0*PI*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
00908
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
00928
00929
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
00963
00964 nfsft_precompute_x(plan_adjoint_ptr);
00965 if (plan_adjoint_ptr != plan_ptr)
00966 {
00967 nfsft_precompute_x(plan_ptr);
00968 }
00969
00970
00971 t_avg = 0.0;
00972 err_infty_avg = 0.0;
00973 err_2_avg = 0.0;
00974
00975
00976 for (i = 0; i < 1; i++)
00977 {
00978
00979
00980
00981
00982
00983
00984 t = nfft_second();
00985
00986
00987
00988
00989
00990 d = 0;
00991 for (k = 0; k < m_theta; k++)
00992 {
00993 for (n = 0; n < m_phi; n++)
00994 {
00995
00996
00997
00998 f_grid[d] *= w[k];
00999 d++;
01000 }
01001 }
01002
01003 t_avg += nfft_second() - t;
01004
01005 nfft_free(w);
01006
01007 t = nfft_second();
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020 if (use_nfsft != NO)
01021 {
01022
01023 nfsft_adjoint(plan_adjoint_ptr);
01024 }
01025 else
01026 {
01027
01028 ndsft_adjoint(plan_adjoint_ptr);
01029 }
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044 if (use_nfsft != NO)
01045 {
01046
01047 nfsft_trafo(plan_ptr);
01048 }
01049 else
01050 {
01051
01052 ndsft_trafo(plan_ptr);
01053 }
01054
01055 t_avg += nfft_second() - t;
01056
01057
01058
01059
01060 nfsft_finalize(plan_adjoint_ptr);
01061 if (plan_ptr != plan_adjoint_ptr)
01062 {
01063 nfsft_finalize(plan_ptr);
01064 }
01065
01066
01067 nfft_free(f_hat);
01068 nfft_free(x_grid);
01069
01070 err_infty_avg += nfft_error_l_infty_complex(f, f_compare, m_compare);
01071 err_2_avg += nfft_error_l_2_complex(f, f_compare, m_compare);
01072
01073 nfft_free(f_grid);
01074
01075 if (mode == RANDOM)
01076 {
01077 }
01078 else
01079 {
01080 nfft_free(f_compare);
01081 }
01082
01083
01084
01085
01086
01087
01088 }
01089
01090
01091
01092
01093 t_avg = t_avg/((double)repetitions);
01094
01095
01096 err_infty_avg = err_infty_avg/((double)repetitions);
01097
01098
01099 err_2_avg = err_2_avg/((double)repetitions);
01100
01101
01102 fprintf(stdout,"%+le %+le %+le\n", t_avg, err_infty_avg, err_2_avg);
01103 fprintf(stderr,"%d: %4d %4d %+le %+le %+le\n", tc, NQ[iNQ], SQ[iNQ],
01104 t_avg, err_infty_avg, err_2_avg);
01105 }
01106 }
01107
01108 fprintf(stderr,"\n");
01109
01110
01111 nfsft_forget();
01112
01113
01114 nfft_free(NQ);
01115 nfft_free(SQ);
01116 if (testmode == TIMING)
01117 {
01118 nfft_free(RQ);
01119 }
01120
01121 if (mode == RANDOM)
01122 {
01123 nfft_free(x_compare);
01124 nfft_free(f_compare);
01125 nfft_free(f);
01126 }
01127
01128 if (testmode == TIMING)
01129 {
01130
01131 nfft_free(f_hat);
01132 nfft_free(f);
01133 nfft_free(x_grid);
01134 }
01135
01136 }
01137
01138
01139 return EXIT_SUCCESS;
01140 }
01141