00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00027
00028 #include <stdio.h>
00029 #include <stdlib.h>
00030 #include <math.h>
00031 #include <float.h>
00032 #include <complex.h>
00033
00034
00035 #include "nfft3.h"
00036
00037
00038 #include "nfft3util.h"
00039
00040
00041 #define SYMBOL_ABEL_POISSON(k,h) (pow(h,k))
00042
00043
00044 #define SYMBOL_SINGULARITY(k,h) ((2.0/(2*k+1))*pow(h,k))
00045
00046
00047
00049 #define KT_ABEL_POISSON (0)
00050
00051 #define KT_SINGULARITY (1)
00052
00053 #define KT_LOC_SUPP (2)
00054
00055 #define KT_GAUSSIAN (3)
00056
00058 enum pvalue {NO = 0, YES = 1, BOTH = 2};
00059
00074 static inline double innerProduct(const double phi1, const double theta1,
00075 const double phi2, const double theta2)
00076 {
00077 double pi2theta1 = PI2*theta1, pi2theta2 = PI2*theta2;
00078 return (cos(pi2theta1)*cos(pi2theta2)
00079 + sin(pi2theta1)*sin(pi2theta2)*cos(PI2*(phi1-phi2)));
00080 }
00081
00093 static inline double poissonKernel(const double x, const double h)
00094 {
00095 return (1.0/(PI4))*((1.0-h)*(1.0+h))/pow(sqrt(1.0-2.0*h*x+h*h),3.0);
00096 }
00097
00109 static inline double singularityKernel(const double x, const double h)
00110 {
00111 return (1.0/(PI2))/sqrt(1.0-2.0*h*x+h*h);
00112 }
00113
00127 static inline double locallySupportedKernel(const double x, const double h,
00128 const double lambda)
00129 {
00130 return (x<=h)?(0.0):(pow((x-h),lambda));
00131 }
00132
00145 static inline double gaussianKernel(const double x, const double sigma)
00146 {
00147 return exp(2.0*sigma*(x-1.0));
00148 }
00149
00160 int main (int argc, char **argv)
00161 {
00162 double **p;
00163
00164 int *m;
00165 int **ld;
00166
00167 int ip;
00168 int im;
00169 int ild;
00170 int ipp;
00171 int ip_max;
00172 int im_max;
00173 int ild_max;
00174 int ipp_max;
00175 int tc_max;
00176 int m_max;
00177
00178 int l_max;
00179
00180 int d_max;
00181
00182 long ld_max_prec;
00183
00184 long l_max_prec;
00185
00186 int tc;
00187 int kt;
00188 int cutoff;
00189 double threshold;
00190 double t_d;
00191 double t_dp;
00192
00193 double t_fd;
00194 double t_f;
00195 double temp;
00196 double err_f;
00197 double err_fd;
00198 double t;
00199 int precompute = NO;
00200 fftw_complex *ptr;
00201 double* steed;
00202 fftw_complex *b;
00203 fftw_complex *f_hat;
00204 fftw_complex *a;
00205 double *xi;
00206 double *eta;
00207 fftw_complex *f_m;
00208 fftw_complex *f;
00209 fftw_complex *prec = NULL;
00210 nfsft_plan plan;
00211 nfsft_plan plan_adjoint;
00212 int i;
00213 int k;
00214 int n;
00215 int d;
00216 int l;
00217 int use_nfsft;
00218 int use_nfft;
00219 int use_fpt;
00220 int rinc;
00221 double constant;
00222
00223
00224 fscanf(stdin,"testcases=%d\n",&tc_max);
00225 fprintf(stdout,"%d\n",tc_max);
00226
00227
00228 for (tc = 0; tc < tc_max; tc++)
00229 {
00230
00231 fscanf(stdin,"nfsft=%d\n",&use_nfsft);
00232 fprintf(stdout,"%d\n",use_nfsft);
00233 if (use_nfsft != NO)
00234 {
00235
00236 fscanf(stdin,"nfft=%d\n",&use_nfft);
00237 fprintf(stdout,"%d\n",use_nfft);
00238 if (use_nfft != NO)
00239 {
00240
00241 fscanf(stdin,"cutoff=%d\n",&cutoff);
00242 fprintf(stdout,"%d\n",cutoff);
00243 }
00244 else
00245 {
00246
00247
00248 cutoff = 1;
00249 }
00250
00251 fscanf(stdin,"fpt=%d\n",&use_fpt);
00252 fprintf(stdout,"%d\n",use_fpt);
00253
00254 fscanf(stdin,"threshold=%lf\n",&threshold);
00255 fprintf(stdout,"%lf\n",threshold);
00256 }
00257 else
00258 {
00259
00260
00261 cutoff = 3;
00262 threshold = 1000000000000.0;
00263 }
00264
00265
00266 m_max = 0;
00267
00268 l_max = 0;
00269
00270 d_max = 0;
00271
00272 l_max_prec = 0;
00273
00274 ld_max_prec = 0;
00275
00276
00277
00278 fscanf(stdin,"kernel=%d\n",&kt);
00279 fprintf(stdout,"%d\n",kt);
00280
00281
00282 fscanf(stdin,"parameter_sets=%d\n",&ip_max);
00283 fprintf(stdout,"%d\n",ip_max);
00284
00285
00286 p = (double**) nfft_malloc(ip_max*sizeof(double*));
00287
00288
00289
00290
00291 fscanf(stdin,"parameters=%d\n",&ipp_max);
00292 fprintf(stdout,"%d\n",ipp_max);
00293
00294 for (ip = 0; ip < ip_max; ip++)
00295 {
00296
00297 p[ip] = (double*) nfft_malloc(ipp_max*sizeof(double));
00298
00299
00300 for (ipp = 0; ipp < ipp_max; ipp++)
00301 {
00302
00303 fscanf(stdin,"%lf\n",&p[ip][ipp]);
00304 fprintf(stdout,"%lf\n",p[ip][ipp]);
00305 }
00306 }
00307
00308
00309 fscanf(stdin,"bandwidths=%d\n",&im_max);
00310 fprintf(stdout,"%d\n",im_max);
00311 m = (int*) nfft_malloc(im_max*sizeof(int));
00312
00313
00314 for (im = 0; im < im_max; im++)
00315 {
00316
00317 fscanf(stdin,"%d\n",&m[im]);
00318 fprintf(stdout,"%d\n",m[im]);
00319 m_max = NFFT_MAX(m_max,m[im]);
00320 }
00321
00322
00323 fscanf(stdin,"node_sets=%d\n",&ild_max);
00324 fprintf(stdout,"%d\n",ild_max);
00325 ld = (int**) nfft_malloc(ild_max*sizeof(int*));
00326
00327
00328 for (ild = 0; ild < ild_max; ild++)
00329 {
00330
00331 ld[ild] = (int*) nfft_malloc(5*sizeof(int));
00332
00333
00334 fscanf(stdin,"L=%d ",&ld[ild][0]);
00335 fprintf(stdout,"%d\n",ld[ild][0]);
00336 l_max = NFFT_MAX(l_max,ld[ild][0]);
00337
00338
00339 fscanf(stdin,"D=%d ",&ld[ild][1]);
00340 fprintf(stdout,"%d\n",ld[ild][1]);
00341 d_max = NFFT_MAX(d_max,ld[ild][1]);
00342
00343
00344 fscanf(stdin,"compare=%d ",&ld[ild][2]);
00345 fprintf(stdout,"%d\n",ld[ild][2]);
00346
00347
00348 if (ld[ild][2] == YES)
00349 {
00350
00351 fscanf(stdin,"precomputed=%d\n",&ld[ild][3]);
00352 fprintf(stdout,"%d\n",ld[ild][3]);
00353
00354
00355
00356 fscanf(stdin,"repetitions=%d\n",&ld[ild][4]);
00357 fprintf(stdout,"%d\n",ld[ild][4]);
00358
00359
00360 if (ld[ild][3] == YES)
00361 {
00362
00363 ld_max_prec = NFFT_MAX(ld_max_prec,ld[ild][0]*ld[ild][1]);
00364
00365 l_max_prec = NFFT_MAX(l_max_prec,ld[ild][0]);
00366
00367 precompute = YES;
00368 }
00369 }
00370 else
00371 {
00372
00373 ld[ild][4] = 1;
00374 }
00375 }
00376
00377
00378 b = (fftw_complex*) nfft_malloc(l_max*sizeof(fftw_complex));
00379 eta = (double*) nfft_malloc(2*l_max*sizeof(double));
00380 f_hat = (fftw_complex*) nfft_malloc(NFSFT_F_HAT_SIZE(m_max)*sizeof(fftw_complex));
00381 a = (fftw_complex*) nfft_malloc((m_max+1)*sizeof(fftw_complex));
00382 xi = (double*) nfft_malloc(2*d_max*sizeof(double));
00383 f_m = (fftw_complex*) nfft_malloc(d_max*sizeof(fftw_complex));
00384 f = (fftw_complex*) nfft_malloc(d_max*sizeof(fftw_complex));
00385
00386
00387 if (precompute == YES)
00388 {
00389 prec = (fftw_complex*) nfft_malloc(ld_max_prec*sizeof(fftw_complex));
00390 }
00391
00392
00393 for (l = 0; l < l_max; l++)
00394 {
00395 b[l] = (((double)rand())/RAND_MAX) - 0.5;
00396 eta[2*l] = (((double)rand())/RAND_MAX) - 0.5;
00397 eta[2*l+1] = acos(2.0*(((double)rand())/RAND_MAX) - 1.0)/(PI2);
00398 }
00399
00400
00401 for (d = 0; d < d_max; d++)
00402 {
00403 xi[2*d] = (((double)rand())/RAND_MAX) - 0.5;
00404 xi[2*d+1] = acos(2.0*(((double)rand())/RAND_MAX) - 1.0)/(PI2);
00405 }
00406
00407
00408 nfsft_precompute(m_max,threshold,
00409 ((use_nfsft==NO)?(NFSFT_NO_FAST_ALGORITHM):(0U)), 0U);
00410
00411
00412 for (ip = 0; ip < ip_max; ip++)
00413 {
00414
00415 switch (kt)
00416 {
00417 case KT_ABEL_POISSON:
00418
00419 for (k = 0; k <= m_max; k++)
00420 a[k] = SYMBOL_ABEL_POISSON(k,p[ip][0]);
00421 break;
00422
00423 case KT_SINGULARITY:
00424
00425
00426 for (k = 0; k <= m_max; k++)
00427 a[k] = SYMBOL_SINGULARITY(k,p[ip][0]);
00428 break;
00429
00430 case KT_LOC_SUPP:
00431
00432
00433 a[0] = 1.0;
00434 if (1 <= m_max)
00435 a[1] = ((p[ip][1]+1+p[ip][0])/(p[ip][1]+2.0))*a[0];
00436 for (k = 2; k <= m_max; k++)
00437 a[k] = (1.0/(k+p[ip][1]+1))*((2*k-1)*p[ip][0]*a[k-1] -
00438 (k-p[ip][1]-2)*a[k-2]);
00439 break;
00440
00441 case KT_GAUSSIAN:
00442
00443 steed = (double*) nfft_malloc((m_max+1)*sizeof(double));
00444 nfft_smbi(2.0*p[ip][0],0.5,m_max+1,2,steed);
00445 for (k = 0; k <= m_max; k++)
00446 a[k] = PI2*(sqrt(PI/p[ip][0]))*steed[k];
00447
00448 nfft_free(steed);
00449 break;
00450 }
00451
00452
00453 for (k = 0; k <= m_max; k++)
00454 a[k] *= (2*k+1)/(PI4);
00455
00456
00457 for (ild = 0; ild < ild_max; ild++)
00458 {
00459
00460 if (ld[ild][2] != NO)
00461 {
00462
00463
00464 if (ld[ild][3] != NO)
00465 {
00466
00467 ptr = prec;
00468
00469 rinc = l_max_prec-ld[ild][0];
00470
00471
00472 for (d = 0; d < ld[ild][1]; d++)
00473 {
00474
00475 for (l = 0; l < ld[ild][0]; l++)
00476 {
00477
00478
00479 temp = innerProduct(eta[2*l],eta[2*l+1],xi[2*d],xi[2*d+1]);
00480
00481
00482 switch (kt)
00483 {
00484 case KT_ABEL_POISSON:
00485
00486 *ptr++ = poissonKernel(temp,p[ip][0]);
00487 break;
00488
00489 case KT_SINGULARITY:
00490
00491
00492 *ptr++ = singularityKernel(temp,p[ip][0]);
00493 break;
00494
00495 case KT_LOC_SUPP:
00496
00497
00498 *ptr++ = locallySupportedKernel(temp,p[ip][0],p[ip][1]);
00499 break;
00500
00501 case KT_GAUSSIAN:
00502
00503
00504 *ptr++ = gaussianKernel(temp,p[ip][0]);
00505 break;
00506 }
00507 }
00508
00509 ptr += rinc;
00510 }
00511
00512
00513 t_dp = 0.0;
00514
00515
00516 t = nfft_second();
00517
00518
00519 for (i = 0; i < ld[ild][4]; i++)
00520 {
00521
00522
00523 ptr = prec;
00524
00525 rinc = l_max_prec-ld[ild][0];
00526
00527
00528 if (kt == KT_LOC_SUPP)
00529 {
00530
00531
00532
00533 constant = ((p[ip][1]+1)/(PI2*pow(1-p[ip][0],p[ip][1]+1)));
00534
00535
00536 for (d = 0; d < ld[ild][1]; d++)
00537 {
00538
00539 f[d] = 0.0;
00540
00541
00542 for (l = 0; l < ld[ild][0]; l++)
00543 f[d] += b[l]*(*ptr++);
00544
00545
00546 f[d] *= constant;
00547
00548
00549 ptr += rinc;
00550 }
00551 }
00552 else
00553 {
00554
00555 for (d = 0; d < ld[ild][1]; d++)
00556 {
00557
00558 f[d] = 0.0;
00559
00560
00561 for (l = 0; l < ld[ild][0]; l++)
00562 f[d] += b[l]*(*ptr++);
00563
00564
00565 ptr += rinc;
00566 }
00567 }
00568 }
00569
00570
00571 t_dp = nfft_second() - t;
00572
00573
00574 t_dp = t_dp/((double)ld[ild][4]);
00575 }
00576 else
00577 {
00578
00579 t_dp = -1.0;
00580 }
00581
00582
00583 t_d = 0.0;
00584
00585
00586 t = nfft_second();
00587
00588
00589 for (i = 0; i < ld[ild][4]; i++)
00590 {
00591
00592 switch (kt)
00593 {
00594 case KT_ABEL_POISSON:
00595
00596
00597 for (d = 0; d < ld[ild][1]; d++)
00598 {
00599
00600 f[d] = 0.0;
00601
00602
00603 for (l = 0; l < ld[ild][0]; l++)
00604 {
00605
00606
00607 temp = innerProduct(eta[2*l],eta[2*l+1],xi[2*d],xi[2*d+1]);
00608
00609
00610
00611 f[d] += b[l]*poissonKernel(temp,p[ip][0]);
00612 }
00613 }
00614 break;
00615
00616 case KT_SINGULARITY:
00617
00618 for (d = 0; d < ld[ild][1]; d++)
00619 {
00620
00621 f[d] = 0.0;
00622
00623
00624 for (l = 0; l < ld[ild][0]; l++)
00625 {
00626
00627
00628 temp = innerProduct(eta[2*l],eta[2*l+1],xi[2*d],xi[2*d+1]);
00629
00630
00631
00632 f[d] += b[l]*singularityKernel(temp,p[ip][0]);
00633 }
00634 }
00635 break;
00636
00637 case KT_LOC_SUPP:
00638
00639 constant = ((p[ip][1]+1)/(PI2*pow(1-p[ip][0],p[ip][1]+1)));
00640
00641
00642 for (d = 0; d < ld[ild][1]; d++)
00643 {
00644
00645 f[d] = 0.0;
00646
00647
00648 for (l = 0; l < ld[ild][0]; l++)
00649 {
00650
00651
00652 temp = innerProduct(eta[2*l],eta[2*l+1],xi[2*d],xi[2*d+1]);
00653
00654
00655
00656 f[d] += b[l]*locallySupportedKernel(temp,p[ip][0],p[ip][1]);
00657 }
00658
00659
00660 f[d] *= constant;
00661 }
00662 break;
00663
00664 case KT_GAUSSIAN:
00665
00666 for (d = 0; d < ld[ild][1]; d++)
00667 {
00668
00669 f[d] = 0.0;
00670
00671
00672 for (l = 0; l < ld[ild][0]; l++)
00673 {
00674
00675
00676 temp = innerProduct(eta[2*l],eta[2*l+1],xi[2*d],xi[2*d+1]);
00677
00678
00679 f[d] += b[l]*gaussianKernel(temp,p[ip][0]);
00680 }
00681 }
00682 break;
00683 }
00684 }
00685
00686
00687 t_d = nfft_second() - t;
00688
00689 t_d = t_d/((double)ld[ild][4]);
00690 }
00691 else
00692 {
00693
00694 t_d = -1.0;
00695 t_dp = -1.0;
00696 }
00697
00698
00699
00700 err_fd = -1.0;
00701 err_f = -1.0;
00702 t_fd = -1.0;
00703 t_f = -1.0;
00704
00705
00706 for (im = 0; im < im_max; im++)
00707 {
00708
00709 nfsft_init_guru(&plan_adjoint, m[im],ld[ild][0],
00710 ((use_nfft!=0)?(0U):(NFSFT_USE_NDFT)) |
00711 ((use_fpt!=0)?(0U):(NFSFT_USE_DPT)),
00712 PRE_PHI_HUT | PRE_PSI | FFTW_INIT |
00713 FFT_OUT_OF_PLACE, cutoff);
00714 nfsft_init_guru(&plan,m[im],ld[ild][1],
00715 ((use_nfft!=0)?(0U):(NFSFT_USE_NDFT)) |
00716 ((use_fpt!=0)?(0U):(NFSFT_USE_DPT)),
00717 PRE_PHI_HUT | PRE_PSI | FFTW_INIT |
00718 FFT_OUT_OF_PLACE,
00719 cutoff);
00720 plan_adjoint.f_hat = f_hat;
00721 plan_adjoint.x = eta;
00722 plan_adjoint.f = b;
00723 plan.f_hat = f_hat;
00724 plan.x = xi;
00725 plan.f = f_m;
00726 nfsft_precompute_x(&plan_adjoint);
00727 nfsft_precompute_x(&plan);
00728
00729
00730 if (use_nfsft == BOTH)
00731 {
00732
00733 t_fd = 0.0;
00734
00735
00736 t = nfft_second();
00737
00738
00739 for (i = 0; i < ld[ild][4]; i++)
00740 {
00741
00742
00743 ndsft_adjoint(&plan_adjoint);
00744
00745
00746 for (k = 0; k <= m[im]; k++)
00747 for (n = -k; n <= k; n++)
00748 f_hat[NFSFT_INDEX(k,n,&plan_adjoint)] *= a[k];
00749
00750
00751 ndsft_trafo(&plan);
00752
00753 }
00754
00755
00756 t_fd = nfft_second() - t;
00757
00758
00759 t_fd = t_fd/((double)ld[ild][4]);
00760
00761
00762 if (ld[ild][2] != NO)
00763 {
00764
00765 err_fd = nfft_error_l_infty_1_complex(f, f_m, ld[ild][1], b,
00766 ld[ild][0]);
00767 }
00768 }
00769
00770
00771 if (use_nfsft != NO)
00772 {
00773
00774 t_f = 0.0;
00775 }
00776 else
00777 {
00778
00779
00780 t_fd = 0.0;
00781 }
00782
00783
00784 t = nfft_second();
00785
00786
00787 for (i = 0; i < ld[ild][4]; i++)
00788 {
00789
00790 if (use_nfsft != NO)
00791 {
00792
00793 nfsft_adjoint(&plan_adjoint);
00794 }
00795 else
00796 {
00797
00798 ndsft_adjoint(&plan_adjoint);
00799 }
00800
00801
00802 for (k = 0; k <= m[im]; k++)
00803 for (n = -k; n <= k; n++)
00804 f_hat[NFSFT_INDEX(k,n,&plan_adjoint)] *= a[k];
00805
00806
00807 if (use_nfsft != NO)
00808 {
00809
00810 nfsft_trafo(&plan);
00811 }
00812 else
00813 {
00814
00815 ndsft_trafo(&plan);
00816 }
00817 }
00818
00819
00820 if (use_nfsft != NO)
00821 {
00822
00823 t_f = nfft_second() - t;
00824 }
00825 else
00826 {
00827
00828 t_fd = nfft_second() - t;
00829 }
00830
00831
00832 if (use_nfsft != NO)
00833 {
00834
00835 t_f = t_f/((double)ld[ild][4]);
00836 }
00837 else
00838 {
00839
00840 t_fd = t_fd/((double)ld[ild][4]);
00841 }
00842
00843
00844 if (ld[ild][2] != NO)
00845 {
00846
00847 if (use_nfsft != NO)
00848 {
00849
00850 err_f = nfft_error_l_infty_1_complex(f, f_m, ld[ild][1], b,
00851 ld[ild][0]);
00852 }
00853 else
00854 {
00855
00856 err_fd = nfft_error_l_infty_1_complex(f, f_m, ld[ild][1], b,
00857 ld[ild][0]);
00858 }
00859 }
00860
00861
00862 fprintf(stdout,"%e\n%e\n%e\n%e\n%e\n%e\n\n",t_d,t_dp,t_fd,t_f,err_fd,
00863 err_f);
00864
00865
00866 nfsft_finalize(&plan_adjoint);
00867 nfsft_finalize(&plan);
00868 }
00869
00870 }
00871 }
00872
00873
00874 nfsft_forget();
00875
00876
00877
00878 if (precompute == YES)
00879 {
00880
00881 nfft_free(prec);
00882 }
00883
00884 nfft_free(f);
00885 nfft_free(f_m);
00886 nfft_free(xi);
00887 nfft_free(eta);
00888 nfft_free(a);
00889 nfft_free(f_hat);
00890 nfft_free(b);
00891
00892
00893 for (ild = 0; ild < ild_max; ild++)
00894 nfft_free(ld[ild]);
00895 nfft_free(ld);
00896
00897
00898 nfft_free(m);
00899
00900
00901 for (ip = 0; ip < ip_max; ip++)
00902 nfft_free(p[ip]);
00903 nfft_free(p);
00904 }
00905
00906
00907 return EXIT_SUCCESS;
00908 }
00909