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

better support for adjoint SHTs

Browse source
This commit is contained in:
Martin Reinecke 2012-12-03 11:08:03 +01:00
parent 76be57342a
commit 0d8d82677d
5 changed files with 48 additions and 27 deletions
Download patch file Download diff file Expand all files Collapse all files

20
libsharp/sharp.c
View file

@ -162,7 +162,7 @@ void sharp_destroy_alm_info (sharp_alm_info *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_a2m, const double *wgt_m2a, sharp_geom_info **geom_info)
const double *wgt, sharp_geom_info **geom_info)
{
sharp_geom_info *info = RALLOC(sharp_geom_info,1);
sharp_ringinfo *infos = RALLOC(sharp_ringinfo,nrings);
@ -177,8 +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].w_a2m = (wgt_a2m != NULL) ? wgt_a2m[m] : 1.;
infos[m].w_m2a = (wgt_m2a != NULL) ? wgt_m2a[m] : 1.;
infos[m].weight = (wgt != NULL) ? wgt[m] : 1.;
infos[m].phi0 = phi0[m];
infos[m].ofs = ofs[m];
infos[m].stride = stride[m];
@ -275,13 +274,13 @@ static void ringhelper_phase2ring (ringhelper *self,
}
#endif
real_plan_backward_c (self->plan, (double *)(self->work));
double wgt = (flags&SHARP_USE_WEIGHTS) ? info->weight : 1.;
if (flags&SHARP_DP)
for (int m=0; m<nph; ++m)
((double *)data)[m*stride+info->ofs]+=creal(self->work[m])*info->w_a2m;
((double *)data)[m*stride+info->ofs]+=creal(self->work[m])*wgt;
else
for (int m=0; m<nph; ++m)
((float *)data)[m*stride+info->ofs] +=
(float)(creal(self->work[m])*info->w_a2m);
((float *)data)[m*stride+info->ofs] += (float)(creal(self->work[m])*wgt);
}
static void ringhelper_ring2phase (ringhelper *self,
@ -296,12 +295,13 @@ static void ringhelper_ring2phase (ringhelper *self,
#endif
ringhelper_update (self, nph, mmax, -info->phi0);
double wgt = (flags&SHARP_USE_WEIGHTS) ? info->weight : 1;
if (flags&SHARP_DP)
for (int m=0; m<nph; ++m)
self->work[m] = ((double *)data)[info->ofs+m*info->stride]*info->w_m2a;
self->work[m] = ((double *)data)[info->ofs+m*info->stride]*wgt;
else
for (int m=0; m<nph; ++m)
self->work[m] = ((float *)data)[info->ofs+m*info->stride]*info->w_m2a;
self->work[m] = ((float *)data)[info->ofs+m*info->stride]*wgt;
real_plan_forward_c (self->plan, (double *)self->work);
@ -598,6 +598,10 @@ static void sharp_build_job_common (sharp_job *job, sharp_jobtype type,
UTIL_ASSERT((ntrans>0)&&(ntrans<=SHARP_MAXTRANS),
"bad number of simultaneous transforms");
if (type==SHARP_ALM2MAP_DERIV1) spin=1;
if (type==SHARP_MAP2ALM) flags|=SHARP_USE_WEIGHTS;
if (type==SHARP_Yt) type=SHARP_MAP2ALM;
if (type==SHARP_WY) { type=SHARP_ALM2MAP; flags|=SHARP_USE_WEIGHTS; }
UTIL_ASSERT((spin>=0)&&(spin<=30), "bad spin");
job->type = type;
job->spin = spin;

5
libsharp/sharp_core_inc2.c
View file

@ -804,6 +804,11 @@ static void Z(inner_loop) (sharp_job *job, const int *ispair,
}
break;
}
default:
{
UTIL_FAIL("must not happen");
break;
}
}
}

5
libsharp/sharp_core_inc3.c
View file

@ -794,6 +794,11 @@ static void Y(inner_loop) (sharp_job *job, const int *ispair,
}
break;
}
default:
{
UTIL_FAIL("must not happen");
break;
}
}
}

13
libsharp/sharp_geomhelpers.c
View file

@ -90,7 +90,7 @@ void sharp_make_weighted_healpix_geom_info (int nside, int stride,
weight_[m]=4.*pi/npix*weight[northring-1];
}
sharp_make_geom_info (nrings, nph, ofs, stride_, phi0, theta, NULL, weight_,
sharp_make_geom_info (nrings, nph, ofs, stride_, phi0, theta, weight_,
geom_info);
DEALLOC(theta);
@ -182,7 +182,7 @@ void sharp_make_gauss_geom_info (int nrings, int nphi, double phi0,
weight[m]*=2*pi/nphi;
}
sharp_make_geom_info (nrings, nph, ofs, stride_, phi0_, theta, NULL, weight,
sharp_make_geom_info (nrings, nph, ofs, stride_, phi0_, theta, weight,
geom_info);
DEALLOC(theta);
@ -229,7 +229,7 @@ void sharp_make_ecp_geom_info (int nrings, int ppring, double phi0,
weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(nrings*nph[m]);
}
sharp_make_geom_info (nrings, nph, ofs, stride_, phi0_, theta, NULL, weight,
sharp_make_geom_info (nrings, nph, ofs, stride_, phi0_, theta, weight,
geom_info);
DEALLOC(theta);
@ -278,7 +278,7 @@ void sharp_make_hw_geom_info (int nrings, int ppring, double phi0,
weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(n*nph[m]);
}
sharp_make_geom_info (nrings, nph, ofs, stride_, phi0_, theta, NULL, weight,
sharp_make_geom_info (nrings, nph, ofs, stride_, phi0_, theta, weight,
geom_info);
DEALLOC(theta);
@ -326,7 +326,7 @@ void sharp_make_fejer2_geom_info (int nrings, int ppring, double phi0,
weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(n*nph[m]);
}
sharp_make_geom_info (nrings, nph, ofs, stride_, phi0_, theta, NULL, weight,
sharp_make_geom_info (nrings, nph, ofs, stride_, phi0_, theta, weight,
geom_info);
DEALLOC(theta);
@ -336,6 +336,7 @@ void sharp_make_fejer2_geom_info (int nrings, int ppring, double phi0,
DEALLOC(ofs);
DEALLOC(stride_);
}
void sharp_make_mw_geom_info (int nrings, int ppring, double phi0,
int stride_lon, int stride_lat, sharp_geom_info **geom_info)
{
@ -357,7 +358,7 @@ void sharp_make_mw_geom_info (int nrings, int ppring, double phi0,
stride_[m]=stride_lon;
}
sharp_make_geom_info (nrings, nph, ofs, stride_, phi0_, theta, NULL, NULL,
sharp_make_geom_info (nrings, nph, ofs, stride_, phi0_, theta, NULL,
geom_info);
DEALLOC(theta);

32
libsharp/sharp_lowlevel.h
View file

@ -42,7 +42,7 @@ extern "C" {
Helper type containing information about a single ring. */
typedef struct
{
double theta, phi0, w_a2m, w_m2a, cth, sth;
double theta, phi0, weight, cth, sth;
ptrdiff_t ofs;
int nph, stride;
} sharp_ringinfo;
@ -127,15 +127,14 @@ 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 wgt_a2m the pixel weight to be used for the ring in alm2map
transforms. Pass NULL to use 1.0 as weight for all rings.
\param wgt_m2a the pixel weight to be used for the ring in map2alm
transforms. Pass NULL to use 1.0 as weight for all rings.
\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
*/
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_a2m, const double *wgt_m2a, sharp_geom_info **geom_info);
const double *wgt, sharp_geom_info **geom_info);
/*! Deallocates the geometry information in \a info. */
void sharp_destroy_geom_info (sharp_geom_info *info);
@ -146,16 +145,23 @@ void sharp_destroy_geom_info (sharp_geom_info *info);
/*! \{ */
/*! Enumeration of SHARP job types. */
typedef enum { SHARP_MAP2ALM, /*!< analysis */
SHARP_ALM2MAP, /*!< synthesis */
SHARP_ALM2MAP_DERIV1 /*!< synthesis of first derivatives */
typedef enum { SHARP_YtW=0, /*!< analysis */
SHARP_MAP2ALM=SHARP_YtW, /*!< analysis */
SHARP_Y=1, /*!< synthesis */
SHARP_ALM2MAP=SHARP_Y, /*!< synthesis */
SHARP_Yt=2, /*!< adjoint synthesis */
SHARP_WY=3, /*!< adjoint analysis */
SHARP_ALM2MAP_DERIV1=4 /*!< synthesis of first derivatives */
} sharp_jobtype;
/*! Job flags */
typedef enum { SHARP_DP = 1<<4, /*!< map and a_lm are in double precision */
SHARP_ADD= 1<<5, /*!< results are added to the output arrays,
instead of overwriting them */
SHARP_NVMAX = (1<<4)-1 /* internal use only */
typedef enum { SHARP_DP = 1<<4,
/*!< map and a_lm are in double precision */
SHARP_ADD = 1<<5,
/*!< results are added to the output arrays, instead of
overwriting them */
SHARP_USE_WEIGHTS = 1<<6, /* internal use only */
SHARP_NVMAX = (1<<4)-1 /* internal use only */
} sharp_jobflags;
/*! Performs a libsharp SHT job. The interface deliberately does not use