Aquila-Consortium/libsharp2
0
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

separate out pocketfft

Browse source
This commit is contained in:
Martin Reinecke 2019-02-27 11:12:44 +01:00
parent 540e7e44f8
commit bf43082182
7 changed files with 22 additions and 2272 deletions
Download patch file Download diff file Expand all files Collapse all files

14
libsharp/sharp.c
View file

@ -77,7 +77,7 @@ typedef struct
double phi0_;
dcmplx *shiftarr;
int s_shift;
rfft_plan plan;
pocketfft_plan_r plan;
int length;
int norot;
} ringhelper;
@ -90,7 +90,7 @@ static void ringhelper_init (ringhelper *self)
static void ringhelper_destroy (ringhelper *self)
{
if (self->plan) destroy_rfft_plan(self->plan);
if (self->plan) pocketfft_delete_plan_r(self->plan);
DEALLOC(self->shiftarr);
ringhelper_init(self);
}
@ -110,11 +110,11 @@ NOINLINE static void ringhelper_update (ringhelper *self, int nph, int mmax, dou
// double *tmp=(double *) self->shiftarr;
// sincos_multi (mmax+1, phi0, &tmp[1], &tmp[0], 2);
}
// if (!self->plan) self->plan=make_rfft_plan(nph);
// if (!self->plan) self->plan=pocketfft_make_plan_r(nph);
if (nph!=(int)self->length)
{
if (self->plan) destroy_rfft_plan(self->plan);
self->plan=make_rfft_plan(nph);
if (self->plan) pocketfft_delete_plan_r(self->plan);
self->plan=pocketfft_make_plan_r(nph);
self->length=nph;
}
}
@ -331,7 +331,7 @@ NOINLINE static void ringhelper_phase2ring (ringhelper *self,
}
}
data[1]=data[0];
rfft_backward (self->plan, &(data[1]), 1.);
pocketfft_backward_r (self->plan, &(data[1]), 1.);
}
NOINLINE static void ringhelper_ring2phase (ringhelper *self,
@ -350,7 +350,7 @@ NOINLINE static void ringhelper_ring2phase (ringhelper *self,
if (flags&SHARP_REAL_HARMONICS)
wgt *= sqrt_two;
rfft_forward (self->plan, &(data[1]), 1.);
pocketfft_forward_r (self->plan, &(data[1]), 1.);
data[0]=data[1];
data[1]=data[nph+1]=0.;

18
libsharp/sharp_geomhelpers.c
View file

@ -155,9 +155,9 @@ void sharp_make_fejer1_geom_info (int nrings, int ppring, double phi0,
weight[2*k ]=2./(1.-4.*k*k)*sin((k*pi)/nrings);
}
if ((nrings&1)==0) weight[nrings-1]=0.;
rfft_plan plan = make_rfft_plan(nrings);
rfft_backward(plan,weight,1.);
destroy_rfft_plan(plan);
pocketfft_plan_r plan = pocketfft_make_plan_r(nrings);
pocketfft_backward_r(plan,weight,1.);
pocketfft_delete_plan_r(plan);
for (int m=0; m<(nrings+1)/2; ++m)
{
@ -202,9 +202,9 @@ void sharp_make_cc_geom_info (int nrings, int ppring, double phi0,
for (int k=1; k<=(n/2-1); ++k)
weight[2*k-1]=2./(1.-4.*k*k) + dw;
weight[2*(n/2)-1]=(n-3.)/(2*(n/2)-1) -1. -dw*((2-(n&1))*n-1);
rfft_plan plan = make_rfft_plan(n);
rfft_backward(plan,weight,1.);
destroy_rfft_plan(plan);
pocketfft_plan_r plan = pocketfft_make_plan_r(n);
pocketfft_backward_r(plan,weight,1.);
pocketfft_delete_plan_r(plan);
weight[n]=weight[0];
for (int m=0; m<(nrings+1)/2; ++m)
@ -250,9 +250,9 @@ void sharp_make_fejer2_geom_info (int nrings, int ppring, double phi0,
for (int k=1; k<=(n/2-1); ++k)
weight[2*k-1]=2./(1.-4.*k*k);
weight[2*(n/2)-1]=(n-3.)/(2*(n/2)-1) -1.;
rfft_plan plan = make_rfft_plan(n);
rfft_backward(plan,weight,1.);
destroy_rfft_plan(plan);
pocketfft_plan_r plan = pocketfft_make_plan_r(n);
pocketfft_backward_r(plan,weight,1.);
pocketfft_delete_plan_r(plan);
for (int m=0; m<nrings; ++m)
weight[m]=weight[m+1];