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

replace dp with flags and introduce SHARP_ALWAYS_WEIGHTS, SHARP_NEVER_WEIGHTS

Browse source
This commit is contained in:
Dag Sverre Seljebotn 2012-11-06 13:33:44 +01:00
parent bea64bd65f
commit 93a01f769c
8 changed files with 83 additions and 48 deletions
Download patch file Download diff file Expand all files Collapse all files

72
libsharp/sharp.c
View file

@ -233,7 +233,7 @@ static int sharp_get_mmax (int *mval, int nm)
static void ringhelper_phase2ring (ringhelper *self,
const sharp_ringinfo *info, void *data, int mmax, const dcmplx *phase,
int pstride, sharp_fde fde)
int pstride, sharp_fde fde, int flags)
{
int nph = info->nph;
int stride = info->stride;
@ -273,17 +273,26 @@ static void ringhelper_phase2ring (ringhelper *self,
}
#endif
real_plan_backward_c (self->plan, (double *)(self->work));
if (fde==DOUBLE)
for (int m=0; m<nph; ++m)
((double *)data)[m*stride+info->ofs] += creal(self->work[m]);
else
for (int m=0; m<nph; ++m)
((float *)data)[m*stride+info->ofs] += (float)creal(self->work[m]);
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;
else
for (int m=0; m<nph; ++m)
((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]);
else
for (int m=0; m<nph; ++m)
((float *)data)[m*stride+info->ofs] += (float)creal(self->work[m]);
}
}
static void ringhelper_ring2phase (ringhelper *self,
const sharp_ringinfo *info, const void *data, int mmax, dcmplx *phase,
int pstride, sharp_fde fde)
int pstride, sharp_fde fde, int flags)
{
int nph = info->nph;
#if 1
@ -293,12 +302,21 @@ static void ringhelper_ring2phase (ringhelper *self,
#endif
ringhelper_update (self, nph, mmax, -info->phi0);
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;
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 {
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);
@ -315,20 +333,20 @@ static void ringhelper_ring2phase (ringhelper *self,
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 pstride, sharp_fde fde, int flags)
{
ringhelper_ring2phase (self, &(pair->r1), data, mmax, phase1, pstride, fde);
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);
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)
const sharp_ringpair *pair, void *data, sharp_fde fde, int flags)
{
ringhelper_phase2ring (self, &(pair->r1), data, mmax, phase1, pstride, fde);
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);
ringhelper_phase2ring (self, &(pair->r2), data, mmax, phase2, pstride, fde, flags);
}
static void fill_map (const sharp_geom_info *ginfo, void *map, double value,
@ -400,7 +418,7 @@ 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->phase[dim2+2*i], &job->phase[dim2+2*i+1], pstride, job->fde, job->flags);
}
ringhelper_destroy(&helper);
} /* end of parallel region */
@ -463,7 +481,7 @@ static void phase2map (sharp_job *job, int mmax, int llim, int ulim)
for (int i=0; i<job->ntrans*job->nmaps; ++i)
ringhelper_phase2pair(&helper,mmax,&job->phase[dim2+2*i],
&job->phase[dim2+2*i+1],pstride,&job->ginfo->pair[ith],job->map[i],
job->fde);
job->fde, job->flags);
}
ringhelper_destroy(&helper);
} /* end of parallel region */
@ -555,9 +573,10 @@ static void sharp_execute_job (sharp_job *job)
static void sharp_build_job_common (sharp_job *job, sharp_jobtype type,
int spin, int add_output, void *alm, void *map,
const sharp_geom_info *geom_info, const sharp_alm_info *alm_info, int ntrans,
int dp, int nv)
int flags, int nv)
{
UTIL_ASSERT((ntrans>0),"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");
if (type==SHARP_ALM2MAP_DERIV1) spin=1;
UTIL_ASSERT((spin>=0)&&(spin<=30), "bad spin");
job->type = type;
@ -574,16 +593,17 @@ static void sharp_build_job_common (sharp_job *job, sharp_jobtype type,
job->ntrans = ntrans;
job->alm=alm;
job->map=map;
job->fde=dp ? DOUBLE : FLOAT;
job->flags = flags;
job->fde=(flags & SHARP_DP) ? DOUBLE : FLOAT;
}
void sharp_execute (sharp_jobtype type, int spin, int add_output, void *alm,
void *map, const sharp_geom_info *geom_info, const sharp_alm_info *alm_info,
int ntrans, int dp, int nv, double *time, unsigned long long *opcnt)
int ntrans, int flags, int nv, double *time, unsigned long long *opcnt)
{
sharp_job job;
sharp_build_job_common (&job, type, spin, add_output, alm, map, geom_info,
alm_info, ntrans, dp, nv);
alm_info, ntrans, flags, nv);
sharp_execute_job (&job);
if (time!=NULL) *time = job.time;

12
libsharp/sharp_acctest.c
View file

@ -121,7 +121,7 @@ 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,1,0,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)
{
@ -132,7 +132,7 @@ 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,1,0,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)
{
@ -150,7 +150,7 @@ static void check_sign_scale(void)
"error");
}
sharp_execute(SHARP_ALM2MAP,2,0,&alm[0],&map[0],tinfo,alms,ntrans,1,0,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)
{
@ -168,7 +168,7 @@ static void check_sign_scale(void)
"error");
}
sharp_execute(SHARP_ALM2MAP_DERIV1,1,0,&alm[0],&map[0],tinfo,alms,ntrans,1,
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)
{
@ -215,9 +215,9 @@ 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,1,nv,
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,1,nv,
sharp_execute(SHARP_MAP2ALM,spin,0,&alm2[0],&map[0],tinfo,alms,ntrans,SHARP_DP,nv,
NULL,NULL);
measure_errors(alm,alm2,nalms,ncomp);

1
libsharp/sharp_internal.h
View file

@ -46,6 +46,7 @@ typedef struct
int spin;
int add_output;
int nmaps, nalm;
int flags;
sharp_fde fde;
void **map;
void **alm;

19
libsharp/sharp_lowlevel.h
View file

@ -144,6 +144,18 @@ typedef enum { SHARP_MAP2ALM, /*!< analysis */
SHARP_ALM2MAP_DERIV1 /*!< synthesis of first derivatives */
} sharp_jobtype;
/*! Job flags */
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 */
/* 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.
\param type the type of SHT
@ -166,8 +178,9 @@ typedef enum { SHARP_MAP2ALM, /*!< analysis */
\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 dp if 0, the \a alm is expected to have the type "complex float **"
and \a map is expected to have the type "float **"; otherwise the expected
\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
types are "complex double **" and "double **", respectively.
\param nv Internally used SHT parameter. Set to 0 unless you know what you are
doing.
@ -177,7 +190,7 @@ typedef enum { SHARP_MAP2ALM, /*!< analysis */
operation count for this SHT will be written here. */
void sharp_execute (sharp_jobtype type, int spin, int add_output, void *alm,
void *map, const sharp_geom_info *geom_info, const sharp_alm_info *alm_info,
int ntrans, int dp, int nv, double *time, unsigned long long *opcnt);
int ntrans, int flags, int nv, double *time, unsigned long long *opcnt);
/*! \} */

4
libsharp/sharp_mpi.c
View file

@ -285,12 +285,12 @@ static void sharp_execute_job_mpi (sharp_job *job, MPI_Comm comm)
void sharp_execute_mpi (MPI_Comm comm, sharp_jobtype type, int spin,
int add_output, void *alm, void *map, const sharp_geom_info *geom_info,
const sharp_alm_info *alm_info, int ntrans, int dp, int nv, double *time,
const sharp_alm_info *alm_info, int ntrans, int flags, int nv, double *time,
unsigned long long *opcnt)
{
sharp_job job;
sharp_build_job_common (&job, type, spin, add_output, alm, map, geom_info,
alm_info, ntrans, dp, nv);
alm_info, ntrans, flags, nv);
sharp_execute_job_mpi (&job, comm);
if (time!=NULL) *time = job.time;

7
libsharp/sharp_mpi.h
View file

@ -64,8 +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 dp if 0, the \a alm is expected to have the type "complex float **"
and \a map is expected to have the type "float **"; otherwise the expected
\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
types are "complex double **" and "double **", respectively.
\param nv Internally used SHT parameter. Set to 0 unless you know what you are
doing.
@ -75,7 +76,7 @@ extern "C" {
operation count for this SHT will be written here. */
void sharp_execute_mpi (MPI_Comm comm, sharp_jobtype type, int spin,
int add_output, void *alm, void *map, const sharp_geom_info *geom_info,
const sharp_alm_info *alm_info, int ntrans, int dp, int nv, double *time,
const sharp_alm_info *alm_info, int ntrans, int flags, int nv, double *time,
unsigned long long *opcnt);
#ifdef __cplusplus

8
libsharp/sharp_test.c
View file

@ -119,7 +119,7 @@ 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,1,0,
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);
@ -130,7 +130,7 @@ 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,1,0,
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);
@ -138,7 +138,7 @@ static void map2alm_iter (sharp_geom_info *tinfo, double **map,
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,1,0,
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);
@ -173,7 +173,7 @@ 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,1,0,
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);

8
libsharp/sharp_test_mpi.c
View file

@ -201,7 +201,7 @@ static void map2alm_iter (sharp_geom_info *tinfo, double **map,
unsigned long long jopcnt;
sharp_execute_mpi(MPI_COMM_WORLD,SHARP_MAP2ALM,spin,0,&alm[0],&map[0],
tinfo,alms,ntrans,1,0,&jtime,&jopcnt);
tinfo,alms,ntrans,SHARP_DP,0,&jtime,&jopcnt);
unsigned long long opcnt=totalops(jopcnt);
double timer=maxTime(jtime);
if (mytask==0) printf("wall time for map2alm: %fs\n",timer);
@ -214,7 +214,7 @@ static void map2alm_iter (sharp_geom_info *tinfo, double **map,
ALLOC2D(map2,double,ncomp,npix);
if (mytask==0) printf ("\niteration %i:\n", iter+1);
sharp_execute_mpi(MPI_COMM_WORLD,SHARP_ALM2MAP,spin,0,&alm[0],&map2[0],
tinfo,alms,ntrans,1,0,&jtime,&jopcnt);
tinfo,alms,ntrans,SHARP_DP,0,&jtime,&jopcnt);
opcnt=totalops(jopcnt);
timer=maxTime(jtime);
if (mytask==0) printf("wall time for alm2map: %fs\n",timer);
@ -224,7 +224,7 @@ static void map2alm_iter (sharp_geom_info *tinfo, double **map,
map2[i][m] = map[i][m]-map2[i][m];
sharp_execute_mpi(MPI_COMM_WORLD,SHARP_MAP2ALM,spin,1,&alm[0],&map2[0],
tinfo,alms,ntrans,1,0,&jtime,&jopcnt);
tinfo,alms,ntrans,SHARP_DP,0,&jtime,&jopcnt);
opcnt=totalops(jopcnt);
timer=maxTime(jtime);
if (mytask==0) printf("wall time for map2alm: %fs\n",wallTime()-timer);
@ -264,7 +264,7 @@ static void check_accuracy (sharp_geom_info *tinfo, ptrdiff_t lmax,
if (mytask==0) printf ("\niteration 0:\n");
sharp_execute_mpi(MPI_COMM_WORLD,SHARP_ALM2MAP,spin,0,&alm[0],&map[0],
tinfo,alms,ntrans,1,0,&jtime,&jopcnt);
tinfo,alms,ntrans,SHARP_DP,0,&jtime,&jopcnt);
unsigned long long opcnt=totalops(jopcnt);
double timer=maxTime(jtime);
if (mytask==0) printf("wall time for alm2map: %fs\n",timer);