45 static void voronoi_weights_S2(R *w, R *xi, INT M)
76 x = (R*)
X(malloc)(M *
sizeof(R));
77 y = (R*)
X(malloc)(M *
sizeof(R));
78 z = (R*)
X(malloc)(M *
sizeof(R));
80 list = (
int*)
X(malloc)((6*M-12+1)*
sizeof(
int));
81 lptr = (
int*)
X(malloc)((6*M-12+1)*
sizeof(
int));
82 lend = (
int*)
X(malloc)((M+1)*
sizeof(
int));
83 near = (
int*)
X(malloc)((M+1)*
sizeof(
int));
84 next = (
int*)
X(malloc)((M+1)*
sizeof(
int));
85 dist = (R*)
X(malloc)((M+1)*
sizeof(R));
86 ltri = (
int*)
X(malloc)((6*M+1)*
sizeof(
int));
87 listc = (
int*)
X(malloc)((6*M-12+1)*
sizeof(
int));
88 xc = (R*)
X(malloc)((2*M-4+1)*
sizeof(R));
89 yc = (R*)
X(malloc)((2*M-4+1)*
sizeof(R));
90 zc = (R*)
X(malloc)((2*M-4+1)*
sizeof(R));
91 rc = (R*)
X(malloc)((2*M-4+1)*
sizeof(R));
92 vr = (R*)
X(malloc)(3*(2*M-4+1)*
sizeof(R));
96 for (k = 0; k < M; k++)
98 x[k] = SIN(K2PI*xi[2*k+1])*COS(K2PI*xi[2*k]);
99 y[k] = SIN(K2PI*xi[2*k+1])*SIN(K2PI*xi[2*k]);
100 z[k] = COS(K2PI*xi[2*k+1]);
104 trmesh_(&Mlocal, x, y, z, list, lptr, lend, &lnew, near, next, dist, &ier);
110 crlist_(&Mlocal, &Mlocal, x, y, z, list, lend, lptr, &lnew, ltri, listc, &nb, xc,
116 for (k = 0; k < M; k++)
134 vr[3*j+1] = yc[kv-1];
135 vr[3*j+2] = zc[kv-1];
140 for (el = 0; el < j; el++)
142 a += areas_(vr, &vr[3*(el+1)],&vr[3*(((el+1)%j)+1)]);
171 enum boolean {NO = 0, YES = 1};
183 int main (
int argc,
char **argv)
209 double _Complex *temp2;
216 double *alpha, *beta, *gamma;
219 fscanf(stdin,
"testcases=%d\n",&T);
220 fprintf(stderr,
"%d\n",T);
223 for (t = 0; t < T; t++)
226 fscanf(stdin,
"nfsft=%d\n",&use_nfsft);
227 fprintf(stderr,
"%d\n",use_nfsft);
231 fscanf(stdin,
"nfft=%d\n",&use_nfft);
232 fprintf(stderr,
"%d\n",use_nfsft);
236 fscanf(stdin,
"cutoff=%d\n",&cutoff);
237 fprintf(stderr,
"%d\n",cutoff);
246 fscanf(stdin,
"fpt=%d\n",&use_fpt);
247 fprintf(stderr,
"%d\n",use_fpt);
251 fscanf(stdin,
"threshold=%lf\n",&threshold);
252 fprintf(stderr,
"%lf\n",threshold);
272 fscanf(stdin,
"bandwidth=%d\n",&N);
273 fprintf(stderr,
"%d\n",N);
276 nfsft_precompute(N,threshold,
277 ((use_nfsft==NO)?(NFSFT_NO_FAST_ALGORITHM):(0U)), 0U);
280 fscanf(stdin,
"nodes=%d\n",&M);
281 fprintf(stderr,
"%d\n",M);
286 X(next_power_of_2_exp)(N, &npt, &npt_exp);
287 fprintf(stderr,
"npt = %d, npt_exp = %d\n", npt, npt_exp);
288 fprintf(stderr,
"Optimal interpolation!\n");
291 temp = (
double*)
nfft_malloc((2*N+1)*
sizeof(double));
292 temp2 = (
double _Complex*)
nfft_malloc((N+1)*
sizeof(
double _Complex));
295 for (j = 0; j <= N; j++)
297 xs = 2.0 + (2.0*j)/(N+1);
298 ys[j] = (2.0-((j == 0)?(1.0):(0.0)))*4.0*nfft_bspline(4,xs,scratch);
303 for (j = 0; j <= N; j++)
309 qweights = (
double*)
nfft_malloc(qlength*
sizeof(
double));
311 fplan = fftw_plan_r2r_1d(N+1, qweights, qweights, FFTW_REDFT00, 0U);
312 for (j = 0; j < N+1; j++)
314 qweights[j] = -2.0/(4*j*j-1);
319 for (j = 0; j < N+1; j++)
321 qweights[j] *= 1.0/(2.0*N+1.0);
322 qweights[2*N+1-1-j] = qweights[j];
325 fplan = fftw_plan_r2r_1d(2*N+1, temp, temp, FFTW_REDFT00, 0U);
326 for (j = 0; j <= N; j++)
328 temp[j] = ((j==0 || j == 2*N)?(1.0):(0.5))*ys[j];
330 for (j = N+1; j < 2*N+1; j++)
336 for (j = 0; j < 2*N+1; j++)
338 temp[j] *= qweights[j];
343 for (j = 0; j < 2*N+1; j++)
345 temp[j] *= ((j==0 || j == 2*N)?(1.0):(0.5));
352 set = fpt_init(1, npt_exp, 0U);
354 alpha = (
double*)
nfft_malloc((N+2)*
sizeof(double));
355 beta = (
double*)
nfft_malloc((N+2)*
sizeof(double));
356 gamma = (
double*)
nfft_malloc((N+2)*
sizeof(double));
358 alpha_al_row(alpha, N, 0);
359 beta_al_row(beta, N, 0);
360 gamma_al_row(gamma, N, 0);
362 fpt_precompute(set, 0, alpha, beta, gamma, 0, 1000.0);
364 fpt_transposed(set,0, temp2, temp2, N, 0U);
372 fftw_destroy_plan(fplan);
381 nfsft_init_guru(&plan, N, M,
382 ((use_nfft!=0)?(0U):(NFSFT_USE_NDFT)) |
383 ((use_fpt!=0)?(0U):(NFSFT_USE_DPT)) | NFSFT_MALLOC_F | NFSFT_MALLOC_X |
384 NFSFT_MALLOC_F_HAT | NFSFT_NORMALIZED | NFSFT_ZERO_F_HAT,
385 PRE_PHI_HUT | PRE_PSI | FFTW_INIT |
395 solver_init_advanced_complex(&iplan, (
nfft_mv_plan_complex*)(&plan), CGNR | PRECOMPUTE_WEIGHT | PRECOMPUTE_DAMP);
399 for (j = 0; j < M; j++)
401 fscanf(stdin,
"%le %le %le %le\n",&plan.
x[2*j+1],&plan.
x[2*j],&re,&im);
402 plan.
x[2*j+1] = plan.
x[2*j+1]/(2.0*KPI);
403 plan.
x[2*j] = plan.
x[2*j]/(2.0*KPI);
404 if (plan.
x[2*j] >= 0.5)
406 plan.
x[2*j] = plan.
x[2*j] - 1;
408 iplan.
y[j] = re + _Complex_I * im;
409 fprintf(stderr,
"%le %le %le %le\n",plan.
x[2*j+1],plan.
x[2*j],
410 creal(iplan.
y[j]),cimag(iplan.
y[j]));
414 fscanf(stdin,
"nodes_eval=%d\n",&M2);
415 fprintf(stderr,
"%d\n",M2);
418 nfsft_init_guru(&plan2, N, M2,
419 ((use_nfft!=0)?(0U):(NFSFT_USE_NDFT)) |
420 ((use_fpt!=0)?(0U):(NFSFT_USE_DPT)) | NFSFT_MALLOC_F | NFSFT_MALLOC_X |
421 NFSFT_MALLOC_F_HAT | NFSFT_NORMALIZED | NFSFT_ZERO_F_HAT,
422 PRE_PHI_HUT | PRE_PSI | FFTW_INIT |
427 for (j = 0; j < M2; j++)
429 fscanf(stdin,
"%le %le\n",&plan2.
x[2*j+1],&plan2.
x[2*j]);
430 plan2.
x[2*j+1] = plan2.
x[2*j+1]/(2.0*KPI);
431 plan2.
x[2*j] = plan2.
x[2*j]/(2.0*KPI);
432 if (plan2.
x[2*j] >= 0.5)
434 plan2.
x[2*j] = plan2.
x[2*j] - 1;
436 fprintf(stderr,
"%le %le\n",plan2.
x[2*j+1],plan2.
x[2*j]);
439 nfsft_precompute_x(&plan);
441 nfsft_precompute_x(&plan2);
458 for (j = 0; j < plan.
N_total; j++)
460 iplan.
w_hat[j] = 0.0;
463 for (k = 0; k <= N; k++)
465 for (j = -k; j <= k; j++)
467 iplan.
w_hat[NFSFT_INDEX(k,j,&plan)] = 1.0/(pow(k+1.0,2.0)); ;
473 for (j = 0; j < plan.
N_total; j++)
475 iplan.
w_hat[j] = 0.0;
478 for (k = 0; k <= N; k++)
480 for (j = -k; j <= k; j++)
482 iplan.
w_hat[NFSFT_INDEX(k,j,&plan)] = 1/(pow(k+1.0,2.5));
487 voronoi_weights_S2(iplan.
w, plan.
x, M);
491 for (j = 0; j < plan.
M_total; j++)
493 fprintf(stderr,
"%le\n",iplan.
w[j]);
496 fprintf(stderr,
"sum = %le\n",a);
499 fprintf(stderr,
"N_total = %d\n", plan.
N_total);
500 fprintf(stderr,
"M_total = %d\n", plan.
M_total);
503 for (k = 0; k < plan.
N_total; k++)
509 solver_before_loop_complex(&iplan);
516 for (m = 0; m < 29; m++)
518 fprintf(stderr,
"Residual ||r||=%e,\n",sqrt(iplan.
dot_r_iter));
519 solver_loop_one_step_complex(&iplan);
541 for (k = 0; k < plan2.
M_total; k++)
543 fprintf(stdout,
"%le\n",cabs(plan2.
f[k]));
546 solver_finalize_complex(&iplan);
548 nfsft_finalize(&plan);
550 nfsft_finalize(&plan2);
double * x
the nodes for ,
int main(int argc, char **argv)
The main program.
fftw_complex * f_hat
Fourier coefficients.
double * w
weighting factors
data structure for an NFSFT (nonequispaced fast spherical Fourier transform) plan with double precisi...
double dot_r_iter
weighted dotproduct of r_iter
Holds data for a set of cascade summations.
NFFT_INT M_total
Total number of samples.
#define X(name)
Include header for C99 complex datatype.
void * nfft_malloc(size_t n)
NFFT_INT N_total
Total number of Fourier coefficients.
fftw_complex * y
right hand side, samples
data structure for an inverse NFFT plan with double precision
#define CSWAP(x, y)
Swap two vectors.
double * w_hat
damping factors
fftw_complex * f_hat_iter
iterative solution