Skip to content

GitLab

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

cleanup

parent 8d10272f
Hide whitespace changes
Inline Side-by-side
libsharp2/sharp.cc
View file @ b6b06fd8
......@@ -130,11 +130,11 @@ static int ringpair_compare (const void *xa, const void *xb)
void sharp_make_general_alm_info (int lmax, int nm, int stride, const int *mval,
const ptrdiff_t *mstart, int flags, sharp_alm_info **alm_info)
{
sharp_alm_info *info = RALLOC(sharp_alm_info,1);
sharp_alm_info *info = new sharp_alm_info;
info->lmax = lmax;
info->nm = nm;
info->mval = RALLOC(int,nm);
info->mvstart = RALLOC(ptrdiff_t,nm);
info->mval.resize(nm);
info->mvstart.resize(nm);
info->stride = stride;
info->flags = flags;
for (int mi=0; mi<nm; ++mi)
......@@ -176,21 +176,19 @@ ptrdiff_t sharp_alm_count(const sharp_alm_info *self)
void sharp_destroy_alm_info (sharp_alm_info *info)
{
DEALLOC (info->mval);
DEALLOC (info->mvstart);
DEALLOC (info);
delete info;
}
void 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 **geom_info)
{
sharp_geom_info *info = RALLOC(sharp_geom_info,1);
sharp_geom_info *info = new sharp_geom_info;
vector<sharp_ringinfo> infos(nrings);
int pos=0;
info->pair=RALLOC(sharp_ringpair,nrings);
info->npairs=0;
info->pair.resize(nrings);
int npairs=0;
info->nphmax=0;
*geom_info = info;
......@@ -209,31 +207,32 @@ void sharp_make_geom_info (int nrings, const int *nph, const ptrdiff_t *ofs,
qsort(infos.data(),nrings,sizeof(sharp_ringinfo),ringinfo_compare);
while (pos<nrings)
{
info->pair[info->npairs].r1=infos[pos];
info->pair[npairs].r1=infos[pos];
if ((pos<nrings-1) && FAPPROX(infos[pos].cth,-infos[pos+1].cth,1e-12))
{
if (infos[pos].cth>0) // make sure northern ring is in r1
info->pair[info->npairs].r2=infos[pos+1];
info->pair[npairs].r2=infos[pos+1];
else
{
info->pair[info->npairs].r1=infos[pos+1];
info->pair[info->npairs].r2=infos[pos];
info->pair[npairs].r1=infos[pos+1];
info->pair[npairs].r2=infos[pos];
}
++pos;
}
else
info->pair[info->npairs].r2.nph=-1;
info->pair[npairs].r2.nph=-1;
++pos;
++info->npairs;
++npairs;
}
qsort(info->pair,info->npairs,sizeof(sharp_ringpair),ringpair_compare);
qsort(info->pair.data(),npairs,sizeof(sharp_ringpair),ringpair_compare);
info->pair.resize(npairs);
}
ptrdiff_t sharp_map_size(const sharp_geom_info *info)
{
ptrdiff_t result = 0;
for (int m=0; m<info->npairs; ++m)
for (int m=0; m<info->pair.size(); ++m)
{
result+=info->pair[m].r1.nph;
result+=(info->pair[m].r2.nph>=0) ? (info->pair[m].r2.nph) : 0;
......@@ -243,16 +242,16 @@ ptrdiff_t sharp_map_size(const sharp_geom_info *info)
void sharp_destroy_geom_info (sharp_geom_info *geom_info)
{
DEALLOC (geom_info->pair);
DEALLOC (geom_info);
delete geom_info;
}
/* This currently requires all m values from 0 to nm-1 to be present.
It might be worthwhile to relax this criterion such that holes in the m
distribution are permissible. */
static int sharp_get_mmax (int *mval, int nm)
static int sharp_get_mmax (const vector<int> &mval)
{
//FIXME: if gaps are allowed, we have to search the maximum m in the array
auto nm=mval.size();
vector<int> mcheck(nm,0);
for (int i=0; i<nm; ++i)
{
......@@ -377,7 +376,7 @@ MRUTIL_NOINLINE static void clear_map (const sharp_geom_info *ginfo, void *map,
{
if (flags & SHARP_NO_FFT)
{
for (int j=0;j<ginfo->npairs;++j)
for (int j=0;j<ginfo->pair.size();++j)
{
if (flags&SHARP_DP)
{
......@@ -399,7 +398,7 @@ MRUTIL_NOINLINE static void clear_map (const sharp_geom_info *ginfo, void *map,
{
if (flags&SHARP_DP)
{
for (int j=0;j<ginfo->npairs;++j)
for (int j=0;j<ginfo->pair.size();++j)
{
double *dmap=(double *)map;
if (ginfo->pair[j].r1.stride==1)
......@@ -418,7 +417,7 @@ MRUTIL_NOINLINE static void clear_map (const sharp_geom_info *ginfo, void *map,
}
else
{
for (int j=0;j<ginfo->npairs;++j)
for (int j=0;j<ginfo->pair.size();++j)
{
for (ptrdiff_t i=0;i<ginfo->pair[j].r1.nph;++i)
((float *)map)[ginfo->pair[j].r1.ofs+i*ginfo->pair[j].r1.stride]=0;
......@@ -639,7 +638,7 @@ MRUTIL_NOINLINE static void almtmp2alm (sharp_job *job, int lmax, int mi)
#undef COPY_LOOP
}
MRUTIL_NOINLINE static void ringtmp2ring (sharp_job *job, sharp_ringinfo *ri,
MRUTIL_NOINLINE static void ringtmp2ring (sharp_job *job, const sharp_ringinfo *ri,
const vector<double> &ringtmp, int rstride)
{
if (job->flags & SHARP_DP)
......@@ -671,7 +670,7 @@ MRUTIL_NOINLINE static void ringtmp2ring (sharp_job *job, sharp_ringinfo *ri,
}
}
MRUTIL_NOINLINE static void ring2ringtmp (sharp_job *job, sharp_ringinfo *ri,
MRUTIL_NOINLINE static void ring2ringtmp (sharp_job *job, const sharp_ringinfo *ri,
vector<double> &ringtmp, int rstride)
{
if (job->flags & SHARP_DP)
......@@ -691,7 +690,7 @@ MRUTIL_NOINLINE static void ring2ringtmp (sharp_job *job, sharp_ringinfo *ri,
ringtmp[i*rstride+m+1] = ((float *)(job->map[i]))[ri->ofs+m*ri->stride];
}
static void ring2phase_direct (sharp_job *job, sharp_ringinfo *ri, int mmax,
static void ring2phase_direct (sharp_job *job, const sharp_ringinfo *ri, int mmax,
dcmplx *phase)
{
if (ri->nph<0)
......@@ -713,7 +712,7 @@ static void ring2phase_direct (sharp_job *job, sharp_ringinfo *ri, int mmax,
((fcmplx *)(job->map[i]))[ri->ofs+m*ri->stride]*float(wgt);
}
}
static void phase2ring_direct (sharp_job *job, sharp_ringinfo *ri, int mmax,
static void phase2ring_direct (sharp_job *job, const sharp_ringinfo *ri, int mmax,
dcmplx *phase)
{
if (ri->nph<0) return;
......@@ -821,7 +820,7 @@ MRUTIL_NOINLINE static void sharp_execute_job (sharp_job *job)
mr::timers::SimpleTimer timer;
job->opcnt=0;
int lmax = job->ainfo->lmax,
mmax=sharp_get_mmax(job->ainfo->mval, job->ainfo->nm);
mmax=sharp_get_mmax(job->ainfo->mval);
job->norm_l = (job->type==SHARP_ALM2MAP_DERIV1) ?
sharp_Ylmgen::get_d1norm (lmax) :
......@@ -831,7 +830,7 @@ MRUTIL_NOINLINE static void sharp_execute_job (sharp_job *job)
init_output (job);
int nchunks, chunksize;
get_chunk_info(job->ginfo->npairs,sharp_veclen()*sharp_max_nvec(job->spin),
get_chunk_info(job->ginfo->pair.size(),sharp_veclen()*sharp_max_nvec(job->spin),
&nchunks,&chunksize);
//FIXME: needs to be changed to "nm"
alloc_phase (job,mmax+1,chunksize);
......@@ -840,7 +839,7 @@ MRUTIL_NOINLINE static void sharp_execute_job (sharp_job *job)
/* chunk loop */
for (int chunk=0; chunk<nchunks; ++chunk)
{
int llim=chunk*chunksize, ulim=min(llim+chunksize,job->ginfo->npairs);
int llim=chunk*chunksize, ulim=min<int>(llim+chunksize,job->ginfo->pair.size());
vector<int> ispair(ulim-llim);
vector<int> mlim(ulim-llim);
vector<double> cth(ulim-llim), sth(ulim-llim);
......
libsharp2/sharp.h
View file @ b6b06fd8
......@@ -28,7 +28,8 @@
#ifndef SHARP_SHARP_H
#define SHARP_SHARP_H
#include <stddef.h>
#include <cstddef>
#include <vector>
/*! \internal
Helper type containing information about a single ring. */
......@@ -51,8 +52,8 @@ typedef struct
Type holding all required information about a map geometry. */
typedef struct
{
sharp_ringpair *pair;
int npairs, nphmax;
std::vector<sharp_ringpair> pair;
int nphmax;
} sharp_geom_info;
/*! \defgroup almgroup Helpers for dealing with a_lm */
......@@ -67,12 +68,12 @@ typedef struct
/*! Number of different \a m values in this object */
int nm;
/*! Array with \a nm entries containing the individual m values */
int *mval;
std::vector<int> mval;
/*! Combination of flags from sharp_almflags */
int flags;
/*! Array with \a nm entries containing the (hypothetical) indices of
the coefficients with quantum numbers 0,\a mval[i] */
ptrdiff_t *mvstart;
std::vector<ptrdiff_t> mvstart;
/*! Stride between a_lm and a_(l+1),m */
ptrdiff_t stride;
} sharp_alm_info;
......
libsharp2/sharp_almhelpers.cc
View file @ b6b06fd8
......@@ -31,11 +31,11 @@
void sharp_make_triangular_alm_info (int lmax, int mmax, int stride,
sharp_alm_info **alm_info)
{
sharp_alm_info *info = RALLOC(sharp_alm_info,1);
sharp_alm_info *info = new sharp_alm_info;
info->lmax = lmax;
info->nm = mmax+1;
info->mval = RALLOC(int,mmax+1);
info->mvstart = RALLOC(ptrdiff_t,mmax+1);
info->mval.resize(mmax+1);
info->mvstart.resize(mmax+1);
info->stride = stride;
info->flags = 0;
ptrdiff_t tval = 2*lmax+1;
......@@ -50,11 +50,11 @@ void sharp_make_triangular_alm_info (int lmax, int mmax, int stride,
void sharp_make_rectangular_alm_info (int lmax, int mmax, int stride,
sharp_alm_info **alm_info)
{
sharp_alm_info *info = RALLOC(sharp_alm_info,1);
sharp_alm_info *info = new sharp_alm_info;
info->lmax = lmax;
info->nm = mmax+1;
info->mval = RALLOC(int,mmax+1);
info->mvstart = RALLOC(ptrdiff_t,mmax+1);
info->mval.resize(mmax+1);
info->mvstart.resize(mmax+1);
info->stride = stride;
info->flags = 0;
for (ptrdiff_t m=0; m<=mmax; ++m)
......@@ -70,11 +70,11 @@ void sharp_make_mmajor_real_packed_alm_info (int lmax, int stride,
{
ptrdiff_t idx;
int f;
sharp_alm_info *info = RALLOC(sharp_alm_info,1);
sharp_alm_info *info = new sharp_alm_info;
info->lmax = lmax;
info->nm = nm;
info->mval = RALLOC(int,nm);
info->mvstart = RALLOC(ptrdiff_t,nm);
info->mval.resize(nm);
info->mvstart.resize(nm);
info->stride = stride;
info->flags = SHARP_PACKED | SHARP_REAL_HARMONICS;
idx = 0; /* tracks the number of 'consumed' elements so far; need to correct by m */
......
test/sharp2_testsuite.cc
View file @ b6b06fd8
......@@ -186,7 +186,7 @@ static void reduce_geom_info(sharp_geom_info *ginfo)
{
int npairsnew=0;
ptrdiff_t ofs = 0;
for (int i=mytask; i<ginfo->npairs; i+=ntasks,++npairsnew)
for (int i=mytask; i<ginfo->pair.size(); i+=ntasks,++npairsnew)
{
ginfo->pair[npairsnew]=ginfo->pair[i];
ginfo->pair[npairsnew].r1.ofs=ofs;
......@@ -194,7 +194,7 @@ static void reduce_geom_info(sharp_geom_info *ginfo)
ginfo->pair[npairsnew].r2.ofs=ofs;
if (ginfo->pair[npairsnew].r2.nph>0) ofs+=ginfo->pair[npairsnew].r2.nph;
}
ginfo->npairs=npairsnew;
ginfo.pair.resize(npairsnew);
}
#endif
......@@ -209,7 +209,7 @@ 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 (int i=0; i<ginfo->npairs; ++i)
for (int 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;
......@@ -217,9 +217,9 @@ static ptrdiff_t get_npix(const sharp_geom_info *ginfo)
return res;
}
static double *get_sqsum_and_invert (dcmplx **alm, ptrdiff_t nalms, int ncomp)
static vector<double> get_sqsum_and_invert (dcmplx **alm, ptrdiff_t nalms, int ncomp)
{
double *sqsum=RALLOC(double,ncomp);
vector<double> sqsum(ncomp);
for (int i=0; i<ncomp; ++i)
{
sqsum[i]=0;
......@@ -232,16 +232,16 @@ static double *get_sqsum_and_invert (dcmplx **alm, ptrdiff_t nalms, int ncomp)
return sqsum;
}
static void get_errors (dcmplx **alm, ptrdiff_t nalms, int ncomp, double *sqsum,
double **err_abs, double **err_rel)
static void get_errors (dcmplx **alm, ptrdiff_t nalms, int ncomp, const vector<double> &sqsum,
vector<double> &err_abs, vector<double> &err_rel)
{
long nalms_tot=nalms;
#ifdef USE_MPI
MPI_Allreduce(&nalms,&nalms_tot,1,MPI_LONG,MPI_SUM,MPI_COMM_WORLD);
#endif
*err_abs=RALLOC(double,ncomp);
*err_rel=RALLOC(double,ncomp);
err_abs.resize(ncomp);
err_rel.resize(ncomp);
for (int i=0; i<ncomp; ++i)
{
double sum=0, maxdiff=0, sumtot, sqsumtot, maxdifftot;
......@@ -264,8 +264,8 @@ static void get_errors (dcmplx **alm, ptrdiff_t nalms, int ncomp, double *sqsum,
#endif
sumtot=sqrt(sumtot/nalms_tot);
sqsumtot=sqrt(sqsumtot/nalms_tot);
(*err_abs)[i]=maxdifftot;
(*err_rel)[i]=sumtot/sqsumtot;
err_abs[i]=maxdifftot;
err_rel[i]=sumtot/sqsumtot;
}
}
......@@ -343,7 +343,7 @@ static void get_infos (const char *gname, int lmax, int *mmax, int *gpar1,
sharp_make_gauss_geom_info (nlat, nlon, 0., 1, nlon, ginfo);
ptrdiff_t npix_o=get_npix(*ginfo);
size_t ofs=0;
for (int i=0; i<(*ginfo)->npairs; ++i)
for (int i=0; i<(*ginfo)->pair.size(); ++i)
{
sharp_ringpair *pair=&((*ginfo)->pair[i]);
int pring=1+2*sharp_get_mlim(lmax,0,pair->r1.sth,pair->r1.cth);
......@@ -387,7 +387,7 @@ static void check_sign_scale(void)
sharp_make_gauss_geom_info (nrings, ppring, 0., 1, ppring, &tinfo);
/* flip theta to emulate the "old" Gaussian grid geometry */
for (int i=0; i<tinfo->npairs; ++i)
for (int i=0; i<tinfo->pair.size(); ++i)
{
const double pi=3.141592653589793238462643383279502884197;
tinfo->pair[i].r1.cth=-tinfo->pair[i].r1.cth;
......@@ -471,7 +471,7 @@ static void check_sign_scale(void)
}
static void do_sht (sharp_geom_info *ginfo, sharp_alm_info *ainfo,
int spin, double **err_abs, double **err_rel,
int spin, vector<double> &err_abs, vector<double> &err_rel,
double *t_a2m, double *t_m2a, unsigned long long *op_a2m,
unsigned long long *op_m2a, size_t ntrans)
{
......@@ -506,7 +506,7 @@ static void do_sht (sharp_geom_info *ginfo, sharp_alm_info *ainfo,
if (t_a2m!=NULL) *t_a2m+=maxTime(tta2m);
if (op_a2m!=NULL) *op_a2m+=totalops(toa2m);
}
double *sqsum=get_sqsum_and_invert(alm,nalms,ntrans*ncomp);
auto sqsum=get_sqsum_and_invert(alm,nalms,ntrans*ncomp);
if (t_m2a!=NULL) *t_m2a=0;
if (op_m2a!=NULL) *op_m2a=0;
for (size_t itrans=0; itrans<ntrans; ++itrans)
......@@ -523,7 +523,6 @@ static void do_sht (sharp_geom_info *ginfo, sharp_alm_info *ainfo,
}
get_errors(alm, nalms, ntrans*ncomp, sqsum, err_abs, err_rel);
DEALLOC(sqsum);
DEALLOC2D(map);
DEALLOC2D(alm);
}
......@@ -532,13 +531,11 @@ static void check_accuracy (sharp_geom_info *ginfo, sharp_alm_info *ainfo,
int spin)
{
int ncomp = (spin==0) ? 1 : 2;
double *err_abs, *err_rel;
do_sht (ginfo, ainfo, spin, &err_abs, &err_rel, NULL, NULL,
vector<double> err_abs, err_rel;
do_sht (ginfo, ainfo, spin, err_abs, err_rel, NULL, NULL,
NULL, NULL, 1);
for (int i=0; i<ncomp; ++i)
MR_assert((err_rel[i]<1e-10) && (err_abs[i]<1e-10),"error");
DEALLOC(err_rel);
DEALLOC(err_abs);
}
static void run(int lmax, int mmax, int nlat, int nlon, int spin)
......@@ -596,7 +593,7 @@ static void sharp_test (int argc, const char **argv)
int ncomp = (spin==0) ? 1 : 2;
double t_a2m=1e30, t_m2a=1e30;
unsigned long long op_a2m, op_m2a;
double *err_abs,*err_rel;
vector<double> err_abs, err_rel;
double t_acc=0;
int nrpt=0;
......@@ -604,7 +601,7 @@ static void sharp_test (int argc, const char **argv)
{
++nrpt;
double ta2m2, tm2a2;
do_sht (ginfo, ainfo, spin, &err_abs, &err_rel, &ta2m2, &tm2a2,
do_sht (ginfo, ainfo, spin, err_abs, err_rel, &ta2m2, &tm2a2,
&op_a2m, &op_m2a, ntrans);
if (ta2m2<t_a2m) t_a2m=ta2m2;
if (tm2a2<t_m2a) t_m2a=tm2a2;
......@@ -614,8 +611,6 @@ static void sharp_test (int argc, const char **argv)
if (mytask==0) printf("Best of %d runs\n",nrpt);
break;
}
DEALLOC(err_abs);
DEALLOC(err_rel);
}
if (mytask==0) printf("wall time for alm2map: %fs\n",t_a2m);
......@@ -659,9 +654,6 @@ static void sharp_test (int argc, const char **argv)
getenv("HOST"),argv[2],spin,sharp_veclen(),nomp,ntasks,lmax,mmax,gpar1,gpar2,
t_a2m,1e-9*op_a2m/t_a2m,t_m2a,1e-9*op_m2a/t_m2a,tmem/(1<<20),
100.*(1.-iosize/tmem),maxerel,maxeabs);
DEALLOC(err_abs);
DEALLOC(err_rel);
}
int main(int argc, const char **argv)
......
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