00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00029
00030 #include <stdio.h>
00031 #include <stdlib.h>
00032 #include <math.h>
00033 #include <complex.h>
00034
00035
00036 #include "nfft3util.h"
00037
00038 #include "nfft3.h"
00039
00040 #include "legendre.h"
00041
00043 enum boolean {NO = 0, YES = 1};
00044
00055 int main (int argc, char **argv)
00056 {
00057 int T;
00058 int N;
00059 int M;
00060 int M2;
00061
00062 int t;
00063 double time;
00064 double err_f;
00065 nfsft_plan plan;
00066 nfsft_plan plan2;
00067 solver_plan_complex iplan;
00068 int j;
00069 int k;
00070 int m;
00071 int n;
00072 int use_nfsft;
00073 int use_nfft;
00074 int use_fpt;
00075 int cutoff;
00076 double threshold;
00077 double re;
00078 double im;
00079 double a,b;
00080 double *scratch;
00081 double xs;
00082 double *ys;
00083 double *temp;
00084 double _Complex *temp2;
00085 int qlength;
00086 double *qweights;
00087 fftw_plan fplan;
00088 fpt_set set;
00089 int npt;
00090 int npt_exp;
00091 double *alpha, *beta, *gamma;
00092
00093
00094 fscanf(stdin,"testcases=%d\n",&T);
00095 fprintf(stderr,"%d\n",T);
00096
00097
00098 for (t = 0; t < T; t++)
00099 {
00100
00101 fscanf(stdin,"nfsft=%d\n",&use_nfsft);
00102 fprintf(stderr,"%d\n",use_nfsft);
00103 if (use_nfsft != NO)
00104 {
00105
00106 fscanf(stdin,"nfft=%d\n",&use_nfft);
00107 fprintf(stderr,"%d\n",use_nfsft);
00108 if (use_nfft != NO)
00109 {
00110
00111 fscanf(stdin,"cutoff=%d\n",&cutoff);
00112 fprintf(stderr,"%d\n",cutoff);
00113 }
00114 else
00115 {
00116
00117
00118 cutoff = 1;
00119 }
00120
00121 fscanf(stdin,"fpt=%d\n",&use_fpt);
00122 fprintf(stderr,"%d\n",use_fpt);
00123 if (use_fpt != NO)
00124 {
00125
00126 fscanf(stdin,"threshold=%lf\n",&threshold);
00127 fprintf(stderr,"%lf\n",threshold);
00128 }
00129 else
00130 {
00131
00132
00133 threshold = 1000.0;
00134 }
00135 }
00136 else
00137 {
00138
00139
00140 use_nfft = NO;
00141 use_fpt = NO;
00142 cutoff = 3;
00143 threshold = 1000.0;
00144 }
00145
00146
00147 fscanf(stdin,"bandwidth=%d\n",&N);
00148 fprintf(stderr,"%d\n",N);
00149
00150
00151 nfsft_precompute(N,threshold,
00152 ((use_nfsft==NO)?(NFSFT_NO_FAST_ALGORITHM):(0U)), 0U);
00153
00154
00155 fscanf(stdin,"nodes=%d\n",&M);
00156 fprintf(stderr,"%d\n",M);
00157
00158
00159 if ((N+1)*(N+1) > M)
00160 {
00161 nfft_next_power_of_2_exp(N, &npt, &npt_exp);
00162 fprintf(stderr, "npt = %d, npt_exp = %d\n", npt, npt_exp);
00163 fprintf(stderr,"Optimal interpolation!\n");
00164 scratch = (double*) nfft_malloc(4*sizeof(double));
00165 ys = (double*) nfft_malloc((N+1)*sizeof(double));
00166 temp = (double*) nfft_malloc((2*N+1)*sizeof(double));
00167 temp2 = (double _Complex*) nfft_malloc((N+1)*sizeof(double _Complex));
00168
00169 a = 0.0;
00170 for (j = 0; j <= N; j++)
00171 {
00172 xs = 2.0 + (2.0*j)/(N+1);
00173 ys[j] = (2.0-((j == 0)?(1.0):(0.0)))*4.0*nfft_bspline(4,xs,scratch);
00174
00175 a += ys[j];
00176 }
00177
00178 for (j = 0; j <= N; j++)
00179 {
00180 ys[j] *= 1.0/a;
00181 }
00182
00183 qlength = 2*N+1;
00184 qweights = (double*) nfft_malloc(qlength*sizeof(double));
00185
00186 fplan = fftw_plan_r2r_1d(N+1, qweights, qweights, FFTW_REDFT00, 0U);
00187 for (j = 0; j < N+1; j++)
00188 {
00189 qweights[j] = -2.0/(4*j*j-1);
00190 }
00191 fftw_execute(fplan);
00192 qweights[0] *= 0.5;
00193
00194 for (j = 0; j < N+1; j++)
00195 {
00196 qweights[j] *= 1.0/(2.0*N+1.0);
00197 qweights[2*N+1-1-j] = qweights[j];
00198 }
00199
00200 fplan = fftw_plan_r2r_1d(2*N+1, temp, temp, FFTW_REDFT00, 0U);
00201 for (j = 0; j <= N; j++)
00202 {
00203 temp[j] = ((j==0 || j == 2*N)?(1.0):(0.5))*ys[j];
00204 }
00205 for (j = N+1; j < 2*N+1; j++)
00206 {
00207 temp[j] = 0.0;
00208 }
00209 fftw_execute(fplan);
00210
00211 for (j = 0; j < 2*N+1; j++)
00212 {
00213 temp[j] *= qweights[j];
00214 }
00215
00216 fftw_execute(fplan);
00217
00218 for (j = 0; j < 2*N+1; j++)
00219 {
00220 temp[j] *= ((j==0 || j == 2*N)?(1.0):(0.5));
00221 if (j <= N)
00222 {
00223 temp2[j] = temp[j];
00224 }
00225 }
00226
00227 set = fpt_init(1, npt_exp, 0U);
00228
00229 alpha = (double*) nfft_malloc((N+2)*sizeof(double));
00230 beta = (double*) nfft_malloc((N+2)*sizeof(double));
00231 gamma = (double*) nfft_malloc((N+2)*sizeof(double));
00232
00233 alpha_al_row(alpha, N, 0);
00234 beta_al_row(beta, N, 0);
00235 gamma_al_row(gamma, N, 0);
00236
00237 fpt_precompute(set, 0, alpha, beta, gamma, 0, 1000.0);
00238
00239 fpt_transposed(set,0, temp2, temp2, N, 0U);
00240
00241 fpt_finalize(set);
00242
00243 nfft_free(alpha);
00244 nfft_free(beta);
00245 nfft_free(gamma);
00246
00247 fftw_destroy_plan(fplan);
00248
00249 nfft_free(scratch);
00250 nfft_free(qweights);
00251 nfft_free(ys);
00252 nfft_free(temp);
00253 }
00254
00255
00256 nfsft_init_guru(&plan, N, M,
00257 ((use_nfft!=0)?(0U):(NFSFT_USE_NDFT)) |
00258 ((use_fpt!=0)?(0U):(NFSFT_USE_DPT)) | NFSFT_MALLOC_F | NFSFT_MALLOC_X |
00259 NFSFT_MALLOC_F_HAT | NFSFT_NORMALIZED | NFSFT_ZERO_F_HAT,
00260 PRE_PHI_HUT | PRE_PSI | FFTW_INIT |
00261 FFT_OUT_OF_PLACE,
00262 cutoff);
00263
00264 if ((N+1)*(N+1) > M)
00265 {
00266 solver_init_advanced_complex(&iplan, (mv_plan_complex*)(&plan), CGNE | PRECOMPUTE_DAMP);
00267 }
00268 else
00269 {
00270 solver_init_advanced_complex(&iplan, (mv_plan_complex*)(&plan), CGNR | PRECOMPUTE_WEIGHT | PRECOMPUTE_DAMP);
00271 }
00272
00273
00274 for (j = 0; j < M; j++)
00275 {
00276 fscanf(stdin,"%le %le %le %le\n",&plan.x[2*j+1],&plan.x[2*j],&re,&im);
00277 plan.x[2*j+1] = plan.x[2*j+1]/(2.0*PI);
00278 plan.x[2*j] = plan.x[2*j]/(2.0*PI);
00279 if (plan.x[2*j] >= 0.5)
00280 {
00281 plan.x[2*j] = plan.x[2*j] - 1;
00282 }
00283 iplan.y[j] = re + _Complex_I * im;
00284 fprintf(stderr,"%le %le %le %le\n",plan.x[2*j+1],plan.x[2*j],
00285 creal(iplan.y[j]),cimag(iplan.y[j]));
00286 }
00287
00288
00289 fscanf(stdin,"nodes_eval=%d\n",&M2);
00290 fprintf(stderr,"%d\n",M2);
00291
00292
00293 nfsft_init_guru(&plan2, N, M2,
00294 ((use_nfft!=0)?(0U):(NFSFT_USE_NDFT)) |
00295 ((use_fpt!=0)?(0U):(NFSFT_USE_DPT)) | NFSFT_MALLOC_F | NFSFT_MALLOC_X |
00296 NFSFT_MALLOC_F_HAT | NFSFT_NORMALIZED | NFSFT_ZERO_F_HAT,
00297 PRE_PHI_HUT | PRE_PSI | FFTW_INIT |
00298 FFT_OUT_OF_PLACE,
00299 cutoff);
00300
00301
00302 for (j = 0; j < M2; j++)
00303 {
00304 fscanf(stdin,"%le %le\n",&plan2.x[2*j+1],&plan2.x[2*j]);
00305 plan2.x[2*j+1] = plan2.x[2*j+1]/(2.0*PI);
00306 plan2.x[2*j] = plan2.x[2*j]/(2.0*PI);
00307 if (plan2.x[2*j] >= 0.5)
00308 {
00309 plan2.x[2*j] = plan2.x[2*j] - 1;
00310 }
00311 fprintf(stderr,"%le %le\n",plan2.x[2*j+1],plan2.x[2*j]);
00312 }
00313
00314 nfsft_precompute_x(&plan);
00315
00316 nfsft_precompute_x(&plan2);
00317
00318
00319 if ((N+1)*(N+1) > M)
00320 {
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333 for (j = 0; j < plan.N_total; j++)
00334 {
00335 iplan.w_hat[j] = 0.0;
00336 }
00337
00338 for (k = 0; k <= N; k++)
00339 {
00340 for (j = -k; j <= k; j++)
00341 {
00342 iplan.w_hat[NFSFT_INDEX(k,j,&plan)] = 1.0/(pow(k+1.0,2.0)); ;
00343 }
00344 }
00345 }
00346 else
00347 {
00348 for (j = 0; j < plan.N_total; j++)
00349 {
00350 iplan.w_hat[j] = 0.0;
00351 }
00352
00353 for (k = 0; k <= N; k++)
00354 {
00355 for (j = -k; j <= k; j++)
00356 {
00357 iplan.w_hat[NFSFT_INDEX(k,j,&plan)] = 1/(pow(k+1.0,2.5));
00358 }
00359 }
00360
00361
00362 nfft_voronoi_weights_S2(iplan.w, plan.x, M);
00363
00364
00365 a = 0.0;
00366 for (j = 0; j < plan.M_total; j++)
00367 {
00368 fprintf(stderr,"%le\n",iplan.w[j]);
00369 a += iplan.w[j];
00370 }
00371 fprintf(stderr,"sum = %le\n",a);
00372 }
00373
00374 fprintf(stderr, "N_total = %d\n", plan.N_total);
00375 fprintf(stderr, "M_total = %d\n", plan.M_total);
00376
00377
00378 for (k = 0; k < plan.N_total; k++)
00379 {
00380 iplan.f_hat_iter[k] = 0.0;
00381 }
00382
00383
00384 solver_before_loop_complex(&iplan);
00385
00386
00387
00388
00389
00390
00391 for (m = 0; m < 29; m++)
00392 {
00393 fprintf(stderr,"Residual ||r||=%e,\n",sqrt(iplan.dot_r_iter));
00394 solver_loop_one_step_complex(&iplan);
00395 }
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413 NFFT_SWAP_complex(iplan.f_hat_iter, plan2.f_hat);
00414 nfsft_trafo(&plan2);
00415 NFFT_SWAP_complex(iplan.f_hat_iter, plan2.f_hat);
00416 for (k = 0; k < plan2.M_total; k++)
00417 {
00418 fprintf(stdout,"%le\n",cabs(plan2.f[k]));
00419 }
00420
00421 solver_finalize_complex(&iplan);
00422
00423 nfsft_finalize(&plan);
00424
00425 nfsft_finalize(&plan2);
00426
00427
00428 nfsft_forget();
00429
00430 if ((N+1)*(N+1) > M)
00431 {
00432 nfft_free(temp2);
00433 }
00434
00435 }
00436
00437
00438 return EXIT_SUCCESS;
00439 }
00440