102 trmesh_(&Mlocal, x, y, z, list, lptr, lend, &lnew, near, next, dist, &ier);
108 crlist_(&Mlocal, &Mlocal, x, y, z, list, lend, lptr, &lnew, ltri, listc, &nb, xc,
181int main (
int argc,
char **argv)
206 double _Complex *temp2;
213 double *alpha, *beta, *gamma;
216 fscanf(stdin,
"testcases=%d\n",&T);
217 fprintf(stderr,
"%d\n",T);
220 for (t = 0; t < T; t++)
223 fscanf(stdin,
"nfsft=%d\n",&use_nfsft);
224 fprintf(stderr,
"%d\n",use_nfsft);
228 fscanf(stdin,
"nfft=%d\n",&use_nfft);
229 fprintf(stderr,
"%d\n",use_nfsft);
233 fscanf(stdin,
"cutoff=%d\n",&cutoff);
234 fprintf(stderr,
"%d\n",cutoff);
243 fscanf(stdin,
"fpt=%d\n",&use_fpt);
244 fprintf(stderr,
"%d\n",use_fpt);
248 fscanf(stdin,
"threshold=%lf\n",&threshold);
249 fprintf(stderr,
"%lf\n",threshold);
269 fscanf(stdin,
"bandwidth=%d\n",&N);
270 fprintf(stderr,
"%d\n",N);
277 fscanf(stdin,
"nodes=%d\n",&M);
278 fprintf(stderr,
"%d\n",M);
283 X(next_power_of_2_exp_int)(N, &npt, &npt_exp);
284 fprintf(stderr,
"npt = %d, npt_exp = %d\n", npt, npt_exp);
285 fprintf(stderr,
"Optimal interpolation!\n");
287 temp = (
double*)
nfft_malloc((2*N+1)*
sizeof(double));
288 temp2 = (
double _Complex*)
nfft_malloc((N+1)*
sizeof(
double _Complex));
291 for (j = 0; j <= N; j++)
293 xs = 2.0 + (2.0*j)/(N+1);
294 ys[j] = (2.0-((j == 0)?(1.0):(0.0)))*4.0*nfft_bsplines(4,xs);
299 for (j = 0; j <= N; j++)
305 qweights = (
double*)
nfft_malloc(qlength*
sizeof(
double));
307 fplan = fftw_plan_r2r_1d(N+1, qweights, qweights, FFTW_REDFT00, 0U);
308 for (j = 0; j < N+1; j++)
310 qweights[j] = -2.0/(4*j*j-1);
315 for (j = 0; j < N+1; j++)
317 qweights[j] *= 1.0/(2.0*N+1.0);
318 qweights[2*N+1-1-j] = qweights[j];
321 fplan = fftw_plan_r2r_1d(2*N+1, temp, temp, FFTW_REDFT00, 0U);
322 for (j = 0; j <= N; j++)
324 temp[j] = ((j==0 || j == 2*N)?(1.0):(0.5))*ys[j];
326 for (j = N+1; j < 2*N+1; j++)
332 for (j = 0; j < 2*N+1; j++)
334 temp[j] *= qweights[j];
339 for (j = 0; j < 2*N+1; j++)
341 temp[j] *= ((j==0 || j == 2*N)?(1.0):(0.5));
348 set = fpt_init(1, npt_exp, 0U);
350 alpha = (
double*)
nfft_malloc((N+2)*
sizeof(double));
351 beta = (
double*)
nfft_malloc((N+2)*
sizeof(double));
352 gamma = (
double*)
nfft_malloc((N+2)*
sizeof(double));
354 alpha_al_row(alpha, N, 0);
355 beta_al_row(beta, N, 0);
356 gamma_al_row(gamma, N, 0);
358 fpt_precompute(set, 0, alpha, beta, gamma, 0, 1000.0);
360 fpt_transposed(set,0, temp2, temp2, N, 0U);
368 fftw_destroy_plan(fplan);
376 nfsft_init_guru(&plan, N, M,
394 for (j = 0; j < M; j++)
396 fscanf(stdin,
"%le %le %le %le\n",&plan.
x[2*j+1],&plan.
x[2*j],&re,&im);
397 plan.
x[2*j+1] = plan.
x[2*j+1]/(2.0*KPI);
398 plan.
x[2*j] = plan.
x[2*j]/(2.0*KPI);
399 if (plan.
x[2*j] >= 0.5)
401 plan.
x[2*j] = plan.
x[2*j] - 1;
403 iplan.
y[j] = re + _Complex_I * im;
404 fprintf(stderr,
"%le %le %le %le\n",plan.
x[2*j+1],plan.
x[2*j],
405 creal(iplan.
y[j]),cimag(iplan.
y[j]));
409 fscanf(stdin,
"nodes_eval=%d\n",&M2);
410 fprintf(stderr,
"%d\n",M2);
413 nfsft_init_guru(&plan2, N, M2,
422 for (j = 0; j < M2; j++)
424 fscanf(stdin,
"%le %le\n",&plan2.
x[2*j+1],&plan2.
x[2*j]);
425 plan2.
x[2*j+1] = plan2.
x[2*j+1]/(2.0*KPI);
426 plan2.
x[2*j] = plan2.
x[2*j]/(2.0*KPI);
427 if (plan2.
x[2*j] >= 0.5)
429 plan2.
x[2*j] = plan2.
x[2*j] - 1;
431 fprintf(stderr,
"%le %le\n",plan2.
x[2*j+1],plan2.
x[2*j]);
434 nfsft_precompute_x(&plan);
436 nfsft_precompute_x(&plan2);
453 for (j = 0; j < plan.
N_total; j++)
455 iplan.
w_hat[j] = 0.0;
458 for (k = 0; k <= N; k++)
460 for (j = -k; j <= k; j++)
468 for (j = 0; j < plan.
N_total; j++)
470 iplan.
w_hat[j] = 0.0;
473 for (k = 0; k <= N; k++)
475 for (j = -k; j <= k; j++)
482 voronoi_weights_S2(iplan.
w, plan.
x, M);
486 for (j = 0; j < plan.
M_total; j++)
488 fprintf(stderr,
"%le\n",iplan.
w[j]);
491 fprintf(stderr,
"sum = %le\n",a);
494 fprintf(stderr,
"N_total = %d\n", plan.
N_total);
495 fprintf(stderr,
"M_total = %d\n", plan.
M_total);
498 for (k = 0; k < plan.
N_total; k++)
504 solver_before_loop_complex(&iplan);
511 for (m = 0; m < 29; m++)
513 fprintf(stderr,
"Residual ||r||=%e,\n",sqrt(iplan.
dot_r_iter));
514 solver_loop_one_step_complex(&iplan);
536 for (k = 0; k < plan2.
M_total; k++)
538 fprintf(stdout,
"%le\n",cabs(plan2.
f[k]));
541 solver_finalize_complex(&iplan);