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

first try; broken

Browse source
This commit is contained in:
Martin Reinecke 2018-12-18 13:45:25 +01:00
parent d257e92e54
commit 202507ecaa
3 changed files with 78 additions and 4 deletions
Download patch file Download diff file Expand all files Collapse all files

54
libsharp/sharp_core.c
View file

@ -280,6 +280,58 @@ NOINLINE static void calc_alm2map (sharp_job * restrict job,
alm2map_kernel(d, rf, alm, l, lmax, nv2);
}
NOINLINE static void calc_alm2map_alt (sharp_job * restrict job,
const sharp_Ylmgen_C * restrict gen, s0data_v * restrict d, int nth)
{
int l=gen->m, lmax=gen->lmax;
int nv2 = (nth+VLEN-1)/VLEN;
job->opcnt += (lmax+1-l) * 6*nth;
const double inv_sqrt4pi = 0.2820947917738781434740397257803862929220;
Tv mfac = vload(gen->mfac[gen->m]);
for (int i=0; i<nv2; ++i)
{
d->lam1[i]=vzero;
mypow(d->sth[i],l,gen->powlimit,&d->lam2[i],&d->scale[i]);
d->lam2[i] *= mfac;
Tvnormalize(&d->lam2[i],&d->scale[i],sharp_ftol);
}
const dcmplx * restrict alm=job->almtmp;
dcmplx * restrict alm2=RALLOC(dcmplx, gen->lmax+5);
{
for (int il=0, l=gen->m; l<=gen->lmax; ++il,l+=2)
{
dcmplx al = alm[l];
dcmplx al1 = (l+1>gen->lmax) ? 0. : alm[l+1];
dcmplx al2 = (l+2>gen->lmax) ? 0. : alm[l+2];
alm2[l ] = gen->alpha[il]*(gen->eps[l+1]*al + gen->eps[l+2]*al2);
alm2[l+1] = gen->alpha[il]*al1;
}
}
for (int i=0; i<nv2; ++i)
getCorfac(d->scale[i], &d->corfac[i], gen->cf);
for (int il=0, l=gen->m; l<=lmax; ++il, l+=2)
{
Tv ar1=vload(creal(alm2[l ])), ai1=vload(cimag(alm2[l ]));
Tv ar2=vload(creal(alm2[l+1])), ai2=vload(cimag(alm2[l+1]));
for (int i=0; i<nv2; ++i)
{
d->p1r[i] += d->lam2[i]*d->corfac[i]*ar1;
d->p1i[i] += d->lam2[i]*d->corfac[i]*ai1;
d->p2r[i] += d->cth[i]*d->lam2[i]*d->corfac[i]*ar2;
d->p2i[i] += d->cth[i]*d->lam2[i]*d->corfac[i]*ai2;
Tv tmp = (gen->a[il]*d->cth[i]*d->cth[i] + gen->b[il])*d->lam2[i] + d->lam1[i];
d->lam1[i] = d->lam2[i];
d->lam2[i] = tmp;
if (rescale(&d->lam1[i], &d->lam2[i], &d->scale[i], vload(sharp_ftol)))
getCorfac(d->scale[i], &d->corfac[i], gen->cf);
}
}
DEALLOC(alm2);
}
NOINLINE static void map2alm_kernel(s0data_v * restrict d,
const sharp_ylmgen_dbl2 * restrict rf, dcmplx * restrict alm, int l,
int lmax, int nv2)
@ -812,7 +864,7 @@ NOINLINE static void inner_loop_a2m(sharp_job *job, const int *ispair,
d.s.sth[i]=d.s.sth[nth-1];
d.s.p1r[i]=d.s.p1i[i]=d.s.p2r[i]=d.s.p2i[i]=0.;
}
calc_alm2map (job, gen, &d.v, nth);
calc_alm2map_alt (job, gen, &d.v, nth);
for (int i=0; i<nth; ++i)
{
int tgt=itgt[i];

26
libsharp/sharp_ylmgen_c.c
View file

@ -74,13 +74,17 @@ void sharp_Ylmgen_init (sharp_Ylmgen_C *gen, int l_max, int m_max, int spin)
gen->mfac[0] = inv_sqrt4pi;
for (int m=1; m<=gen->mmax; ++m)
gen->mfac[m] = gen->mfac[m-1]*sqrt((2*m+1.)/(2*m));
gen->root = RALLOC(double,2*gen->lmax+5);
gen->iroot = RALLOC(double,2*gen->lmax+5);
for (int m=0; m<2*gen->lmax+5; ++m)
gen->root = RALLOC(double,2*gen->lmax+6);
gen->iroot = RALLOC(double,2*gen->lmax+6);
for (int m=0; m<2*gen->lmax+6; ++m)
{
gen->root[m] = sqrt(m);
gen->iroot[m] = (m==0) ? 0. : 1./gen->root[m];
}
gen->eps=RALLOC(double, gen->lmax+10);
gen->alpha=RALLOC(double, gen->lmax/2+10);
gen->a=RALLOC(double, gen->lmax/2+10);
gen->b=RALLOC(double, gen->lmax/2+10);
}
else
{
@ -137,6 +141,10 @@ void sharp_Ylmgen_destroy (sharp_Ylmgen_C *gen)
DEALLOC(gen->mfac);
DEALLOC(gen->root);
DEALLOC(gen->iroot);
DEALLOC(gen->eps);
DEALLOC(gen->alpha);
DEALLOC(gen->a);
DEALLOC(gen->b);
}
else
{
@ -165,6 +173,18 @@ void sharp_Ylmgen_prepare (sharp_Ylmgen_C *gen, int m)
gen->rf[l].f[0] = tmp*gen->root[2*l+1];
gen->rf[l].f[1] = tmp*gen->root[l+m]*gen->root[l-m]*gen->iroot[2*l-1];
}
for (int l=m; l<gen->lmax+10; ++l)
gen->eps[l] = sqrt((l*l-m*m)/(4.*l*l-1));
gen->alpha[0] = 1./gen->eps[m+1];
gen->alpha[1] = gen->eps[m+1]/(gen->eps[m+2]*gen->eps[m+3]);
for (int il=1, l=m+2; l<gen->lmax+5; ++il, l+=2)
gen->alpha[il+1]= ((il&1) ? -1 : 1)/(gen->eps[l+2]*gen->eps[l+3]*gen->alpha[il]);
for (int il=0, l=m; l<gen->lmax+5; ++il, l+=2)
{
gen->a[il] = ((il&1) ? -1 : 1)*gen->alpha[il]*gen->alpha[il];
double t1 = gen->eps[l+2], t2 = gen->eps[l+1];
gen->b[il] = -gen->a[il]*(t1*t1+t2*t2);
}
}
else
{

2
libsharp/sharp_ylmgen_c.h
View file

@ -58,6 +58,8 @@ typedef struct
double *mfac;
sharp_ylmgen_dbl2 *rf;
double *eps, *alpha, *a, *b;
/* used if s!=0 */
int sinPow, cosPow, preMinus_p, preMinus_m;
double *prefac;