From 0d8d82677df283038118e9078204306c081c78eb Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Mon, 3 Dec 2012 11:08:03 +0100 Subject: [PATCH] better support for adjoint SHTs --- libsharp/sharp.c | 20 ++++++++++++-------- libsharp/sharp_core_inc2.c | 5 +++++ libsharp/sharp_core_inc3.c | 5 +++++ libsharp/sharp_geomhelpers.c | 13 +++++++------ libsharp/sharp_lowlevel.h | 32 +++++++++++++++++++------------- 5 files changed, 48 insertions(+), 27 deletions(-) diff --git a/libsharp/sharp.c b/libsharp/sharp.c index 8e2cf3a..883eb17 100644 --- a/libsharp/sharp.c +++ b/libsharp/sharp.c @@ -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; mofs]+=creal(self->work[m])*info->w_a2m; + ((double *)data)[m*stride+info->ofs]+=creal(self->work[m])*wgt; else for (int m=0; mofs] += - (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; mwork[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; mwork[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; diff --git a/libsharp/sharp_core_inc2.c b/libsharp/sharp_core_inc2.c index 7672eff..d980232 100644 --- a/libsharp/sharp_core_inc2.c +++ b/libsharp/sharp_core_inc2.c @@ -804,6 +804,11 @@ static void Z(inner_loop) (sharp_job *job, const int *ispair, } break; } + default: + { + UTIL_FAIL("must not happen"); + break; + } } } diff --git a/libsharp/sharp_core_inc3.c b/libsharp/sharp_core_inc3.c index e74bd51..c46780e 100644 --- a/libsharp/sharp_core_inc3.c +++ b/libsharp/sharp_core_inc3.c @@ -794,6 +794,11 @@ static void Y(inner_loop) (sharp_job *job, const int *ispair, } break; } + default: + { + UTIL_FAIL("must not happen"); + break; + } } } diff --git a/libsharp/sharp_geomhelpers.c b/libsharp/sharp_geomhelpers.c index 1653534..a8f2054 100644 --- a/libsharp/sharp_geomhelpers.c +++ b/libsharp/sharp_geomhelpers.c @@ -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); diff --git a/libsharp/sharp_lowlevel.h b/libsharp/sharp_lowlevel.h index 45fe0d3..ac28022 100644 --- a/libsharp/sharp_lowlevel.h +++ b/libsharp/sharp_lowlevel.h @@ -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