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

cleanups and fixes

Browse source
This commit is contained in:
Martin Reinecke 2018-12-20 13:34:48 +01:00
parent 3890bf174b
commit 04fbe296b2
3 changed files with 43 additions and 37 deletions
Download patch file Download diff file Expand all files Collapse all files

12
libsharp/sharp_core.c
View file

@ -219,15 +219,11 @@ NOINLINE static void alm2map_kernel(s0data_v * restrict d,
Tv a2=vload(ab[il+1].f[0]), b2=vload(ab[il+1].f[1]);
for (int i=0; i<nv2; ++i)
{
d->p1r[i] += d->lam2[i]*ar1;
d->p1i[i] += d->lam2[i]*ai1;
d->p2r[i] += d->lam2[i]*ar2;
d->p2i[i] += d->lam2[i]*ai2;
d->lam1[i] = (a1*d->csq[i] + b1)*d->lam2[i] + d->lam1[i];
d->p1r[i] += d->lam1[i]*ar3;
d->p1i[i] += d->lam1[i]*ai3;
d->p2r[i] += d->lam1[i]*ar4;
d->p2i[i] += d->lam1[i]*ai4;
d->p1r[i] += d->lam2[i]*ar1 + d->lam1[i]*ar3;
d->p1i[i] += d->lam2[i]*ai1 + d->lam1[i]*ai3;
d->p2r[i] += d->lam2[i]*ar2 + d->lam1[i]*ar4;
d->p2i[i] += d->lam2[i]*ai2 + d->lam1[i]*ai4;
d->lam2[i] = (a2*d->csq[i] + b2)*d->lam1[i] + d->lam2[i];
}
}

49
libsharp/sharp_testsuite.c
View file

@ -241,13 +241,14 @@ static int good_fft_size(int n)
}
static void get_infos (const char *gname, int lmax, int *mmax, int *gpar1,
int *gpar2, sharp_geom_info **ginfo, sharp_alm_info **ainfo)
int *gpar2, sharp_geom_info **ginfo, sharp_alm_info **ainfo, int verbose)
{
UTIL_ASSERT(lmax>=0,"lmax must not be negative");
if (*mmax<0) *mmax=lmax;
UTIL_ASSERT(*mmax<=lmax,"mmax larger than lmax");
if (mytask==0) printf ("lmax: %d, mmax: %d\n",lmax,*mmax);
verbose &= (mytask==0);
if (verbose) printf ("lmax: %d, mmax: %d\n",lmax,*mmax);
sharp_make_triangular_alm_info(lmax,*mmax,1,ainfo);
#ifdef USE_MPI
@ -259,14 +260,14 @@ static void get_infos (const char *gname, int lmax, int *mmax, int *gpar1,
if (*gpar1<1) *gpar1=lmax/2;
if (*gpar1==0) ++(*gpar1);
sharp_make_healpix_geom_info (*gpar1, 1, ginfo);
if (mytask==0) printf ("HEALPix grid, nside=%d\n",*gpar1);
if (verbose) printf ("HEALPix grid, nside=%d\n",*gpar1);
}
else if (strcmp(gname,"gauss")==0)
{
if (*gpar1<1) *gpar1=lmax+1;
if (*gpar2<1) *gpar2=2*(*mmax)+1;
sharp_make_gauss_geom_info (*gpar1, *gpar2, 0., 1, *gpar2, ginfo);
if (mytask==0)
if (verbose)
printf ("Gauss-Legendre grid, nlat=%d, nlon=%d\n",*gpar1,*gpar2);
}
else if (strcmp(gname,"fejer1")==0)
@ -274,21 +275,21 @@ static void get_infos (const char *gname, int lmax, int *mmax, int *gpar1,
if (*gpar1<1) *gpar1=2*lmax+1;
if (*gpar2<1) *gpar2=2*(*mmax)+1;
sharp_make_fejer1_geom_info (*gpar1, *gpar2, 0., 1, *gpar2, ginfo);
if (mytask==0) printf ("Fejer1 grid, nlat=%d, nlon=%d\n",*gpar1,*gpar2);
if (verbose) printf ("Fejer1 grid, nlat=%d, nlon=%d\n",*gpar1,*gpar2);
}
else if (strcmp(gname,"fejer2")==0)
{
if (*gpar1<1) *gpar1=2*lmax+1;
if (*gpar2<1) *gpar2=2*(*mmax)+1;
sharp_make_fejer2_geom_info (*gpar1, *gpar2, 0., 1, *gpar2, ginfo);
if (mytask==0) printf ("Fejer2 grid, nlat=%d, nlon=%d\n",*gpar1,*gpar2);
if (verbose) printf ("Fejer2 grid, nlat=%d, nlon=%d\n",*gpar1,*gpar2);
}
else if (strcmp(gname,"cc")==0)
{
if (*gpar1<1) *gpar1=2*lmax+1;
if (*gpar2<1) *gpar2=2*(*mmax)+1;
sharp_make_cc_geom_info (*gpar1, *gpar2, 0., 1, *gpar2, ginfo);
if (mytask==0)
if (verbose)
printf("Clenshaw-Curtis grid, nlat=%d, nlon=%d\n",*gpar1,*gpar2);
}
else if (strcmp(gname,"smallgauss")==0)
@ -318,7 +319,7 @@ static void get_infos (const char *gname, int lmax, int *mmax, int *gpar1,
ofs+=pring;
}
}
if (mytask==0)
if (verbose)
{
ptrdiff_t npix=get_npix(*ginfo);
printf("Small Gauss grid, nlat=%d, npix=%ld, savings=%.2f%%\n",
@ -485,6 +486,16 @@ static void check_accuracy (sharp_geom_info *ginfo, sharp_alm_info *ainfo,
DEALLOC(err_abs);
}
static void run(int lmax, int mmax, int nlat, int nlon, int spin)
{
sharp_geom_info *ginfo;
sharp_alm_info *ainfo;
get_infos ("gauss", lmax, &mmax, &nlat, &nlon, &ginfo, &ainfo, 0);
check_accuracy(ginfo,ainfo,spin);
sharp_destroy_alm_info(ainfo);
sharp_destroy_geom_info(ginfo);
}
static void sharp_acctest(void)
{
if (mytask==0) sharp_module_startup("sharp_acctest",1,1,"",1);
@ -495,17 +506,15 @@ static void sharp_acctest(void)
if (mytask==0) printf("Testing map analysis accuracy.\n");
sharp_geom_info *ginfo;
sharp_alm_info *ainfo;
int lmax=127, mmax=127, nlat=128, nlon=256;
get_infos ("gauss", lmax, &mmax, &nlat, &nlon, &ginfo, &ainfo);
check_accuracy(ginfo,ainfo,0);
check_accuracy(ginfo,ainfo,1);
check_accuracy(ginfo,ainfo,2);
check_accuracy(ginfo,ainfo,3);
check_accuracy(ginfo,ainfo,30);
sharp_destroy_alm_info(ainfo);
sharp_destroy_geom_info(ginfo);
run(127, 127, 128, 256, 0);
run(127, 127, 128, 256, 1);
run(127, 127, 128, 256, 2);
run(127, 127, 128, 256, 3);
run(127, 127, 128, 256, 30);
run(5, 0, 6, 1, 0);
run(5, 0, 7, 2, 0);
run(8, 8, 9, 17, 0);
run(8, 8, 9, 17, 2);
if (mytask==0) printf("Passed.\n\n");
}
@ -524,7 +533,7 @@ static void sharp_test (int argc, const char **argv)
sharp_geom_info *ginfo;
sharp_alm_info *ainfo;
get_infos (argv[2], lmax, &mmax, &gpar1, &gpar2, &ginfo, &ainfo);
get_infos (argv[2], lmax, &mmax, &gpar1, &gpar2, &ginfo, &ainfo, 1);
int ncomp = (spin==0) ? 1 : 2;
double t_a2m=1e30, t_m2a=1e30;

19
libsharp/sharp_ylmgen_c.c
View file

@ -73,16 +73,16 @@ 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+7);
gen->iroot = RALLOC(double,2*gen->lmax+7);
for (int m=0; m<2*gen->lmax+7; ++m)
gen->root = RALLOC(double,2*gen->lmax+8);
gen->iroot = RALLOC(double,2*gen->lmax+8);
for (int m=0; m<2*gen->lmax+8; ++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->ab=RALLOC(sharp_ylmgen_dbl2, gen->lmax/2+10);
gen->eps=RALLOC(double, gen->lmax+4);
gen->alpha=RALLOC(double, gen->lmax/2+2);
gen->ab=RALLOC(sharp_ylmgen_dbl2, gen->lmax/2+2);
}
else
{
@ -161,15 +161,16 @@ void sharp_Ylmgen_prepare (sharp_Ylmgen_C *gen, int m)
if (gen->s==0)
{
for (int l=m; l<gen->lmax+10; ++l)
gen->eps[m] = 0.;
for (int l=m+1; l<gen->lmax+4; ++l)
gen->eps[l] = gen->root[l+m]*gen->root[l-m]
*gen->iroot[2*l+1]*gen->iroot[2*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)
for (int il=1, l=m+2; l<gen->lmax+1; ++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)
for (int il=0, l=m; l<gen->lmax+2; ++il, l+=2)
{
gen->ab[il].f[0] = ((il&1) ? -1 : 1)*gen->alpha[il]*gen->alpha[il];
double t1 = gen->eps[l+2], t2 = gen->eps[l+1];