00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00026 #include "config.h"
00027
00028
00029 #include <math.h>
00030 #include <stdlib.h>
00031 #include <stdio.h>
00032 #include <string.h>
00033 #include <time.h>
00034 #ifdef HAVE_COMPLEX_H
00035 #include <complex.h>
00036 #endif
00037
00038
00039 #include "nfft3util.h"
00040
00041
00042 #include "nfft3.h"
00043
00044 #include "infft.h"
00045
00047 enum boolean {NO = 0, YES = 1};
00048
00050 enum testtype {ERROR = 0, TIMING = 1};
00051
00053 enum gridtype {GRID_GAUSS_LEGENDRE = 0, GRID_CLENSHAW_CURTIS = 1,
00054 GRID_HEALPIX = 2, GRID_EQUIDISTRIBUTION = 3, GRID_EQUIDISTRIBUTION_UNIFORM = 4};
00055
00057 enum functiontype {FUNCTION_RANDOM_BANDLIMITED = 0, FUNCTION_F1 = 1,
00058 FUNCTION_F2 = 2, FUNCTION_F3 = 3, FUNCTION_F4 = 4, FUNCTION_F5 = 5,
00059 FUNCTION_F6 = 6};
00060
00062 enum modes {USE_GRID = 0, RANDOM = 1};
00063
00072 int main (int argc, char **argv)
00073 {
00074 int tc;
00075 int tc_max;
00077 int *NQ;
00079 int NQ_max;
00081 int *SQ;
00083 int SQ_max;
00084 int *RQ;
00086 int iNQ;
00087 int iNQ_max;
00088 int testfunction;
00089 int N;
00091 int use_nfsft;
00092 int use_nfft;
00093 int use_fpt;
00094 int cutoff;
00095 double threshold;
00097 int gridtype;
00098 int repetitions;
00099 int mode;
00101 double *w;
00102 double *x_grid;
00103 double *x_compare;
00104 double _Complex *f_grid;
00105 double _Complex *f_compare;
00106 double _Complex *f;
00107 double _Complex *f_hat_gen;
00109 double _Complex *f_hat;
00111 nfsft_plan plan_adjoint;
00112 nfsft_plan plan;
00113 nfsft_plan plan_gen;
00115 double t_avg;
00116 double err_infty_avg;
00117 double err_2_avg;
00119 int i;
00120 int k;
00121 int n;
00122 int d;
00124 int m_theta;
00126 int m_phi;
00128 int m_total;
00129 double *theta;
00131 double *phi;
00133 fftw_plan fplan;
00135
00136 int d2;
00137 int M;
00138 double theta_s;
00139 double x1,x2,x3,temp;
00140 int m_compare;
00141 nfsft_plan *plan_adjoint_ptr;
00142 nfsft_plan *plan_ptr;
00143 double *w_temp;
00144 int testmode;
00145 ticks t0, t1;
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*) nfft_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*) nfft_malloc(m_compare*sizeof(double _Complex));
00263 f = (double _Complex*) nfft_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*) nfft_malloc(iNQ_max*sizeof(int));
00277 SQ = (int*) nfft_malloc(iNQ_max*sizeof(int));
00278 if (testmode == TIMING)
00279 {
00280 RQ = (int*) nfft_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*) nfft_malloc(NFSFT_F_HAT_SIZE(NQ_max)*sizeof(double _Complex));
00314 f = (double _Complex*) nfft_malloc(SQ_max*sizeof(double _Complex));
00315 x_grid = (double*) nfft_malloc(2*SQ_max*sizeof(double));
00316 for (d = 0; d < SQ_max; d++)
00317 {
00318 f[d] = (((double)rand())/RAND_MAX)-0.5 + _Complex_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 t0 = getticks();
00348
00349 if (use_nfsft != NO)
00350 {
00351
00352 nfsft_adjoint(&plan);
00353 }
00354 else
00355 {
00356
00357 nfsft_adjoint_direct(&plan);
00358 }
00359
00360 t1 = getticks();
00361 t_avg += nfft_elapsed_seconds(t1,t0);
00362 }
00363
00364 t_avg = t_avg/((double)RQ[iNQ]);
00365
00366 nfsft_finalize(&plan);
00367
00368 fprintf(stdout,"%+le\n", t_avg);
00369 fprintf(stderr,"%d: %4d %4d %+le\n", tc, NQ[iNQ], SQ[iNQ], t_avg);
00370 }
00371 else
00372 {
00373
00374 switch (gridtype)
00375 {
00376 case GRID_GAUSS_LEGENDRE:
00377
00378 m_theta = SQ[iNQ] + 1;
00379 m_phi = 2*SQ[iNQ] + 2;
00380 m_total = m_theta*m_phi;
00381 break;
00382 case GRID_CLENSHAW_CURTIS:
00383
00384 m_theta = 2*SQ[iNQ] + 1;
00385 m_phi = 2*SQ[iNQ] + 2;
00386 m_total = m_theta*m_phi;
00387 break;
00388 case GRID_HEALPIX:
00389 m_theta = 1;
00390 m_phi = 12*SQ[iNQ]*SQ[iNQ];
00391 m_total = m_theta * m_phi;
00392
00393 break;
00394 case GRID_EQUIDISTRIBUTION:
00395 case GRID_EQUIDISTRIBUTION_UNIFORM:
00396 m_theta = 2;
00397
00398 for (k = 1; k < SQ[iNQ]; k++)
00399 {
00400 m_theta += (int)floor((2*PI)/acos((cos(PI/(double)SQ[iNQ])-
00401 cos(k*PI/(double)SQ[iNQ])*cos(k*PI/(double)SQ[iNQ]))/
00402 (sin(k*PI/(double)SQ[iNQ])*sin(k*PI/(double)SQ[iNQ]))));
00403
00404 }
00405
00406 m_phi = 1;
00407 m_total = m_theta * m_phi;
00408 break;
00409 }
00410
00411
00412 w = (double*) nfft_malloc(m_theta*sizeof(double));
00413 x_grid = (double*) nfft_malloc(2*m_total*sizeof(double));
00414
00415
00416
00417 switch (gridtype)
00418 {
00419 case GRID_GAUSS_LEGENDRE:
00420
00421
00422
00423
00424 for (k = 0; k < m_theta; k++)
00425 {
00426 fscanf(stdin,"%le\n",&w[k]);
00427 w[k] *= (2.0*PI)/((double)m_phi);
00428 }
00429
00430
00431
00432
00433 theta = (double*) nfft_malloc(m_theta*sizeof(double));
00434 phi = (double*) nfft_malloc(m_phi*sizeof(double));
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444 for (k = 0; k < m_theta; k++)
00445 {
00446 fscanf(stdin,"%le\n",&theta[k]);
00447 }
00448
00449
00450 for (n = 0; n < m_phi; n++)
00451 {
00452 phi[n] = n/((double)m_phi);
00453 phi[n] -= ((phi[n]>=0.5)?(1.0):(0.0));
00454 }
00455
00456
00457
00458
00459
00460 d = 0;
00461 for (k = 0; k < m_theta; k++)
00462 {
00463 for (n = 0; n < m_phi; n++)
00464 {
00465 x_grid[2*d] = phi[n];
00466 x_grid[2*d+1] = theta[k];
00467 d++;
00468 }
00469 }
00470
00471
00472
00473
00474 nfft_free(theta);
00475 nfft_free(phi);
00476
00477 break;
00478
00479 case GRID_CLENSHAW_CURTIS:
00480
00481
00482 theta = (double*) nfft_malloc(m_theta*sizeof(double));
00483 phi = (double*) nfft_malloc(m_phi*sizeof(double));
00484
00485
00486 for (k = 0; k < m_theta; k++)
00487 {
00488 theta[k] = k/((double)2*(m_theta-1));
00489 }
00490
00491
00492 for (n = 0; n < m_phi; n++)
00493 {
00494 phi[n] = n/((double)m_phi);
00495 phi[n] -= ((phi[n]>=0.5)?(1.0):(0.0));
00496 }
00497
00498
00499 fplan = fftw_plan_r2r_1d(SQ[iNQ]+1, w, w, FFTW_REDFT00, 0U);
00500 for (k = 0; k < SQ[iNQ]+1; k++)
00501 {
00502 w[k] = -2.0/(4*k*k-1);
00503 }
00504 fftw_execute(fplan);
00505 w[0] *= 0.5;
00506
00507 for (k = 0; k < SQ[iNQ]+1; k++)
00508 {
00509 w[k] *= (2.0*PI)/((double)(m_theta-1)*m_phi);
00510 w[m_theta-1-k] = w[k];
00511 }
00512 fftw_destroy_plan(fplan);
00513
00514
00515 d = 0;
00516 for (k = 0; k < m_theta; k++)
00517 {
00518 for (n = 0; n < m_phi; n++)
00519 {
00520 x_grid[2*d] = phi[n];
00521 x_grid[2*d+1] = theta[k];
00522 d++;
00523 }
00524 }
00525
00526
00527 nfft_free(theta);
00528 nfft_free(phi);
00529
00530 break;
00531
00532 case GRID_HEALPIX:
00533
00534 d = 0;
00535 for (k = 1; k <= SQ[iNQ]-1; k++)
00536 {
00537 for (n = 0; n <= 4*k-1; n++)
00538 {
00539 x_grid[2*d+1] = 1 - (k*k)/((double)(3.0*SQ[iNQ]*SQ[iNQ]));
00540 x_grid[2*d] = ((n+0.5)/(4*k));
00541 x_grid[2*d] -= (x_grid[2*d]>=0.5)?(1.0):(0.0);
00542 d++;
00543 }
00544 }
00545
00546 d2 = d-1;
00547
00548 for (k = SQ[iNQ]; k <= 3*SQ[iNQ]; k++)
00549 {
00550 for (n = 0; n <= 4*SQ[iNQ]-1; n++)
00551 {
00552 x_grid[2*d+1] = 2.0/(3*SQ[iNQ])*(2*SQ[iNQ]-k);
00553 x_grid[2*d] = (n+((k%2==0)?(0.5):(0.0)))/(4*SQ[iNQ]);
00554 x_grid[2*d] -= (x_grid[2*d]>=0.5)?(1.0):(0.0);
00555 d++;
00556 }
00557 }
00558
00559 for (k = 1; k <= SQ[iNQ]-1; k++)
00560 {
00561 for (n = 0; n <= 4*k-1; n++)
00562 {
00563 x_grid[2*d+1] = -x_grid[2*d2+1];
00564 x_grid[2*d] = x_grid[2*d2];
00565 d++;
00566 d2--;
00567 }
00568 }
00569
00570 for (d = 0; d < m_total; d++)
00571 {
00572 x_grid[2*d+1] = acos(x_grid[2*d+1])/(2.0*PI);
00573 }
00574
00575 w[0] = (4.0*PI)/(m_total);
00576 break;
00577
00578 case GRID_EQUIDISTRIBUTION:
00579 case GRID_EQUIDISTRIBUTION_UNIFORM:
00580
00581
00582 if (gridtype == GRID_EQUIDISTRIBUTION)
00583 {
00584 w_temp = (double*) nfft_malloc((SQ[iNQ]+1)*sizeof(double));
00585 fplan = fftw_plan_r2r_1d(SQ[iNQ]/2+1, w_temp, w_temp, FFTW_REDFT00, 0U);
00586 for (k = 0; k < SQ[iNQ]/2+1; k++)
00587 {
00588 w_temp[k] = -2.0/(4*k*k-1);
00589 }
00590 fftw_execute(fplan);
00591 w_temp[0] *= 0.5;
00592
00593 for (k = 0; k < SQ[iNQ]/2+1; k++)
00594 {
00595 w_temp[k] *= (2.0*PI)/((double)(SQ[iNQ]));
00596 w_temp[SQ[iNQ]-k] = w_temp[k];
00597 }
00598 fftw_destroy_plan(fplan);
00599 }
00600
00601 d = 0;
00602 x_grid[2*d] = -0.5;
00603 x_grid[2*d+1] = 0.0;
00604 if (gridtype == GRID_EQUIDISTRIBUTION)
00605 {
00606 w[d] = w_temp[0];
00607 }
00608 else
00609 {
00610 w[d] = (4.0*PI)/(m_total);
00611 }
00612 d = 1;
00613 x_grid[2*d] = -0.5;
00614 x_grid[2*d+1] = 0.5;
00615 if (gridtype == GRID_EQUIDISTRIBUTION)
00616 {
00617 w[d] = w_temp[SQ[iNQ]];
00618 }
00619 else
00620 {
00621 w[d] = (4.0*PI)/(m_total);
00622 }
00623 d = 2;
00624
00625 for (k = 1; k < SQ[iNQ]; k++)
00626 {
00627 theta_s = (double)k*PI/(double)SQ[iNQ];
00628 M = (int)floor((2.0*PI)/acos((cos(PI/(double)SQ[iNQ])-
00629 cos(theta_s)*cos(theta_s))/(sin(theta_s)*sin(theta_s))));
00630
00631 for (n = 0; n < M; n++)
00632 {
00633 x_grid[2*d] = (n + 0.5)/M;
00634 x_grid[2*d] -= (x_grid[2*d]>=0.5)?(1.0):(0.0);
00635 x_grid[2*d+1] = theta_s/(2.0*PI);
00636 if (gridtype == GRID_EQUIDISTRIBUTION)
00637 {
00638 w[d] = w_temp[k]/((double)(M));
00639 }
00640 else
00641 {
00642 w[d] = (4.0*PI)/(m_total);
00643 }
00644 d++;
00645 }
00646 }
00647
00648 if (gridtype == GRID_EQUIDISTRIBUTION)
00649 {
00650 nfft_free(w_temp);
00651 }
00652 break;
00653
00654 default:
00655 break;
00656 }
00657
00658
00659 f_grid = (double _Complex*) nfft_malloc(m_total*sizeof(double _Complex));
00660
00661 if (mode == RANDOM)
00662 {
00663 }
00664 else
00665 {
00666 m_compare = m_total;
00667 f_compare = (double _Complex*) nfft_malloc(m_compare*sizeof(double _Complex));
00668 x_compare = x_grid;
00669 f = f_grid;
00670 }
00671
00672
00673
00674 switch (testfunction)
00675 {
00676 case FUNCTION_RANDOM_BANDLIMITED:
00677 f_hat_gen = (double _Complex*) nfft_malloc(NFSFT_F_HAT_SIZE(N)*sizeof(double _Complex));
00678
00679
00680
00681
00682 nfsft_init_guru(&plan_gen,N,m_total, NFSFT_NORMALIZED |
00683 ((use_nfft!=NO)?(0U):(NFSFT_USE_NDFT)) |
00684 ((use_fpt!=NO)?(0U):(NFSFT_USE_DPT)),
00685 ((N>512)?(0U):(PRE_PHI_HUT | PRE_PSI)) | FFTW_INIT |
00686 FFT_OUT_OF_PLACE, cutoff);
00687
00688 plan_gen.f_hat = f_hat_gen;
00689 plan_gen.x = x_grid;
00690 plan_gen.f = f_grid;
00691
00692 nfsft_precompute_x(&plan_gen);
00693
00694 for (k = 0; k < plan_gen.N_total; k++)
00695 {
00696 f_hat_gen[k] = 0.0;
00697 }
00698
00699 for (k = 0; k <= N; k++)
00700 {
00701 for (n = -k; n <= k; n++)
00702 {
00703 f_hat_gen[NFSFT_INDEX(k,n,&plan_gen)] =
00704 (((double)rand())/RAND_MAX)-0.5 + _Complex_I*((((double)rand())/RAND_MAX)-0.5);
00705 }
00706 }
00707
00708 if (use_nfsft != NO)
00709 {
00710
00711 nfsft_trafo(&plan_gen);
00712 }
00713 else
00714 {
00715
00716 nfsft_trafo_direct(&plan_gen);
00717 }
00718
00719 nfsft_finalize(&plan_gen);
00720
00721 if (mode == RANDOM)
00722 {
00723 nfsft_init_guru(&plan_gen,N,m_compare, NFSFT_NORMALIZED |
00724 ((use_nfft!=NO)?(0U):(NFSFT_USE_NDFT)) |
00725 ((use_fpt!=NO)?(0U):(NFSFT_USE_DPT)),
00726 ((N>512)?(0U):(PRE_PHI_HUT | PRE_PSI)) | FFTW_INIT |
00727 FFT_OUT_OF_PLACE, cutoff);
00728
00729 plan_gen.f_hat = f_hat_gen;
00730 plan_gen.x = x_compare;
00731 plan_gen.f = f_compare;
00732
00733 nfsft_precompute_x(&plan_gen);
00734
00735 if (use_nfsft != NO)
00736 {
00737
00738 nfsft_trafo(&plan_gen);
00739 }
00740 else
00741 {
00742
00743 nfsft_trafo_direct(&plan_gen);
00744 }
00745
00746 nfsft_finalize(&plan_gen);
00747 }
00748 else
00749 {
00750 memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
00751 }
00752
00753 nfft_free(f_hat_gen);
00754
00755 break;
00756
00757 case FUNCTION_F1:
00758 for (d = 0; d < m_total; d++)
00759 {
00760 x1 = sin(x_grid[2*d+1]*2.0*PI)*cos(x_grid[2*d]*2.0*PI);
00761 x2 = sin(x_grid[2*d+1]*2.0*PI)*sin(x_grid[2*d]*2.0*PI);
00762 x3 = cos(x_grid[2*d+1]*2.0*PI);
00763 f_grid[d] = x1*x2*x3;
00764 }
00765 if (mode == RANDOM)
00766 {
00767 for (d = 0; d < m_compare; d++)
00768 {
00769 x1 = sin(x_compare[2*d+1]*2.0*PI)*cos(x_compare[2*d]*2.0*PI);
00770 x2 = sin(x_compare[2*d+1]*2.0*PI)*sin(x_compare[2*d]*2.0*PI);
00771 x3 = cos(x_compare[2*d+1]*2.0*PI);
00772 f_compare[d] = x1*x2*x3;
00773 }
00774 }
00775 else
00776 {
00777 memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
00778 }
00779 break;
00780 case FUNCTION_F2:
00781 for (d = 0; d < m_total; d++)
00782 {
00783 x1 = sin(x_grid[2*d+1]*2.0*PI)*cos(x_grid[2*d]*2.0*PI);
00784 x2 = sin(x_grid[2*d+1]*2.0*PI)*sin(x_grid[2*d]*2.0*PI);
00785 x3 = cos(x_grid[2*d+1]*2.0*PI);
00786 f_grid[d] = 0.1*exp(x1+x2+x3);
00787 }
00788 if (mode == RANDOM)
00789 {
00790 for (d = 0; d < m_compare; d++)
00791 {
00792 x1 = sin(x_compare[2*d+1]*2.0*PI)*cos(x_compare[2*d]*2.0*PI);
00793 x2 = sin(x_compare[2*d+1]*2.0*PI)*sin(x_compare[2*d]*2.0*PI);
00794 x3 = cos(x_compare[2*d+1]*2.0*PI);
00795 f_compare[d] = 0.1*exp(x1+x2+x3);
00796 }
00797 }
00798 else
00799 {
00800 memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
00801 }
00802 break;
00803 case FUNCTION_F3:
00804 for (d = 0; d < m_total; d++)
00805 {
00806 x1 = sin(x_grid[2*d+1]*2.0*PI)*cos(x_grid[2*d]*2.0*PI);
00807 x2 = sin(x_grid[2*d+1]*2.0*PI)*sin(x_grid[2*d]*2.0*PI);
00808 x3 = cos(x_grid[2*d+1]*2.0*PI);
00809 temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3);
00810 f_grid[d] = 0.1*temp;
00811 }
00812 if (mode == RANDOM)
00813 {
00814 for (d = 0; d < m_compare; d++)
00815 {
00816 x1 = sin(x_compare[2*d+1]*2.0*PI)*cos(x_compare[2*d]*2.0*PI);
00817 x2 = sin(x_compare[2*d+1]*2.0*PI)*sin(x_compare[2*d]*2.0*PI);
00818 x3 = cos(x_compare[2*d+1]*2.0*PI);
00819 temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3);
00820 f_compare[d] = 0.1*temp;
00821 }
00822 }
00823 else
00824 {
00825 memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
00826 }
00827 break;
00828 case FUNCTION_F4:
00829 for (d = 0; d < m_total; d++)
00830 {
00831 x1 = sin(x_grid[2*d+1]*2.0*PI)*cos(x_grid[2*d]*2.0*PI);
00832 x2 = sin(x_grid[2*d+1]*2.0*PI)*sin(x_grid[2*d]*2.0*PI);
00833 x3 = cos(x_grid[2*d+1]*2.0*PI);
00834 temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3);
00835 f_grid[d] = 1.0/(temp);
00836 }
00837 if (mode == RANDOM)
00838 {
00839 for (d = 0; d < m_compare; d++)
00840 {
00841 x1 = sin(x_compare[2*d+1]*2.0*PI)*cos(x_compare[2*d]*2.0*PI);
00842 x2 = sin(x_compare[2*d+1]*2.0*PI)*sin(x_compare[2*d]*2.0*PI);
00843 x3 = cos(x_compare[2*d+1]*2.0*PI);
00844 temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3);
00845 f_compare[d] = 1.0/(temp);
00846 }
00847 }
00848 else
00849 {
00850 memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
00851 }
00852 break;
00853 case FUNCTION_F5:
00854 for (d = 0; d < m_total; d++)
00855 {
00856 x1 = sin(x_grid[2*d+1]*2.0*PI)*cos(x_grid[2*d]*2.0*PI);
00857 x2 = sin(x_grid[2*d+1]*2.0*PI)*sin(x_grid[2*d]*2.0*PI);
00858 x3 = cos(x_grid[2*d+1]*2.0*PI);
00859 temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3);
00860 f_grid[d] = 0.1*sin(1+temp)*sin(1+temp);
00861 }
00862 if (mode == RANDOM)
00863 {
00864 for (d = 0; d < m_compare; d++)
00865 {
00866 x1 = sin(x_compare[2*d+1]*2.0*PI)*cos(x_compare[2*d]*2.0*PI);
00867 x2 = sin(x_compare[2*d+1]*2.0*PI)*sin(x_compare[2*d]*2.0*PI);
00868 x3 = cos(x_compare[2*d+1]*2.0*PI);
00869 temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3);
00870 f_compare[d] = 0.1*sin(1+temp)*sin(1+temp);
00871 }
00872 }
00873 else
00874 {
00875 memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
00876 }
00877 break;
00878 case FUNCTION_F6:
00879 for (d = 0; d < m_total; d++)
00880 {
00881 if (x_grid[2*d+1] <= 0.25)
00882 {
00883 f_grid[d] = 1.0;
00884 }
00885 else
00886 {
00887 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])));
00888 }
00889 }
00890 if (mode == RANDOM)
00891 {
00892 for (d = 0; d < m_compare; d++)
00893 {
00894 if (x_compare[2*d+1] <= 0.25)
00895 {
00896 f_compare[d] = 1.0;
00897 }
00898 else
00899 {
00900 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])));
00901 }
00902 }
00903 }
00904 else
00905 {
00906 memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
00907 }
00908 break;
00909 default:
00910
00911
00912 for (d = 0; d < m_total; d++)
00913 {
00914 f_grid[d] = 1.0;
00915 }
00916 if (mode == RANDOM)
00917 {
00918 for (d = 0; d < m_compare; d++)
00919 {
00920 f_compare[d] = 1.0;
00921 }
00922 }
00923 else
00924 {
00925 memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
00926 }
00927 break;
00928 }
00929
00930
00931
00932
00933 nfsft_init_guru(&plan_adjoint,NQ[iNQ],m_total, NFSFT_NORMALIZED |
00934 ((use_nfft!=NO)?(0U):(NFSFT_USE_NDFT)) |
00935 ((use_fpt!=NO)?(0U):(NFSFT_USE_DPT)),
00936 ((NQ[iNQ]>512)?(0U):(PRE_PHI_HUT | PRE_PSI)) | FFTW_INIT |
00937 FFT_OUT_OF_PLACE, cutoff);
00938
00939 plan_adjoint_ptr = &plan_adjoint;
00940
00941 if (mode == RANDOM)
00942 {
00943 nfsft_init_guru(&plan,NQ[iNQ],m_compare, NFSFT_NORMALIZED |
00944 ((use_nfft!=NO)?(0U):(NFSFT_USE_NDFT)) |
00945 ((use_fpt!=NO)?(0U):(NFSFT_USE_DPT)),
00946 ((NQ[iNQ]>512)?(0U):(PRE_PHI_HUT | PRE_PSI)) | FFTW_INIT |
00947 FFT_OUT_OF_PLACE, cutoff);
00948 plan_ptr = &plan;
00949 }
00950 else
00951 {
00952 plan_ptr = &plan_adjoint;
00953 }
00954
00955 f_hat = (double _Complex*) nfft_malloc(NFSFT_F_HAT_SIZE(NQ[iNQ])*sizeof(double _Complex));
00956
00957 plan_adjoint_ptr->f_hat = f_hat;
00958 plan_adjoint_ptr->x = x_grid;
00959 plan_adjoint_ptr->f = f_grid;
00960
00961 plan_ptr->f_hat = f_hat;
00962 plan_ptr->x = x_compare;
00963 plan_ptr->f = f;
00964
00965
00966
00967 nfsft_precompute_x(plan_adjoint_ptr);
00968 if (plan_adjoint_ptr != plan_ptr)
00969 {
00970 nfsft_precompute_x(plan_ptr);
00971 }
00972
00973
00974 t_avg = 0.0;
00975 err_infty_avg = 0.0;
00976 err_2_avg = 0.0;
00977
00978
00979 for (i = 0; i < 1; i++)
00980 {
00981
00982
00983
00984
00985
00986
00987 t0 = getticks();
00988
00989
00990
00991
00992
00993 d = 0;
00994 for (k = 0; k < m_theta; k++)
00995 {
00996 for (n = 0; n < m_phi; n++)
00997 {
00998
00999
01000
01001 f_grid[d] *= w[k];
01002 d++;
01003 }
01004 }
01005
01006 t1 = getticks();
01007 t_avg += nfft_elapsed_seconds(t1,t0);
01008
01009 nfft_free(w);
01010
01011 t0 = getticks();
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024 if (use_nfsft != NO)
01025 {
01026
01027 nfsft_adjoint(plan_adjoint_ptr);
01028 }
01029 else
01030 {
01031
01032 nfsft_adjoint_direct(plan_adjoint_ptr);
01033 }
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048 if (use_nfsft != NO)
01049 {
01050
01051 nfsft_trafo(plan_ptr);
01052 }
01053 else
01054 {
01055
01056 nfsft_trafo_direct(plan_ptr);
01057 }
01058
01059 t1 = getticks();
01060 t_avg += nfft_elapsed_seconds(t1,t0);
01061
01062
01063
01064
01065 nfsft_finalize(plan_adjoint_ptr);
01066 if (plan_ptr != plan_adjoint_ptr)
01067 {
01068 nfsft_finalize(plan_ptr);
01069 }
01070
01071
01072 nfft_free(f_hat);
01073 nfft_free(x_grid);
01074
01075 err_infty_avg += X(error_l_infty_complex)(f, f_compare, m_compare);
01076 err_2_avg += X(error_l_2_complex)(f, f_compare, m_compare);
01077
01078 nfft_free(f_grid);
01079
01080 if (mode == RANDOM)
01081 {
01082 }
01083 else
01084 {
01085 nfft_free(f_compare);
01086 }
01087
01088
01089
01090
01091
01092
01093 }
01094
01095
01096
01097
01098 t_avg = t_avg/((double)repetitions);
01099
01100
01101 err_infty_avg = err_infty_avg/((double)repetitions);
01102
01103
01104 err_2_avg = err_2_avg/((double)repetitions);
01105
01106
01107 fprintf(stdout,"%+le %+le %+le\n", t_avg, err_infty_avg, err_2_avg);
01108 fprintf(stderr,"%d: %4d %4d %+le %+le %+le\n", tc, NQ[iNQ], SQ[iNQ],
01109 t_avg, err_infty_avg, err_2_avg);
01110 }
01111 }
01112
01113 fprintf(stderr,"\n");
01114
01115
01116 nfsft_forget();
01117
01118
01119 nfft_free(NQ);
01120 nfft_free(SQ);
01121 if (testmode == TIMING)
01122 {
01123 nfft_free(RQ);
01124 }
01125
01126 if (mode == RANDOM)
01127 {
01128 nfft_free(x_compare);
01129 nfft_free(f_compare);
01130 nfft_free(f);
01131 }
01132
01133 if (testmode == TIMING)
01134 {
01135
01136 nfft_free(f_hat);
01137 nfft_free(f);
01138 nfft_free(x_grid);
01139 }
01140
01141 }
01142
01143
01144 return EXIT_SUCCESS;
01145 }
01146