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

Merge branch 'dagss' of git://github.com/dagss/libsharp

Browse source 
add the option to specify NULL as a ring weight pointer
add flag to allow weight application also for alm2map
Conflicts:
	libsharp/sharp.c
This commit is contained in:
Dag Sverre Seljebotn 2012-11-06 18:58:29 +01:00 • committed by Martin Reinecke
commit dbe11e8484
8 changed files with 87 additions and 50 deletions
Download patch file Download diff file Expand all files Collapse all files

74
libsharp/sharp.c
View file

@ -177,7 +177,7 @@ void sharp_make_geom_info (int nrings, const int *nph, const ptrdiff_t *ofs,
infos[m].theta = theta[m];
infos[m].cth = cos(theta[m]);
infos[m].sth = sin(theta[m]);
infos[m].weight = weight[m];
infos[m].weight = (weight != NULL) ? weight[m] : 1.;
infos[m].phi0 = phi0[m];
infos[m].ofs = ofs[m];
infos[m].stride = stride[m];
@ -234,7 +234,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;
@ -274,17 +274,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
@ -294,12 +303,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);
@ -316,20 +334,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,
@ -401,7 +419,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 */
@ -499,7 +517,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 */
@ -591,10 +609,11 @@ 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)&&(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");
if (type==SHARP_ALM2MAP_DERIV1) spin=1;
UTIL_ASSERT((spin>=0)&&(spin<=30), "bad spin");
job->type = type;
@ -611,16 +630,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

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

23
libsharp/sharp_lowlevel.h
View file

@ -123,7 +123,9 @@ void sharp_destroy_alm_info (sharp_alm_info *info);
\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 weight the pixel weight to be used for the ring
\param weight the pixel weight to be used for the ring. Pass NULL to use
1.0 as weight for all rings. By default weights are used for alm2map
but not map2alm, but the execution flags can override this.
\param geom_info will hold a pointer to the newly created data structure
*/
void sharp_make_geom_info (int nrings, const int *nph, const ptrdiff_t *ofs,
@ -144,6 +146,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 +180,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 +192,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);
void sharp_set_chunksize_min(int new_chunksize_min);
void sharp_set_nchunks_max(int new_nchunks_max);

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);