Skip to content

GitLab

Commit ee586c1f authored by Martin Reinecke's avatar Martin Reinecke
Browse files

better design for sharp_geom_info

parent 441c17ed
Hide whitespace changes
Inline Side-by-side
libsharp2/sharp.cc
View file @ ee586c1f
......@@ -45,9 +45,6 @@ using namespace mr;
using dcmplx = complex<double>;
using fcmplx = complex<float>;
static constexpr double sqrt_one_half = 0.707106781186547572737310929369;
static constexpr double sqrt_two = 1.414213562373095145474621858739;
static int chunksize_min=500, nchunks_max=10;
static void get_chunk_info (int ndata, int nmult, int &nchunks, int &chunksize)
......@@ -231,55 +228,53 @@ ptrdiff_t sharp_alm_info::index (int l, int mi)
}
sharp_geom_info::sharp_geom_info(size_t nrings, const size_t *nph, const ptrdiff_t *ofs,
const ptrdiff_t *stride, const double *phi0, const double *theta, const double *wgt)
ptrdiff_t stride_, const double *phi0, const double *theta, const double *wgt)
: ring(nrings), stride(stride_)
{
vector<Tring> infos(nrings);
size_t pos=0;
nphmax=0;
for (size_t m=0; m<nrings; ++m)
{
infos[m].theta = theta[m];
infos[m].cth = cos(theta[m]);
infos[m].sth = sin(theta[m]);
infos[m].weight = (wgt != nullptr) ? wgt[m] : 1.;
infos[m].phi0 = phi0[m];
infos[m].ofs = ofs[m];
infos[m].stride = stride[m];
infos[m].nph = nph[m];
ring[m].theta = theta[m];
ring[m].cth = cos(theta[m]);
ring[m].sth = sin(theta[m]);
ring[m].weight = (wgt != nullptr) ? wgt[m] : 1.;
ring[m].phi0 = phi0[m];
ring[m].ofs = ofs[m];
ring[m].nph = nph[m];
if (nphmax<nph[m]) nphmax=nph[m];
}
sort(infos.begin(), infos.end(),[](const Tring &a, const Tring &b)
sort(ring.begin(), ring.end(),[](const Tring &a, const Tring &b)
{ return (a.sth<b.sth); });
while (pos<nrings)
{
pair.push_back(Tpair());
pair.back().r1=infos[pos];
if ((pos<nrings-1) && approx(infos[pos].cth,-infos[pos+1].cth,1e-12))
pair.back().r1=pos;
if ((pos<nrings-1) && approx(ring[pos].cth,-ring[pos+1].cth,1e-12))
{
if (infos[pos].cth>0) // make sure northern ring is in r1
pair.back().r2=infos[pos+1];
if (ring[pos].cth>0) // make sure northern ring is in r1
pair.back().r2=pos+1;
else
{
pair.back().r1=infos[pos+1];
pair.back().r2=infos[pos];
pair.back().r1=pos+1;
pair.back().r2=pos;
}
++pos;
}
else
pair.back().r2.nph=0;
pair.back().r2=size_t(~0);
++pos;
}
sort(pair.begin(), pair.end(), [] (const Tpair &a, const Tpair &b)
sort(pair.begin(), pair.end(), [this] (const Tpair &a, const Tpair &b)
{
if (a.r1.nph==b.r1.nph)
return (a.r1.phi0 < b.r1.phi0) ? true :
((a.r1.phi0 > b.r1.phi0) ? false :
(a.r1.cth>b.r1.cth));
return a.r1.nph<b.r1.nph;
if (ring[a.r1].nph==ring[b.r1].nph)
return (ring[a.r1].phi0 < ring[b.r1].phi0) ? true :
((ring[a.r1].phi0 > ring[b.r1].phi0) ? false :
(ring[a.r1].cth>ring[b.r1].cth));
return ring[a.r1].nph<ring[b.r1].nph;
});
}
......@@ -302,34 +297,24 @@ static size_t sharp_get_mmax (const vector<size_t> &mval)
MRUTIL_NOINLINE void sharp_geom_info::clear_map (double *map) const
{
for (const auto &p: pair)
for (const auto &r: ring)
{
if (p.r1.stride==1)
memset(&map[p.r1.ofs],0,p.r1.nph*sizeof(double));
if (stride==1)
memset(&map[r.ofs],0,r.nph*sizeof(double));
else
for (size_t i=0;i<p.r1.nph;++i)
map[p.r1.ofs+i*p.r1.stride]=0;
if ((p.r2.nph>0)&&(p.r2.stride==1))
memset(&map[p.r2.ofs],0,p.r2.nph*sizeof(double));
else
for (size_t i=0;i<p.r2.nph;++i)
map[p.r2.ofs+i*p.r2.stride]=0;
for (size_t i=0;i<r.nph;++i)
map[r.ofs+i*stride]=0;
}
}
MRUTIL_NOINLINE void sharp_geom_info::clear_map (float *map) const
{
for (const auto &p: pair)
for (const auto &r: ring)
{
if (p.r1.stride==1)
memset(&map[p.r1.ofs],0,p.r1.nph*sizeof(float));
else
for (size_t i=0;i<p.r1.nph;++i)
map[p.r1.ofs+i*p.r1.stride]=0;
if ((p.r2.nph>0)&&(p.r2.stride==1))
memset(&map[p.r2.ofs],0,p.r2.nph*sizeof(float));
if (stride==1)
memset(&map[r.ofs],0,r.nph*sizeof(float));
else
for (size_t i=0;i<p.r2.nph;++i)
map[p.r2.ofs+i*p.r2.stride]=0;
for (size_t i=0;i<r.nph;++i)
map[r.ofs+i*stride]=0;
}
}
......@@ -481,56 +466,55 @@ MRUTIL_NOINLINE void sharp_job::almtmp2alm (size_t lmax, size_t mi)
}
}
MRUTIL_NOINLINE void sharp_job::ringtmp2ring (const sharp_geom_info::Tring &ri,
//virtual
void sharp_geom_info::add_ring(size_t iring, const double *ringtmp, double *map) const
{
double *MRUTIL_RESTRICT p1=&map[ring[iring].ofs];
for (size_t m=0; m<ring[iring].nph; ++m)
p1[m*stride] += ringtmp[m];
}
//virtual
void sharp_geom_info::add_ring(size_t iring, const double *ringtmp, float *map) const
{
float *MRUTIL_RESTRICT p1=&map[ring[iring].ofs];
for (size_t m=0; m<ring[iring].nph; ++m)
p1[m*stride] += float(ringtmp[m]);
}
//virtual
void sharp_geom_info::get_ring(size_t iring, const double *map, double *ringtmp) const
{
const double *MRUTIL_RESTRICT p1=&map[ring[iring].ofs];
for (size_t m=0; m<ring[iring].nph; ++m)
ringtmp[m] = p1[m*stride];
}
//virtual
void sharp_geom_info::get_ring(size_t iring, const float *map, double *ringtmp) const
{
const float *MRUTIL_RESTRICT p1=&map[ring[iring].ofs];
for (size_t m=0; m<ring[iring].nph; ++m)
ringtmp[m] = p1[m*stride];
}
MRUTIL_NOINLINE void sharp_job::ringtmp2ring (const sharp_geom_info &ginfo, size_t iring,
const vector<double> &ringtmp, ptrdiff_t rstride)
{
if (flags & SHARP_DP)
{
double **dmap = (double **)map;
for (size_t i=0; i<nmaps; ++i)
{
double *MRUTIL_RESTRICT p1=&dmap[i][ri.ofs];
const double *MRUTIL_RESTRICT p2=&ringtmp[i*rstride+1];
if (ri.stride==1)
{
if (flags&SHARP_ADD)
for (int m=0; m<ri.nph; ++m)
p1[m] += p2[m];
else
memcpy(p1,p2,ri.nph*sizeof(double));
}
else
for (size_t m=0; m<ri.nph; ++m)
p1[m*ri.stride] += p2[m];
}
}
ginfo.add_ring(iring, &ringtmp[i*rstride+1], ((double **)map)[i]);
else
{
float **fmap = (float **)map;
for (size_t i=0; i<nmaps; ++i)
for (size_t m=0; m<ri.nph; ++m)
fmap[i][ri.ofs+m*ri.stride] += (float)ringtmp[i*rstride+m+1];
}
ginfo.add_ring(iring, &ringtmp[i*rstride+1], ((float **)map)[i]);
}
MRUTIL_NOINLINE void sharp_job::ring2ringtmp (const sharp_geom_info::Tring &ri,
MRUTIL_NOINLINE void sharp_job::ring2ringtmp (const sharp_geom_info &ginfo, size_t iring,
vector<double> &ringtmp, ptrdiff_t rstride)
{
if (flags & SHARP_DP)
for (size_t i=0; i<nmaps; ++i)
{
double *MRUTIL_RESTRICT p1=&ringtmp[i*rstride+1],
*MRUTIL_RESTRICT p2=&(((double *)(map[i]))[ri.ofs]);
if (ri.stride==1)
memcpy(p1,p2,ri.nph*sizeof(double));
else
for (size_t m=0; m<ri.nph; ++m)
p1[m] = p2[m*ri.stride];
}
ginfo.get_ring(iring, ((double **)map)[i], &ringtmp[i*rstride+1]);
else
for (size_t i=0; i<nmaps; ++i)
for (size_t m=0; m<ri.nph; ++m)
ringtmp[i*rstride+m+1] = ((float *)(map[i]))[ri.ofs+m*ri.stride];
ginfo.get_ring(iring, ((float **)map)[i], &ringtmp[i*rstride+1]);
}
//FIXME: set phase to zero if not SHARP_MAP2ALM?
......@@ -547,15 +531,15 @@ MRUTIL_NOINLINE void sharp_job::map2phase (size_t mmax, size_t llim, size_t ulim
while (auto rng=sched.getNext()) for(auto ith=rng.lo+llim; ith<rng.hi+llim; ++ith)
{
int dim2 = s_th*(ith-llim);
ring2ringtmp(ginfo->pair[ith].r1,ringtmp,rstride);
ring2ringtmp(*ginfo, ginfo->pair[ith].r1,ringtmp,rstride);
for (size_t i=0; i<nmaps; ++i)
helper.ring2phase (ginfo->pair[ith].r1,
helper.ring2phase (ginfo->ring[ginfo->pair[ith].r1],
&ringtmp[i*rstride],mmax,&phase[dim2+2*i],pstride,flags);
if (ginfo->pair[ith].r2.nph>0)
if (ginfo->pair[ith].r2!=~size_t(0))
{
ring2ringtmp(ginfo->pair[ith].r2,ringtmp,rstride);
ring2ringtmp(*ginfo, ginfo->pair[ith].r2,ringtmp,rstride);
for (size_t i=0; i<nmaps; ++i)
helper.ring2phase (ginfo->pair[ith].r2,
helper.ring2phase (ginfo->ring[ginfo->pair[ith].r2],
&ringtmp[i*rstride],mmax,&phase[dim2+2*i+1],pstride,flags);
}
}
......@@ -576,15 +560,15 @@ MRUTIL_NOINLINE void sharp_job::phase2map (size_t mmax, size_t llim, size_t ulim
{
int dim2 = s_th*(ith-llim);
for (size_t i=0; i<nmaps; ++i)
helper.phase2ring (ginfo->pair[ith].r1,
helper.phase2ring (ginfo->ring[ginfo->pair[ith].r1],
&ringtmp[i*rstride],mmax,&phase[dim2+2*i],pstride,flags);
ringtmp2ring(ginfo->pair[ith].r1,ringtmp,rstride);
if (ginfo->pair[ith].r2.nph>0)
ringtmp2ring(*ginfo, ginfo->pair[ith].r1,ringtmp,rstride);
if (ginfo->pair[ith].r2!=~size_t(0))
{
for (size_t i=0; i<nmaps; ++i)
helper.phase2ring (ginfo->pair[ith].r2,
helper.phase2ring (ginfo->ring[ginfo->pair[ith].r2],
&ringtmp[i*rstride],mmax,&phase[dim2+2*i+1],pstride,flags);
ringtmp2ring(ginfo->pair[ith].r2,ringtmp,rstride);
ringtmp2ring(*ginfo, ginfo->pair[ith].r2,ringtmp,rstride);
}
}
}); /* end of parallel region */
......@@ -621,9 +605,9 @@ MRUTIL_NOINLINE void sharp_job::execute()
vector<double> cth(ulim-llim), sth(ulim-llim);
for (size_t i=0; i<ulim-llim; ++i)
{
ispair[i] = ginfo->pair[i+llim].r2.nph>0;
cth[i] = ginfo->pair[i+llim].r1.cth;
sth[i] = ginfo->pair[i+llim].r1.sth;
ispair[i] = ginfo->pair[i+llim].r2!=~size_t(0);
cth[i] = ginfo->ring[ginfo->pair[i+llim].r1].cth;
sth[i] = ginfo->ring[ginfo->pair[i+llim].r1].sth;
mlim[i] = sharp_get_mlim(lmax, spin, sth[i], cth[i]);
}
......
libsharp2/sharp.h
View file @ ee586c1f
......@@ -32,23 +32,23 @@
#include <vector>
#include <memory>
/*! \internal
Type holding all required information about a map geometry. */
struct sharp_geom_info
class sharp_geom_info
{
struct Tring
{
double theta, phi0, weight, cth, sth;
ptrdiff_t ofs;
size_t nph;
public:
struct Tring
{
double theta, phi0, weight, cth, sth;
ptrdiff_t ofs;
size_t nph;
};
struct Tpair
{
size_t r1, r2;
};
std::vector<Tring> ring;
std::vector<Tpair> pair;
ptrdiff_t stride;
};
struct Tpair
{
Tring r1,r2;
};
std::vector<Tpair> pair;
size_t nphmax;
size_t nphmax;
/*! Creates a geometry information from a set of ring descriptions.
All arrays passed to this function must have \a nrings elements.
\param nrings the number of rings in the map
......@@ -61,11 +61,15 @@ struct sharp_geom_info
and adjoint map2alm transforms.
Pass nullptr to use 1.0 as weight for all rings.
*/
sharp_geom_info(size_t nrings, const size_t *nph, const ptrdiff_t *ofs,
const ptrdiff_t *stride, const double *phi0, const double *theta,
const double *wgt);
void clear_map(double *map) const;
void clear_map(float *map) const;
sharp_geom_info(size_t nrings, const size_t *nph, const ptrdiff_t *ofs,
ptrdiff_t stride_, const double *phi0, const double *theta, const double *wgt);
virtual ~sharp_geom_info() {}
virtual void clear_map(double *map) const;
virtual void clear_map(float *map) const;
virtual void get_ring(size_t iring, const double *map, double *ringtmp) const;
virtual void get_ring(size_t iring, const float *map, double *ringtmp) const;
virtual void add_ring(size_t iring, const double *ringtmp, double *map) const;
virtual void add_ring(size_t iring, const double *ringtmp, float *map) const;
};
/*! \defgroup almgroup Helpers for dealing with a_lm */
......
libsharp2/sharp_geomhelpers.cc
View file @ ee586c1f
......@@ -43,13 +43,12 @@ unique_ptr<sharp_geom_info> sharp_make_subset_healpix_geom_info (size_t nside, p
vector<double> theta(nrings), weight_(nrings), phi0(nrings);
vector<size_t> nph(nrings);
vector<ptrdiff_t> stride_(nrings), ofs(nrings);
vector<ptrdiff_t> ofs(nrings);
ptrdiff_t curofs=0, checkofs; /* checkofs used for assertion introduced when adding rings arg */
for (size_t m=0; m<nrings; ++m)
{
auto ring = (rings==nullptr)? (m+1) : rings[m];
size_t northring = (ring>2*nside) ? 4*nside-ring : ring;
stride_[m] = stride;
if (northring < nside)
{
theta[m] = 2*asin(northring/(sqrt(6.)*nside));
......@@ -84,7 +83,7 @@ unique_ptr<sharp_geom_info> sharp_make_subset_healpix_geom_info (size_t nside, p
curofs+=nph[m];
}
return make_unique<sharp_geom_info>(nrings, nph.data(), ofs.data(), stride_.data(), phi0.data(), theta.data(), weight_.data());
return make_unique<sharp_geom_info>(nrings, nph.data(), ofs.data(), stride, phi0.data(), theta.data(), weight_.data());
}
unique_ptr<sharp_geom_info> sharp_make_weighted_healpix_geom_info (size_t nside, ptrdiff_t stride,
......@@ -99,7 +98,6 @@ unique_ptr<sharp_geom_info> sharp_make_gauss_geom_info (size_t nrings, size_t np
const double pi=3.141592653589793238462643383279502884197;
vector<size_t> nph(nrings);
vector<ptrdiff_t> stride_(nrings);
vector<double> phi0_(nrings);
vector<ptrdiff_t> ofs(nrings);
......@@ -112,11 +110,10 @@ unique_ptr<sharp_geom_info> sharp_make_gauss_geom_info (size_t nrings, size_t np
nph[m]=nphi;
phi0_[m]=phi0;
ofs[m]=ptrdiff_t(m*stride_lat);
stride_[m]=stride_lon;
weight[m]*=2*pi/nphi;
}
return make_unique<sharp_geom_info> (nrings, nph.data(), ofs.data(), stride_.data(), phi0_.data(), theta.data(), weight.data());
return make_unique<sharp_geom_info> (nrings, nph.data(), ofs.data(), stride_lon, phi0_.data(), theta.data(), weight.data());
}
/* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */
......@@ -127,7 +124,7 @@ unique_ptr<sharp_geom_info> sharp_make_fejer1_geom_info (size_t nrings, size_t p
vector<double> theta(nrings), weight(nrings), phi0_(nrings);
vector<size_t> nph(nrings);
vector<ptrdiff_t> stride_(nrings), ofs(nrings);
vector<ptrdiff_t> ofs(nrings);
weight[0]=2.;
for (size_t k=1; k<=(nrings-1)/2; ++k)
......@@ -146,11 +143,10 @@ unique_ptr<sharp_geom_info> sharp_make_fejer1_geom_info (size_t nrings, size_t p
phi0_[m]=phi0_[nrings-1-m]=phi0;
ofs[m]=ptrdiff_t(m*stride_lat);
ofs[nrings-1-m]=ptrdiff_t((nrings-1-m)*stride_lat);
stride_[m]=stride_[nrings-1-m]=stride_lon;
weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(nrings*nph[m]);
}
return make_unique<sharp_geom_info> (nrings, nph.data(), ofs.data(), stride_.data(), phi0_.data(), theta.data(), weight.data());
return make_unique<sharp_geom_info> (nrings, nph.data(), ofs.data(), stride_lon, phi0_.data(), theta.data(), weight.data());
}
/* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */
......@@ -161,7 +157,7 @@ unique_ptr<sharp_geom_info> sharp_make_cc_geom_info (size_t nrings, size_t pprin
vector<double> theta(nrings), weight(nrings,0.), phi0_(nrings);
vector<size_t> nph(nrings);
vector<ptrdiff_t> stride_(nrings), ofs(nrings);
vector<ptrdiff_t> ofs(nrings);
size_t n=nrings-1;
double dw=-1./(n*n-1.+(n&1));
......@@ -181,11 +177,10 @@ unique_ptr<sharp_geom_info> sharp_make_cc_geom_info (size_t nrings, size_t pprin
phi0_[m]=phi0_[nrings-1-m]=phi0;
ofs[m]=ptrdiff_t(m*stride_lat);
ofs[nrings-1-m]=ptrdiff_t((nrings-1-m)*stride_lat);
stride_[m]=stride_[nrings-1-m]=stride_lon;
weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(n*nph[m]);
}
return make_unique<sharp_geom_info> (nrings, nph.data(), ofs.data(), stride_.data(), phi0_.data(), theta.data(), weight.data());
return make_unique<sharp_geom_info> (nrings, nph.data(), ofs.data(), stride_lon, phi0_.data(), theta.data(), weight.data());
}
/* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */
......@@ -196,7 +191,7 @@ unique_ptr<sharp_geom_info> sharp_make_fejer2_geom_info (size_t nrings, size_t p
vector<double> theta(nrings), weight(nrings+1, 0.), phi0_(nrings);
vector<size_t> nph(nrings);
vector<ptrdiff_t> stride_(nrings), ofs(nrings);
vector<ptrdiff_t> ofs(nrings);
size_t n=nrings+1;
weight[0]=2.;
......@@ -215,11 +210,10 @@ unique_ptr<sharp_geom_info> sharp_make_fejer2_geom_info (size_t nrings, size_t p
phi0_[m]=phi0_[nrings-1-m]=phi0;
ofs[m]=ptrdiff_t(m*stride_lat);
ofs[nrings-1-m]=ptrdiff_t((nrings-1-m)*stride_lat);
stride_[m]=stride_[nrings-1-m]=stride_lon;
weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(n*nph[m]);
}
return make_unique<sharp_geom_info> (nrings, nph.data(), ofs.data(), stride_.data(), phi0_.data(), theta.data(), weight.data());
return make_unique<sharp_geom_info> (nrings, nph.data(), ofs.data(), stride_lon, phi0_.data(), theta.data(), weight.data());
}
unique_ptr<sharp_geom_info> sharp_make_mw_geom_info (size_t nrings, size_t ppring, double phi0,
......@@ -229,7 +223,7 @@ unique_ptr<sharp_geom_info> sharp_make_mw_geom_info (size_t nrings, size_t pprin
vector<double> theta(nrings), phi0_(nrings);
vector<size_t> nph(nrings);
vector<ptrdiff_t> stride_(nrings), ofs(nrings);
vector<ptrdiff_t> ofs(nrings);
for (size_t m=0; m<nrings; ++m)
{
......@@ -238,8 +232,7 @@ unique_ptr<sharp_geom_info> sharp_make_mw_geom_info (size_t nrings, size_t pprin
nph[m]=ppring;
phi0_[m]=phi0;
ofs[m]=ptrdiff_t(m*stride_lat);
stride_[m]=stride_lon;
}
return make_unique<sharp_geom_info> (nrings, nph.data(), ofs.data(), stride_.data(), phi0_.data(), theta.data(), nullptr);
return make_unique<sharp_geom_info> (nrings, nph.data(), ofs.data(), stride_lon, phi0_.data(), theta.data(), nullptr);
}
libsharp2/sharp_internal.h
View file @ ee586c1f
......@@ -58,9 +58,9 @@ struct sharp_job
void init_output();
void alm2almtmp (size_t lmax, size_t mi);
void almtmp2alm (size_t lmax, size_t mi);
void ring2ringtmp (const sharp_geom_info::Tring &ri, std::vector<double> &ringtmp,
void ring2ringtmp (const sharp_geom_info &ginfo, size_t iring, std::vector<double> &ringtmp,
ptrdiff_t rstride);
void ringtmp2ring (const sharp_geom_info::Tring &ri, const std::vector<double> &ringtmp, ptrdiff_t rstride);
void ringtmp2ring (const sharp_geom_info &ginfo, size_t iring, const std::vector<double> &ringtmp, ptrdiff_t rstride);
void map2phase (size_t mmax, size_t llim, size_t ulim);
void phase2map (size_t mmax, size_t llim, size_t ulim);
void execute();
......
mr_util/threading.h
View file @ ee586c1f
......@@ -41,6 +41,7 @@ struct Range
class Scheduler
{
public:
virtual ~Scheduler() {}
virtual size_t num_threads() const = 0;
virtual size_t thread_num() const = 0;
virtual Range getNext() = 0;
......
test/sharp2_testsuite.cc
View file @ ee586c1f
......@@ -137,11 +137,8 @@ static ptrdiff_t get_nalms(const sharp_alm_info &ainfo)
static ptrdiff_t get_npix(const sharp_geom_info &ginfo)
{
ptrdiff_t res=0;
for (size_t i=0; i<ginfo.pair.size(); ++i)
{
res += ginfo.pair[i].r1.nph;
if (ginfo.pair[i].r2.nph>0) res += ginfo.pair[i].r2.nph;
}
for (const auto &r: ginfo.ring)
res += r.nph;
return res;
}
......@@ -261,23 +258,15 @@ static void get_infos (const string &gname, int lmax, int &mmax, int &gpar1,
ginfo=sharp_make_gauss_geom_info (nlat, nlon, 0., 1, nlon);
ptrdiff_t npix_o=get_npix(*ginfo);
size_t ofs=0;
for (size_t i=0; i<ginfo->pair.size(); ++i)
for (auto &r:ginfo->ring)
{
auto &pair(ginfo->pair[i]);
int pring=1+2*sharp_get_mlim(lmax,0,pair.r1.sth,pair.r1.cth);
int pring=1+2*sharp_get_mlim(lmax,0,r.sth,r.cth);
if (pring>nlon) pring=nlon;
pring=good_fft_size(pring);
pair.r1.nph=pring;
pair.r1.weight*=nlon*1./pring;
pair.r1.ofs=ofs;
r.nph=pring;
r.weight*=nlon*1./pring;
r.ofs=ofs;
ofs+=pring;
if (pair.r2.nph>0)
{
pair.r2.nph=pring;
pair.r2.weight*=nlon*1./pring;
pair.r2.ofs=ofs;
ofs+=pring;
}
}
if (verbose)
{
......@@ -300,13 +289,11 @@ static void check_sign_scale(void)
auto tinfo = sharp_make_gauss_geom_info (nrings, ppring, 0., 1, ppring);
/* flip theta to emulate the "old" Gaussian grid geometry */
for (size_t i=0; i<tinfo->pair.size(); ++i)
for (auto &r: tinfo->ring)
{
const double pi=3.141592653589793238462643383279502884197;
tinfo->pair[i].r1.cth=-tinfo->pair[i].r1.cth;
tinfo->pair[i].r2.cth=-tinfo->pair[i].r2.cth;
tinfo->pair[i].r1.theta=pi-tinfo->pair[i].r1.theta;
tinfo->pair[i].r2.theta=pi-tinfo->pair[i].r2.theta;
r.cth=-r.cth;
r.theta=pi-r.theta;
}
auto alms = sharp_make_triangular_alm_info(lmax,mmax,1);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment