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

update to new flags in a few missing places; reformatting

Browse source
This commit is contained in:
Martin Reinecke 2012-11-06 19:10:34 +01:00
parent dbe11e8484
commit d377d60d4f
8 changed files with 64 additions and 50 deletions
Download patch file Download diff file Expand all files Collapse all files

35
libsharp/sharp.c
View file

@ -274,14 +274,18 @@ static void ringhelper_phase2ring (ringhelper *self,
}
#endif
real_plan_backward_c (self->plan, (double *)(self->work));
if (flags & SHARP_ALM2MAP_USE_WEIGHTS) {
if (flags & SHARP_ALM2MAP_USE_WEIGHTS)
{
if (fde==DOUBLE)
for (int m=0; m<nph; ++m)
((double *)data)[m*stride+info->ofs] += creal(self->work[m]) * info->weight;
((double *)data)[m*stride+info->ofs]+=creal(self->work[m])*info->weight;
else
for (int m=0; m<nph; ++m)
((float *)data)[m*stride+info->ofs] += (float)creal(self->work[m]) * info->weight;
} else {
((float *)data)[m*stride+info->ofs] +=
(float)(creal(self->work[m])*info->weight);
}
else
{
if (fde==DOUBLE)
for (int m=0; m<nph; ++m)
((double *)data)[m*stride+info->ofs] += creal(self->work[m]);
@ -303,21 +307,24 @@ static void ringhelper_ring2phase (ringhelper *self,
#endif
ringhelper_update (self, nph, mmax, -info->phi0);
if (flags & SHARP_MAP2ALM_IGNORE_WEIGHTS) {
if (flags & SHARP_MAP2ALM_IGNORE_WEIGHTS)
{
if (fde==DOUBLE)
for (int m=0; m<nph; ++m)
self->work[m] = ((double *)data)[info->ofs+m*info->stride];
else
for (int m=0; m<nph; ++m)
self->work[m] = ((float *)data)[info->ofs+m*info->stride];
} else {
}
else
{
if (fde==DOUBLE)
for (int m=0; m<nph; ++m)
self->work[m] = ((double *)data)[info->ofs+m*info->stride]*info->weight;
else
for (int m=0; m<nph; ++m)
self->work[m] = ((float *)data)[info->ofs+m*info->stride]*info->weight;
}
}
real_plan_forward_c (self->plan, (double *)self->work);
@ -336,18 +343,18 @@ static void ringhelper_pair2phase (ringhelper *self, int mmax,
const sharp_ringpair *pair, const void *data, dcmplx *phase1, dcmplx *phase2,
int pstride, sharp_fde fde, int flags)
{
ringhelper_ring2phase (self, &(pair->r1), data, mmax, phase1, pstride, fde, flags);
ringhelper_ring2phase (self,&(pair->r1),data,mmax,phase1,pstride,fde,flags);
if (pair->r2.nph>0)
ringhelper_ring2phase (self, &(pair->r2), data, mmax, phase2, pstride, fde, flags);
ringhelper_ring2phase (self,&(pair->r2),data,mmax,phase2,pstride,fde,flags);
}
static void ringhelper_phase2pair (ringhelper *self, int mmax,
const dcmplx *phase1, const dcmplx *phase2, int pstride,
const sharp_ringpair *pair, void *data, sharp_fde fde, int flags)
{
ringhelper_phase2ring (self, &(pair->r1), data, mmax, phase1, pstride, fde, flags);
ringhelper_phase2ring (self,&(pair->r1),data,mmax,phase1,pstride,fde,flags);
if (pair->r2.nph>0)
ringhelper_phase2ring (self, &(pair->r2), data, mmax, phase2, pstride, fde, flags);
ringhelper_phase2ring (self,&(pair->r2),data,mmax,phase2,pstride,fde,flags);
}
static void fill_map (const sharp_geom_info *ginfo, void *map, double value,
@ -419,7 +426,8 @@ static void map2phase (sharp_job *job, int mmax, int llim, int ulim)
int dim2 = pstride*(ith-llim)*(mmax+1);
for (int i=0; i<job->ntrans*job->nmaps; ++i)
ringhelper_pair2phase(&helper,mmax,&job->ginfo->pair[ith], job->map[i],
&job->phase[dim2+2*i], &job->phase[dim2+2*i+1], pstride, job->fde, job->flags);
&job->phase[dim2+2*i], &job->phase[dim2+2*i+1], pstride, job->fde,
job->flags);
}
ringhelper_destroy(&helper);
} /* end of parallel region */
@ -613,7 +621,8 @@ static void sharp_build_job_common (sharp_job *job, sharp_jobtype type,
{
UTIL_ASSERT((ntrans>0)&&(ntrans<=SHARP_MAXTRANS),
"bad number of simultaneous transforms");
UTIL_ASSERT((flags>0),"passed -1 for old 'dp' argument which is now replaced by 'flags', please pass SHARP_DP");
UTIL_ASSERT((flags>0), "passed -1 for old 'dp' argument which is now\n"
"replaced by 'flags', please pass SHARP_DP");
if (type==SHARP_ALM2MAP_DERIV1) spin=1;
UTIL_ASSERT((spin>=0)&&(spin<=30), "bad spin");
job->type = type;

24
libsharp/sharp_acctest.c
View file

@ -121,8 +121,8 @@ static void check_sign_scale(void)
for (int j=0; j<nalms; ++j)
alm[i][j]=1.+_Complex_I;
sharp_execute(SHARP_ALM2MAP,0,0,&alm[0],&map[0],tinfo,alms,ntrans,SHARP_DP,0,NULL,
NULL);
sharp_execute(SHARP_ALM2MAP,0,0,&alm[0],&map[0],tinfo,alms,ntrans,SHARP_DP,
0,NULL,NULL);
for (int it=0; it<ntrans; ++it)
{
UTIL_ASSERT(FAPPROX(map[it][0 ], 3.588246976618616912e+00,1e-12),
@ -132,8 +132,8 @@ static void check_sign_scale(void)
UTIL_ASSERT(FAPPROX(map[it][npix-1],-1.234675107554816442e+01,1e-12),
"error");
}
sharp_execute(SHARP_ALM2MAP,1,0,&alm[0],&map[0],tinfo,alms,ntrans,SHARP_DP,0,NULL,
NULL);
sharp_execute(SHARP_ALM2MAP,1,0,&alm[0],&map[0],tinfo,alms,ntrans,SHARP_DP,
0,NULL,NULL);
for (int it=0; it<ntrans; ++it)
{
UTIL_ASSERT(FAPPROX(map[2*it ][0 ], 2.750897760535633285e+00,1e-12),
@ -150,8 +150,8 @@ static void check_sign_scale(void)
"error");
}
sharp_execute(SHARP_ALM2MAP,2,0,&alm[0],&map[0],tinfo,alms,ntrans,SHARP_DP,0,NULL,
NULL);
sharp_execute(SHARP_ALM2MAP,2,0,&alm[0],&map[0],tinfo,alms,ntrans,SHARP_DP,
0,NULL,NULL);
for (int it=0; it<ntrans; ++it)
{
UTIL_ASSERT(FAPPROX(map[2*it ][0 ],-1.398186224727334448e+00,1e-12),
@ -168,8 +168,8 @@ static void check_sign_scale(void)
"error");
}
sharp_execute(SHARP_ALM2MAP_DERIV1,1,0,&alm[0],&map[0],tinfo,alms,ntrans,SHARP_DP,
0,NULL,NULL);
sharp_execute(SHARP_ALM2MAP_DERIV1,1,0,&alm[0],&map[0],tinfo,alms,ntrans,
SHARP_DP,0,NULL,NULL);
for (int it=0; it<ntrans; ++it)
{
UTIL_ASSERT(FAPPROX(map[2*it ][0 ],-6.859393905369091105e-01,1e-11),
@ -215,10 +215,10 @@ static void check_accuracy (sharp_geom_info *tinfo, ptrdiff_t lmax,
dcmplx **alm2;
ALLOC2D(alm2,dcmplx,ncomp,nalms);
sharp_execute(SHARP_ALM2MAP,spin,0,&alm[0],&map[0],tinfo,alms,ntrans,SHARP_DP,nv,
NULL,NULL);
sharp_execute(SHARP_MAP2ALM,spin,0,&alm2[0],&map[0],tinfo,alms,ntrans,SHARP_DP,nv,
NULL,NULL);
sharp_execute(SHARP_ALM2MAP,spin,0,&alm[0],&map[0],tinfo,alms,ntrans,
SHARP_DP,nv,NULL,NULL);
sharp_execute(SHARP_MAP2ALM,spin,0,&alm2[0],&map[0],tinfo,alms,ntrans,
SHARP_DP,nv,NULL,NULL);
measure_errors(alm,alm2,nalms,ncomp);
DEALLOC2D(map);

4
libsharp/sharp_bench.c
View file

@ -73,8 +73,8 @@ static void bench_sht (int spin, int nv, sharp_jobtype type,
{
double jtime;
unsigned long long jopcnt;
sharp_execute(type,spin,0,&alm[0],&map[0],tinfo,alms,ntrans,1,nv,&jtime,
&jopcnt);
sharp_execute(type,spin,0,&alm[0],&map[0],tinfo,alms,ntrans,SHARP_DP,nv,
&jtime,&jopcnt);
if (jopcnt<*opcnt) *opcnt=jopcnt;
if (jtime<*time) *time=jtime;

2
libsharp/sharp_bench2.c
View file

@ -189,7 +189,7 @@ int main(int argc, char **argv)
double ltime;
unsigned long long lopcnt;
sharp_execute_mpi(MPI_COMM_WORLD,jtype,spin,0,&alm[0],&map[0],
tinfo,alms,ntrans,1,0,&ltime,&lopcnt);
tinfo,alms,ntrans,SHARP_DP,0,&ltime,&lopcnt);
ltime=maxTime(ltime);
if (ltime<time) { time=ltime; opcnt=totalops(lopcnt); }

4
libsharp/sharp_cxx.h
View file

@ -78,10 +78,10 @@ class sharp_base
template<typename T> struct cxxjobhelper__ {};
template<> struct cxxjobhelper__<double>
{ enum {val=1}; };
{ enum {val=SHARP_DP}; };
template<> struct cxxjobhelper__<float>
{ enum {val=0}; };
{ enum {val=SHARP_SP}; };
template<typename T> class sharp_cxxjob: public sharp_base

23
libsharp/sharp_lowlevel.h
View file

@ -150,16 +150,21 @@ typedef enum { SHARP_MAP2ALM, /*!< analysis */
typedef enum { SHARP_SP = 0, /*!< map and alm is in single precision */
SHARP_DP = 1 << 1, /*!< map and alm is in double precision */
SHARP_ALM2MAP_USE_WEIGHTS = 1 << 2, /*!< apply ring weights for alm2map */
SHARP_MAP2ALM_IGNORE_WEIGHTS = 1 << 3, /*!< do not use ring weights for map2alm */
SHARP_ALM2MAP_USE_WEIGHTS = 1 << 2,
/*!< apply ring weights for alm2map */
SHARP_MAP2ALM_IGNORE_WEIGHTS = 1 << 3,
/*!< do not use ring weights for map2alm */
/* convenience flag combinations (stable API even if the default changes) */
SHARP_USE_WEIGHTS = SHARP_ALM2MAP_USE_WEIGHTS, /*!< use ring weights for both map2alm and alm2map */
SHARP_IGNORE_WEIGHTS = SHARP_MAP2ALM_IGNORE_WEIGHTS /*!< do not use ring weights for either map2alm or map2alm */
/* convenience flag combinations
(stable API even if the default changes) */
SHARP_USE_WEIGHTS = SHARP_ALM2MAP_USE_WEIGHTS,
/*!< use ring weights for both map2alm and alm2map */
SHARP_IGNORE_WEIGHTS = SHARP_MAP2ALM_IGNORE_WEIGHTS
/*!< do not use ring weights for either map2alm or map2alm */
} sharp_jobflags;
/*! Performs a libsharp SHT job. The interface deliberately does not use
the C99 "complex" data type, in order to be callable from C.
the C99 "complex" data type, in order to be callable from C89 and C++.
\param type the type of SHT
\param spin the spin of the quantities to be transformed
\param add_output if 0, the output arrays will be overwritten,
@ -180,9 +185,9 @@ typedef enum { SHARP_SP = 0, /*!< map and alm is in single precision */
\a alm arrays. All \c m values from 0 to some \c mmax<=lmax must be present
exactly once.
\param ntrans the number of simultaneous SHTs
\param flags See sharp_jobflags. In particular, if SHARP_SP is set, the \a alm is
expected to have the type "complex float **"
and \a map is expected to have the type "float **"; if SHARP_DP is set, the expected
\param flags See sharp_jobflags. In particular, if SHARP_SP is set, then
\a alm is expected to have the type "complex float **" and \a map is
expected to have the type "float **"; if SHARP_DP is set, the expected
types are "complex double **" and "double **", respectively.
\param nv Internally used SHT parameter. Set to 0 unless you know what you are
doing.

6
libsharp/sharp_mpi.h
View file

@ -64,9 +64,9 @@ extern "C" {
exactly once in the union of all \a alm_info objects over the participating
MPI tasks.
\param ntrans the number of simultaneous SHTs
\param flags See sharp_jobflags. In particular, if SHARP_SP is set, the \a alm is
expected to have the type "complex float **"
and \a map is expected to have the type "float **"; if SHARP_DP is set, the expected
\param flags See sharp_jobflags. In particular, if SHARP_SP is set, then
\a alm is expected to have the type "complex float **" and \a map is
expected to have the type "float **"; if SHARP_DP is set, the expected
types are "complex double **" and "double **", respectively.
\param nv Internally used SHT parameter. Set to 0 unless you know what you are
doing.

16
libsharp/sharp_test.c
View file

@ -119,8 +119,8 @@ static void map2alm_iter (sharp_geom_info *tinfo, double **map,
double time;
unsigned long long opcnt;
sharp_execute(SHARP_MAP2ALM,spin,0,&alm[0],&map[0],tinfo,alms,ntrans,SHARP_DP,0,
&time,&opcnt);
sharp_execute(SHARP_MAP2ALM,spin,0,&alm[0],&map[0],tinfo,alms,ntrans,
SHARP_DP,0,&time,&opcnt);
printf("wall time for map2alm: %fs\n",time);
printf("Performance: %fGFLOPs/s\n",1e-9*opcnt/time);
measure_errors(alm_orig,alm,nalms,ncomp);
@ -130,16 +130,16 @@ static void map2alm_iter (sharp_geom_info *tinfo, double **map,
double **map2;
ALLOC2D(map2,double,ncomp,npix);
printf ("\niteration %i:\n", iter+1);
sharp_execute(SHARP_ALM2MAP,spin,0,&alm[0],&map2[0],tinfo,alms,ntrans,SHARP_DP,0,
&time,&opcnt);
sharp_execute(SHARP_ALM2MAP,spin,0,&alm[0],&map2[0],tinfo,alms,ntrans,
SHARP_DP,0,&time,&opcnt);
printf("wall time for alm2map: %fs\n",time);
printf("Performance: %fGFLOPs/s\n",1e-9*opcnt/time);
for (int i=0; i<ncomp; ++i)
for (ptrdiff_t m=0; m<npix; ++m)
map2[i][m] = map[i][m]-map2[i][m];
sharp_execute(SHARP_MAP2ALM,spin,1,&alm[0],&map2[0],tinfo,alms,ntrans,SHARP_DP,0,
&time,&opcnt);
sharp_execute(SHARP_MAP2ALM,spin,1,&alm[0],&map2[0],tinfo,alms,ntrans,
SHARP_DP,0,&time,&opcnt);
printf("wall time for map2alm: %fs\n",time);
printf("Performance: %fGFLOPs/s\n",1e-9*opcnt/time);
DEALLOC2D(map2);
@ -173,8 +173,8 @@ static void check_accuracy (sharp_geom_info *tinfo, ptrdiff_t lmax,
double time;
unsigned long long opcnt;
printf ("\niteration 0:\n");
sharp_execute(SHARP_ALM2MAP,spin,0,&alm[0],&map[0],tinfo,alms,ntrans,SHARP_DP,0,
&time,&opcnt);
sharp_execute(SHARP_ALM2MAP,spin,0,&alm[0],&map[0],tinfo,alms,ntrans,
SHARP_DP,0,&time,&opcnt);
printf("wall time for alm2map: %fs\n",time);
printf("Performance: %fGFLOPs/s\n",1e-9*opcnt/time);