Skip to content

GitLab

Commit 0fb8e5bd authored by Martin Reinecke's avatar Martin Reinecke
Browse files

cleanup

parent b8c81cad
Hide whitespace changes
Inline Side-by-side
libsharp2/sharp.cc
View file @ 0fb8e5bd
......@@ -110,6 +110,110 @@ struct ringhelper
length=nph;
}
}
MRUTIL_NOINLINE void phase2ring (const sharp_ringinfo &info,
double *data, int mmax, const dcmplx *phase, int pstride, int flags)
{
int nph = info.nph;
update (nph, mmax, info.phi0);
double wgt = (flags&SHARP_USE_WEIGHTS) ? info.weight : 1.;
if (flags&SHARP_REAL_HARMONICS)
wgt *= sqrt_one_half;
if (nph>=2*mmax+1)
{
if (norot)
for (int m=0; m<=mmax; ++m)
{
data[2*m]=phase[m*pstride].real()*wgt;
data[2*m+1]=phase[m*pstride].imag()*wgt;
}
else
for (int m=0; m<=mmax; ++m)
{
dcmplx tmp = phase[m*pstride]*shiftarr[m];
data[2*m]=tmp.real()*wgt;
data[2*m+1]=tmp.imag()*wgt;
}
for (int m=2*(mmax+1); m<nph+2; ++m)
data[m]=0.;
}
else
{
data[0]=phase[0].real()*wgt;
fill(data+1,data+nph+2,0.);
int idx1=1, idx2=nph-1;
for (int m=1; m<=mmax; ++m)
{
dcmplx tmp = phase[m*pstride]*wgt;
if(!norot) tmp*=shiftarr[m];
if (idx1<(nph+2)/2)
{
data[2*idx1]+=tmp.real();
data[2*idx1+1]+=tmp.imag();
}
if (idx2<(nph+2)/2)
{
data[2*idx2]+=tmp.real();
data[2*idx2+1]-=tmp.imag();
}
if (++idx1>=nph) idx1=0;
if (--idx2<0) idx2=nph-1;
}
}
data[1]=data[0];
plan->exec(&(data[1]), 1., false);
}
MRUTIL_NOINLINE void ring2phase (const sharp_ringinfo &info,
double *data, int mmax, dcmplx *phase, int pstride, int flags)
{
int nph = info.nph;
#if 1
int maxidx = mmax; /* Enable this for traditional Healpix compatibility */
#else
int maxidx = min(nph-1,mmax);
#endif
update (nph, mmax, -info.phi0);
double wgt = (flags&SHARP_USE_WEIGHTS) ? info.weight : 1;
if (flags&SHARP_REAL_HARMONICS)
wgt *= sqrt_two;
plan->exec (&(data[1]), 1., true);
data[0]=data[1];
data[1]=data[nph+1]=0.;
if (maxidx<=nph/2)
{
if (norot)
for (int m=0; m<=maxidx; ++m)
phase[m*pstride] = dcmplx(data[2*m], data[2*m+1]) * wgt;
else
for (int m=0; m<=maxidx; ++m)
phase[m*pstride] =
dcmplx(data[2*m], data[2*m+1]) * shiftarr[m] * wgt;
}
else
{
for (int m=0; m<=maxidx; ++m)
{
int idx=m%nph;
dcmplx val;
if (idx<(nph-idx))
val = dcmplx(data[2*idx], data[2*idx+1]) * wgt;
else
val = dcmplx(data[2*(nph-idx)], -data[2*(nph-idx)+1]) * wgt;
if (!norot)
val *= shiftarr[m];
phase[m*pstride]=val;
}
}
for (int m=maxidx+1;m<=mmax; ++m)
phase[m*pstride]=0.;
}
};
sharp_alm_info::sharp_alm_info (int lmax, int nm, int stride, const int *mval,
......@@ -141,54 +245,51 @@ ptrdiff_t sharp_alm_info::index (int l, int mi)
return mvstart[mi]+stride*l;
}
unique_ptr<sharp_geom_info> sharp_make_geom_info (int nrings, const int *nph, const ptrdiff_t *ofs,
const int *stride, const double *phi0, const double *theta,
const double *wgt)
sharp_geom_info::sharp_geom_info(int nrings, const int *nph, const ptrdiff_t *ofs,
const int *stride, const double *phi0, const double *theta, const double *wgt)
{
sharp_geom_info *info = new sharp_geom_info;
vector<sharp_ringinfo> infos(nrings);
int pos=0;
info->pair.resize(nrings);
int npairs=0;
info->nphmax=0;
nphmax=0;
for (int 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 != NULL) ? wgt[m] : 1.;
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];
if (info->nphmax<nph[m]) info->nphmax=nph[m];
if (nphmax<nph[m]) nphmax=nph[m];
}
sort(infos.begin(), infos.end(),[](const sharp_ringinfo &a, const sharp_ringinfo &b)
{ return (a.sth<b.sth); });
while (pos<nrings)
{
info->pair[npairs].r1=infos[pos];
pair.push_back(sharp_ringpair());
pair.back().r1=infos[pos];
if ((pos<nrings-1) && approx(infos[pos].cth,-infos[pos+1].cth,1e-12))
{
if (infos[pos].cth>0) // make sure northern ring is in r1
info->pair[npairs].r2=infos[pos+1];
pair.back().r2=infos[pos+1];
else
{
info->pair[npairs].r1=infos[pos+1];
info->pair[npairs].r2=infos[pos];
pair.back().r1=infos[pos+1];
pair.back().r2=infos[pos];
}
++pos;
}
else
info->pair[npairs].r2.nph=-1;
pair.back().r2.nph=-1;
++pos;
++npairs;
}
info->pair.resize(npairs);
sort(info->pair.begin(), info->pair.end(), [] (const sharp_ringpair &a, const sharp_ringpair &b)
sort(pair.begin(), pair.end(), [] (const sharp_ringpair &a, const sharp_ringpair &b)
{
if (a.r1.nph==b.r1.nph)
return (a.r1.phi0 < b.r1.phi0) ? true :
......@@ -196,7 +297,6 @@ unique_ptr<sharp_geom_info> sharp_make_geom_info (int nrings, const int *nph, co
(a.r1.cth>b.r1.cth));
return a.r1.nph<b.r1.nph;
});
return unique_ptr<sharp_geom_info>(info);
}
/* This currently requires all m values from 0 to nm-1 to be present.
......@@ -217,114 +317,6 @@ static int sharp_get_mmax (const vector<int> &mval)
return nm-1;
}
MRUTIL_NOINLINE static void ringhelper_phase2ring (ringhelper &self,
const sharp_ringinfo &info, double *data, int mmax, const dcmplx *phase,
int pstride, int flags)
{
int nph = info.nph;
self.update (nph, mmax, info.phi0);
double wgt = (flags&SHARP_USE_WEIGHTS) ? info.weight : 1.;
if (flags&SHARP_REAL_HARMONICS)
wgt *= sqrt_one_half;
if (nph>=2*mmax+1)
{
if (self.norot)
for (int m=0; m<=mmax; ++m)
{
data[2*m]=phase[m*pstride].real()*wgt;
data[2*m+1]=phase[m*pstride].imag()*wgt;
}
else
for (int m=0; m<=mmax; ++m)
{
dcmplx tmp = phase[m*pstride]*self.shiftarr[m];
data[2*m]=tmp.real()*wgt;
data[2*m+1]=tmp.imag()*wgt;
}
for (int m=2*(mmax+1); m<nph+2; ++m)
data[m]=0.;
}
else
{
data[0]=phase[0].real()*wgt;
fill(data+1,data+nph+2,0.);
int idx1=1, idx2=nph-1;
for (int m=1; m<=mmax; ++m)
{
dcmplx tmp = phase[m*pstride]*wgt;
if(!self.norot) tmp*=self.shiftarr[m];
if (idx1<(nph+2)/2)
{
data[2*idx1]+=tmp.real();
data[2*idx1+1]+=tmp.imag();
}
if (idx2<(nph+2)/2)
{
data[2*idx2]+=tmp.real();
data[2*idx2+1]-=tmp.imag();
}
if (++idx1>=nph) idx1=0;
if (--idx2<0) idx2=nph-1;
}
}
data[1]=data[0];
self.plan->exec(&(data[1]), 1., false);
}
MRUTIL_NOINLINE static void ringhelper_ring2phase (ringhelper &self,
const sharp_ringinfo &info, double *data, int mmax, dcmplx *phase,
int pstride, int flags)
{
int nph = info.nph;
#if 1
int maxidx = mmax; /* Enable this for traditional Healpix compatibility */
#else
int maxidx = min(nph-1,mmax);
#endif
self.update (nph, mmax, -info.phi0);
double wgt = (flags&SHARP_USE_WEIGHTS) ? info.weight : 1;
if (flags&SHARP_REAL_HARMONICS)
wgt *= sqrt_two;
self.plan->exec (&(data[1]), 1., true);
data[0]=data[1];
data[1]=data[nph+1]=0.;
if (maxidx<=nph/2)
{
if (self.norot)
for (int m=0; m<=maxidx; ++m)
phase[m*pstride] = dcmplx(data[2*m], data[2*m+1]) * wgt;
else
for (int m=0; m<=maxidx; ++m)
phase[m*pstride] =
dcmplx(data[2*m], data[2*m+1]) * self.shiftarr[m] * wgt;
}
else
{
for (int m=0; m<=maxidx; ++m)
{
int idx=m%nph;
dcmplx val;
if (idx<(nph-idx))
val = dcmplx(data[2*idx], data[2*idx+1]) * wgt;
else
val = dcmplx(data[2*(nph-idx)], -data[2*(nph-idx)+1]) * wgt;
if (!self.norot)
val *= self.shiftarr[m];
phase[m*pstride]=val;
}
}
for (int m=maxidx+1;m<=mmax; ++m)
phase[m*pstride]=0.;
}
MRUTIL_NOINLINE static void clear_map (const sharp_geom_info *ginfo, void *map,
int flags)
{
......@@ -711,14 +703,14 @@ MRUTIL_NOINLINE static void map2phase (sharp_job &job, int mmax, int llim, int u
int dim2 = job.s_th*(ith-llim);
ring2ringtmp(job,job.ginfo->pair[ith].r1,ringtmp,rstride);
for (int i=0; i<job.nmaps; ++i)
ringhelper_ring2phase (helper,job.ginfo->pair[ith].r1,
helper.ring2phase (job.ginfo->pair[ith].r1,
&ringtmp[i*rstride],mmax,&job.phase[dim2+2*i],pstride,job.flags);
if (job.ginfo->pair[ith].r2.nph>0)
{
ring2ringtmp(job,job.ginfo->pair[ith].r2,ringtmp,rstride);
for (int i=0; i<job.nmaps; ++i)
ringhelper_ring2phase (helper,job.ginfo->pair[ith].r2,
&ringtmp[i*rstride],mmax,&job.phase[dim2+2*i+1],pstride,job.flags);
helper.ring2phase (job.ginfo->pair[ith].r2,
&ringtmp[i*rstride],mmax,&job.phase[dim2+2*i+1],pstride,job.flags);
}
}
}); /* end of parallel region */
......@@ -752,13 +744,13 @@ MRUTIL_NOINLINE static void phase2map (sharp_job &job, int mmax, int llim, int u
{
int dim2 = job.s_th*(ith-llim);
for (int i=0; i<job.nmaps; ++i)
ringhelper_phase2ring (helper,job.ginfo->pair[ith].r1,
helper.phase2ring (job.ginfo->pair[ith].r1,
&ringtmp[i*rstride],mmax,&job.phase[dim2+2*i],pstride,job.flags);
ringtmp2ring(job,job.ginfo->pair[ith].r1,ringtmp,rstride);
if (job.ginfo->pair[ith].r2.nph>0)
{
for (int i=0; i<job.nmaps; ++i)
ringhelper_phase2ring (helper,job.ginfo->pair[ith].r2,
helper.phase2ring (job.ginfo->pair[ith].r2,
&ringtmp[i*rstride],mmax,&job.phase[dim2+2*i+1],pstride,job.flags);
ringtmp2ring(job,job.ginfo->pair[ith].r2,ringtmp,rstride);
}
......@@ -871,8 +863,8 @@ void sharp_execute (sharp_jobtype type, int spin, void *alm, void *map,
flags);
sharp_execute_job (job);
if (time!=NULL) *time = job.time;
if (opcnt!=NULL) *opcnt = job.opcnt;
if (time!=nullptr) *time = job.time;
if (opcnt!=nullptr) *opcnt = job.opcnt;
}
void sharp_set_chunksize_min(int new_chunksize_min)
......
libsharp2/sharp.h
View file @ 0fb8e5bd
......@@ -55,6 +55,21 @@ struct sharp_geom_info
{
std::vector<sharp_ringpair> pair;
int 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
\param nph the number of pixels in each ring
\param ofs the index of the first pixel in each ring in the map array
\param stride the stride between consecutive pixels
\param phi0 the azimuth (in radians) of the first pixel in each ring
\param theta the colatitude (in radians) of each ring
\param wgt the pixel weight to be used for the ring in map2alm
and adjoint map2alm transforms.
Pass nullptr to use 1.0 as weight for all rings.
*/
sharp_geom_info(int nrings, const int *nph, const ptrdiff_t *ofs,
const int *stride, const double *phi0, const double *theta,
const double *wgt);
};
/*! \defgroup almgroup Helpers for dealing with a_lm */
......@@ -130,23 +145,6 @@ enum sharp_almflags { SHARP_PACKED = 1,
/*! \defgroup geominfogroup Functions for dealing with geometry information */
/*! \{ */
/*! 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
\param nph the number of pixels in each ring
\param ofs the index of the first pixel in each ring in the map array
\param stride the stride between consecutive pixels
\param phi0 the azimuth (in radians) of the first pixel in each ring
\param theta the colatitude (in radians) of each ring
\param wgt the pixel weight to be used for the ring in map2alm
and adjoint map2alm transforms.
Pass NULL to use 1.0 as weight for all rings.
\param geom_info will hold a pointer to the newly created data structure
*/
std::unique_ptr<sharp_geom_info> sharp_make_geom_info (int nrings, const int *nph, const ptrdiff_t *ofs,
const int *stride, const double *phi0, const double *theta,
const double *wgt);
/*! \} */
/*! \defgroup lowlevelgroup Low-level libsharp2 SHT interface */
......@@ -198,9 +196,9 @@ enum sharp_jobflags { SHARP_DP = 1<<4,
\a alm is expected to have the type "complex double **" and \a map is
expected to have the type "double **"; otherwise, the expected
types are "complex float **" and "float **", respectively.
\param time If not NULL, the wall clock time required for this SHT
\param time If not nullptr, the wall clock time required for this SHT
(in seconds) will be written here.
\param opcnt If not NULL, a conservative estimate of the total floating point
\param opcnt If not nullptr, a conservative estimate of the total floating point
operation count for this SHT will be written here. */
void sharp_execute (sharp_jobtype type, int spin, void *alm, void *map,
const sharp_geom_info &geom_info, const sharp_alm_info &alm_info,
......
libsharp2/sharp_core.cc
View file @ 0fb8e5bd
......@@ -38,10 +38,10 @@ using t_veclen = int (*) (void);
using t_max_nvec = int (*) (int spin);
using t_architecture = const char *(*) (void);
static t_inner_loop inner_loop_ = NULL;
static t_veclen veclen_ = NULL;
static t_max_nvec max_nvec_ = NULL;
static t_architecture architecture_ = NULL;
static t_inner_loop inner_loop_ = nullptr;
static t_veclen veclen_ = nullptr;
static t_max_nvec max_nvec_ = nullptr;
static t_architecture architecture_ = nullptr;
#ifdef MULTIARCH
......
libsharp2/sharp_geomhelpers.cc
View file @ 0fb8e5bd
......@@ -47,7 +47,7 @@ unique_ptr<sharp_geom_info> sharp_make_subset_healpix_geom_info (int nside, int
ptrdiff_t curofs=0, checkofs; /* checkofs used for assertion introduced when adding rings arg */
for (int m=0; m<nrings; ++m)
{
int ring = (rings==NULL)? (m+1) : rings[m];
int ring = (rings==nullptr)? (m+1) : rings[m];
ptrdiff_t northring = (ring>2*nside) ? 4*nside-ring : ring;
stride_[m] = stride;
if (northring < nside)
......@@ -76,21 +76,21 @@ unique_ptr<sharp_geom_info> sharp_make_subset_healpix_geom_info (int nside, int
checkofs = (npix - nph[m])*stride - checkofs;
ofs[m] = curofs;
}
weight_[m]=4.*pi/npix*((weight==NULL) ? 1. : weight[northring-1]);
if (rings==NULL) {
weight_[m]=4.*pi/npix*((weight==nullptr) ? 1. : weight[northring-1]);
if (rings==nullptr) {
MR_assert(curofs==checkofs, "Bug in computing ofs[m]");
}
ofs[m] = curofs;
curofs+=nph[m];
}
return sharp_make_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_.data(), phi0.data(), theta.data(), weight_.data());
}
unique_ptr<sharp_geom_info> sharp_make_weighted_healpix_geom_info (int nside, int stride,
const double *weight)
{
return sharp_make_subset_healpix_geom_info(nside, stride, 4 * nside - 1, NULL, weight);
return sharp_make_subset_healpix_geom_info(nside, stride, 4 * nside - 1, nullptr, weight);
}
unique_ptr<sharp_geom_info> sharp_make_gauss_geom_info (int nrings, int nphi, double phi0,
......@@ -115,7 +115,7 @@ unique_ptr<sharp_geom_info> sharp_make_gauss_geom_info (int nrings, int nphi, do
weight[m]*=2*pi/nphi;
}
return sharp_make_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_.data(), phi0_.data(), theta.data(), weight.data());
}
/* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */
......@@ -149,7 +149,7 @@ unique_ptr<sharp_geom_info> sharp_make_fejer1_geom_info (int nrings, int ppring,
weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(nrings*nph[m]);
}
return sharp_make_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_.data(), phi0_.data(), theta.data(), weight.data());
}
/* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */
......@@ -184,7 +184,7 @@ unique_ptr<sharp_geom_info> sharp_make_cc_geom_info (int nrings, int ppring, dou
weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(n*nph[m]);
}
return sharp_make_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_.data(), phi0_.data(), theta.data(), weight.data());
}
/* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */
......@@ -218,7 +218,7 @@ unique_ptr<sharp_geom_info> sharp_make_fejer2_geom_info (int nrings, int ppring,
weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(n*nph[m]);
}
return sharp_make_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_.data(), phi0_.data(), theta.data(), weight.data());
}
unique_ptr<sharp_geom_info> sharp_make_mw_geom_info (int nrings, int ppring, double phi0,
......@@ -240,5 +240,5 @@ unique_ptr<sharp_geom_info> sharp_make_mw_geom_info (int nrings, int ppring, dou
stride_[m]=stride_lon;
}
return sharp_make_geom_info (nrings, nph.data(), ofs.data(), stride_.data(), phi0_.data(), theta.data(), NULL);
return make_unique<sharp_geom_info> (nrings, nph.data(), ofs.data(), stride_.data(), phi0_.data(), theta.data(), nullptr);
}
libsharp2/sharp_geomhelpers.h
View file @ 0fb8e5bd
......@@ -35,8 +35,8 @@
Nside parameter \a nside. \a weight contains the relative ring
weights and must have \a 2*nside entries. The rings array contains
the indices of the rings, with 1 being the first ring at the north
pole; if NULL then we take them to be sequential. Pass 4 * nside - 1
as nrings and NULL to rings to get the full HEALPix grid.
pole; if nullptr then we take them to be sequential. Pass 4 * nside - 1
as nrings and nullptr to rings to get the full HEALPix grid.
\note if \a weight is a null pointer, all weights are assumed to be 1.
\note if \a rings is a null pointer, take all rings
\ingroup geominfogroup */
......@@ -55,7 +55,7 @@ std::unique_ptr<sharp_geom_info> sharp_make_weighted_healpix_geom_info (int nsid
Nside parameter \a nside.
\ingroup geominfogroup */
static inline std::unique_ptr<sharp_geom_info> sharp_make_healpix_geom_info (int nside, int stride)
{ return sharp_make_weighted_healpix_geom_info (nside, stride, NULL); }
{ return sharp_make_weighted_healpix_geom_info (nside, stride, nullptr); }
/*! Creates a geometry information describing a Gaussian map with \a nrings
iso-latitude rings and \a nphi pixels per ring. The azimuth of the first
......
test/sharp2_testsuite.cc
View file @ 0fb8e5bd
......@@ -322,7 +322,7 @@ static void check_sign_scale(void)
vector<dcmplx *>alm({&balm[0], &balm[nalms]});
sharp_execute(SHARP_ALM2MAP,0,alm.data(),map.data(),*tinfo,*alms,SHARP_DP,
NULL,NULL);
nullptr,nullptr);
MR_assert(approx(map[0][0 ], 3.588246976618616912e+00,1e-12),
"error");
MR_assert(approx(map[0][npix/2], 4.042209792157496651e+01,1e-12),
......@@ -331,7 +331,7 @@ static void check_sign_scale(void)
"error");
sharp_execute(SHARP_ALM2MAP,1,alm.data(),map.data(),*tinfo,*alms,SHARP_DP,
NULL,NULL);
nullptr,nullptr);
MR_assert(approx(map[0][0 ], 2.750897760535633285e+00,1e-12),
"error");
MR_assert(approx(map[0][npix/2], 3.137704477368562905e+01,1e-12),
......@@ -346,7 +346,7 @@ static void check_sign_scale(void)
"error");
sharp_execute(SHARP_ALM2MAP,2,alm.data(),map.data(),*tinfo,*alms,SHARP_DP,
NULL,NULL);
nullptr,nullptr);
MR_assert(approx(map[0][0 ],-1.398186224727334448e+00,1e-12),
"error");
MR_assert(approx(map[0][npix/2],-2.456676000884031197e+01,1e-12),
......@@ -361,7 +361,7 @@ static void check_sign_scale(void)
"error");
sharp_execute(SHARP_ALM2MAP_DERIV1,1,alm.data(),map.data(),*tinfo,*alms,
SHARP_DP,NULL,NULL);
SHARP_DP,nullptr,nullptr);
MR_assert(approx(map[0][0 ],-6.859393905369091105e-01,1e-11),
"error");
MR_assert(approx(map[0][npix/2],-2.103947835973212364e+02,1e-12),
......@@ -401,24 +401,24 @@ static void do_sht (sharp_geom_info &ginfo, sharp_alm_info &ainfo,
double tta2m, ttm2a;
unsigned long long toa2m, tom2a;
if (t_a2m!=NULL) *t_a2m=0;
if (op_a2m!=NULL) *op_a2m=0;
if (t_a2m!=nullptr) *t_a2m=0;
if (op_a2m!=nullptr) *op_a2m=0;
for (size_t itrans=0; itrans<ntrans; ++itrans)
{
sharp_execute(SHARP_ALM2MAP,spin,&alm[itrans*ncomp],&map[itrans*ncomp],ginfo,ainfo,
SHARP_DP,&tta2m,&toa2m);
if (t_a2m!=NULL) *t_a2m+=maxTime(tta2m);
if (op_a2m!=NULL) *op_a2m+=totalops(toa2m);
if (t_a2m!=nullptr) *t_a2m+=maxTime(tta2m);
if (op_a2m!=nullptr) *op_a2m+=totalops(toa2m);
}
auto sqsum=get_sqsum_and_invert(alm.data(),nalms,ntrans*ncomp);
if (t_m2a!=NULL) *t_m2a=0;
if (op_m2a!=NULL) *op_m2a=0;
if (t_m2a!=nullptr) *t_m2a=0;
if (op_m2a!=nullptr) *op_m2a=0;
for (size_t itrans=0; itrans<ntrans; ++itrans)
{
sharp_execute(SHARP_MAP2ALM,spin,&alm[itrans*ncomp],&map[itrans*ncomp],ginfo,ainfo,
SHARP_DP|SHARP_ADD,&ttm2a,&tom2a);
if (t_m2a!=NULL) *t_m2a+=maxTime(ttm2a);
if (op_m2a!=NULL) *op_m2a+=totalops(tom2a);
if (t_m2a!=nullptr) *t_m2a+=maxTime(ttm2a);
if (op_m2a!=nullptr) *op_m2a+=totalops(tom2a);
}
get_errors(alm.data(), nalms, ntrans*ncomp, sqsum, err_abs, err_rel);
}
......@@ -428,8 +428,8 @@ static void check_accuracy (sharp_geom_info &ginfo, sharp_alm_info &ainfo,
{
int ncomp = (spin==0) ? 1 : 2;
vector<double> err_abs, err_rel;
do_sht (ginfo, ainfo, spin, err_abs, err_rel, NULL, NULL,
NULL, NULL, 1);
do_sht (ginfo, ainfo, spin, err_abs, err_rel, nullptr, nullptr,
nullptr, nullptr, 1);
for (int i=0; i<ncomp; ++i)