Commit 9c97d150 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

sync: more support for adjoint SHTs

parent 3c24364f
......@@ -25,7 +25,7 @@
*/
/*
* Copyright (C) 2003-2016 Max-Planck-Society
* Copyright (C) 2003-2017 Max-Planck-Society
* Author: Martin Reinecke
*/
......@@ -71,7 +71,7 @@ template void map2alm (const Healpix_Map<double> &map,
bool add_alm);
template<typename T> void alm2map_adjoint (const Healpix_Map<T> &map,
Alm<xcomplex<T> > &alm)
Alm<xcomplex<T> > &alm, bool add_alm)
{
planck_assert (map.Scheme()==RING,
"alm2map_adjoint: map must be in RING scheme");
......@@ -81,13 +81,13 @@ template<typename T> void alm2map_adjoint (const Healpix_Map<T> &map,
sharp_cxxjob<T> job;
job.set_Healpix_geometry (map.Nside());
job.set_triangular_alm_info (alm.Lmax(), alm.Mmax());
job.alm2map_adjoint(&map[0], &alm(0,0), false);
job.alm2map_adjoint(&map[0], &alm(0,0), add_alm);
}
template void alm2map_adjoint (const Healpix_Map<float> &map,
Alm<xcomplex<float> > &alm);
Alm<xcomplex<float> > &alm, bool add_alm);
template void alm2map_adjoint (const Healpix_Map<double> &map,
Alm<xcomplex<double> > &alm);
Alm<xcomplex<double> > &alm, bool add_alm);
template<typename T> void map2alm_iter (const Healpix_Map<T> &map,
Alm<xcomplex<T> > &alm, int num_iter, const arr<double> &weight)
......@@ -171,6 +171,37 @@ template void map2alm_spin
Alm<xcomplex<double> > &alm1, Alm<xcomplex<double> > &alm2,
int spin, const arr<double> &weight, bool add_alm);
template<typename T> void alm2map_spin_adjoint
(const Healpix_Map<T> &map1, const Healpix_Map<T> &map2,
Alm<xcomplex<T> > &alm1, Alm<xcomplex<T> > &alm2,
int spin, bool add_alm)
{
planck_assert (map1.Scheme()==RING,
"alm2map_spin_adjoint: maps must be in RING scheme");
planck_assert (map1.conformable(map2),
"alm2map_spin_adjoint: maps are not conformable");
planck_assert (alm1.conformable(alm1),
"alm2map_spin_adjoint: a_lm are not conformable");
planck_assert (map1.fullyDefined()&&map2.fullyDefined(),
"map contains undefined pixels");
checkLmaxNside(alm1.Lmax(), map1.Nside());
sharp_cxxjob<T> job;
job.set_Healpix_geometry (map1.Nside());
job.set_triangular_alm_info (alm1.Lmax(), alm1.Mmax());
job.alm2map_spin_adjoint(&map1[0],&map2[0],&alm1(0,0),&alm2(0,0),spin,
add_alm);
}
template void alm2map_spin_adjoint
(const Healpix_Map<float> &map1, const Healpix_Map<float> &map2,
Alm<xcomplex<float> > &alm1, Alm<xcomplex<float> > &alm2,
int spin, bool add_alm);
template void alm2map_spin_adjoint
(const Healpix_Map<double> &map1, const Healpix_Map<double> &map2,
Alm<xcomplex<double> > &alm1, Alm<xcomplex<double> > &alm2,
int spin, bool add_alm);
template<typename T> void map2alm_spin_iter2
(const Healpix_Map<T> &map1, const Healpix_Map<T> &map2,
Alm<xcomplex<T> > &alm1, Alm<xcomplex<T> > &alm2,
......@@ -254,6 +285,48 @@ template void map2alm_pol
const arr<double> &weight,
bool add_alm);
template<typename T> void alm2map_pol_adjoint
(const Healpix_Map<T> &mapT,
const Healpix_Map<T> &mapQ,
const Healpix_Map<T> &mapU,
Alm<xcomplex<T> > &almT,
Alm<xcomplex<T> > &almG,
Alm<xcomplex<T> > &almC,
bool add_alm)
{
planck_assert (mapT.Scheme()==RING,
"alm2map_pol_adjoint: maps must be in RING scheme");
planck_assert (mapT.conformable(mapQ) && mapT.conformable(mapU),
"alm2map_pol_adjoint: maps are not conformable");
planck_assert (almT.conformable(almG) && almT.conformable(almC),
"alm2map_pol_adjoint: a_lm are not conformable");
planck_assert (mapT.fullyDefined()&&mapQ.fullyDefined()&&mapU.fullyDefined(),
"map contains undefined pixels");
checkLmaxNside(almT.Lmax(), mapT.Nside());
sharp_cxxjob<T> job;
job.set_Healpix_geometry (mapT.Nside());
job.set_triangular_alm_info (almT.Lmax(), almT.Mmax());
job.alm2map_adjoint(&mapT[0], &almT(0,0), add_alm);
job.alm2map_spin_adjoint(&mapQ[0],&mapU[0],&almG(0,0),&almC(0,0),2,add_alm);
}
template void alm2map_pol_adjoint
(const Healpix_Map<float> &mapT,
const Healpix_Map<float> &mapQ,
const Healpix_Map<float> &mapU,
Alm<xcomplex<float> > &almT,
Alm<xcomplex<float> > &almG,
Alm<xcomplex<float> > &almC,
bool add_alm);
template void alm2map_pol_adjoint
(const Healpix_Map<double> &mapT,
const Healpix_Map<double> &mapQ,
const Healpix_Map<double> &mapU,
Alm<xcomplex<double> > &almT,
Alm<xcomplex<double> > &almG,
Alm<xcomplex<double> > &almC,
bool add_alm);
template<typename T> void map2alm_pol_iter
(const Healpix_Map<T> &mapT,
......@@ -351,24 +424,24 @@ template void map2alm_pol_iter2
template<typename T> void alm2map (const Alm<xcomplex<T> > &alm,
Healpix_Map<T> &map)
Healpix_Map<T> &map, bool add_map)
{
planck_assert (map.Scheme()==RING, "alm2map: map must be in RING scheme");
sharp_cxxjob<T> job;
job.set_Healpix_geometry (map.Nside());
job.set_triangular_alm_info (alm.Lmax(), alm.Mmax());
job.alm2map(&alm(0,0), &map[0], false);
job.alm2map(&alm(0,0), &map[0], add_map);
}
template void alm2map (const Alm<xcomplex<double> > &alm,
Healpix_Map<double> &map);
Healpix_Map<double> &map, bool add_map);
template void alm2map (const Alm<xcomplex<float> > &alm,
Healpix_Map<float> &map);
Healpix_Map<float> &map, bool add_map);
template<typename T> void alm2map_spin
(const Alm<xcomplex<T> > &alm1, const Alm<xcomplex<T> > &alm2,
Healpix_Map<T> &map1, Healpix_Map<T> &map2, int spin)
Healpix_Map<T> &map1, Healpix_Map<T> &map2, int spin, bool add_map)
{
planck_assert (map1.Scheme()==RING,
"alm2map_spin: maps must be in RING scheme");
......@@ -380,15 +453,15 @@ template<typename T> void alm2map_spin
sharp_cxxjob<T> job;
job.set_Healpix_geometry (map1.Nside());
job.set_triangular_alm_info (alm1.Lmax(), alm1.Mmax());
job.alm2map_spin(&alm1(0,0),&alm2(0,0),&map1[0],&map2[0],spin,false);
job.alm2map_spin(&alm1(0,0),&alm2(0,0),&map1[0],&map2[0],spin,add_map);
}
template void alm2map_spin
(const Alm<xcomplex<double> > &alm1, const Alm<xcomplex<double> > &alm2,
Healpix_Map<double> &map, Healpix_Map<double> &map2, int spin);
Healpix_Map<double> &map, Healpix_Map<double> &map2, int spin, bool add_map);
template void alm2map_spin
(const Alm<xcomplex<float> > &alm1, const Alm<xcomplex<float> > &alm2,
Healpix_Map<float> &map, Healpix_Map<float> &map2, int spin);
Healpix_Map<float> &map, Healpix_Map<float> &map2, int spin, bool add_map);
template<typename T> void alm2map_pol
......@@ -397,7 +470,8 @@ template<typename T> void alm2map_pol
const Alm<xcomplex<T> > &almC,
Healpix_Map<T> &mapT,
Healpix_Map<T> &mapQ,
Healpix_Map<T> &mapU)
Healpix_Map<T> &mapU,
bool add_map)
{
planck_assert (mapT.Scheme()==RING,
"alm2map_pol: maps must be in RING scheme");
......@@ -409,8 +483,8 @@ template<typename T> void alm2map_pol
sharp_cxxjob<T> job;
job.set_Healpix_geometry (mapT.Nside());
job.set_triangular_alm_info (almT.Lmax(), almT.Mmax());
job.alm2map(&almT(0,0), &mapT[0], false);
job.alm2map_spin(&almG(0,0), &almC(0,0), &mapQ[0], &mapU[0], 2, false);
job.alm2map(&almT(0,0), &mapT[0], add_map);
job.alm2map_spin(&almG(0,0), &almC(0,0), &mapQ[0], &mapU[0], 2, add_map);
}
template void alm2map_pol (const Alm<xcomplex<double> > &almT,
......@@ -418,14 +492,16 @@ template void alm2map_pol (const Alm<xcomplex<double> > &almT,
const Alm<xcomplex<double> > &almC,
Healpix_Map<double> &mapT,
Healpix_Map<double> &mapQ,
Healpix_Map<double> &mapU);
Healpix_Map<double> &mapU,
bool add_map);
template void alm2map_pol (const Alm<xcomplex<float> > &almT,
const Alm<xcomplex<float> > &almG,
const Alm<xcomplex<float> > &almC,
Healpix_Map<float> &mapT,
Healpix_Map<float> &mapQ,
Healpix_Map<float> &mapU);
Healpix_Map<float> &mapU,
bool add_map);
template<typename T> void alm2map_der1
......
......@@ -108,12 +108,8 @@ template<typename T> class sharp_cxxjob: public sharp_base
private:
static void *conv (T *ptr)
{ return reinterpret_cast<void *>(ptr); }
static void *conv (std::complex<T> *ptr)
{ return reinterpret_cast<void *>(ptr); }
static void *conv (const T *ptr)
{ return const_cast<void *>(reinterpret_cast<const void *>(ptr)); }
static void *conv (const std::complex<T> *ptr)
{ return const_cast<void *>(reinterpret_cast<const void *>(ptr)); }
public:
void alm2map (const T *alm, T *map, bool add) const
......@@ -124,12 +120,7 @@ template<typename T> class sharp_cxxjob: public sharp_base
flags,0,0);
}
void alm2map (const std::complex<T> *alm, T *map, bool add) const
{
void *aptr=conv(alm), *mptr=conv(map);
int flags=cxxjobhelper__<T>::val | (add ? SHARP_ADD : 0);
sharp_execute (SHARP_ALM2MAP, 0, &aptr, &mptr, ginfo, ainfo, 1,
flags,0,0);
}
{ alm2map (reinterpret_cast<const T *>(alm),map,add); }
void alm2map_spin (const T *alm1, const T *alm2,
T *map1, T *map2, int spin, bool add) const
{
......@@ -142,11 +133,8 @@ template<typename T> class sharp_cxxjob: public sharp_base
void alm2map_spin (const std::complex<T> *alm1, const std::complex<T> *alm2,
T *map1, T *map2, int spin, bool add) const
{
void *aptr[2], *mptr[2];
aptr[0]=conv(alm1); aptr[1]=conv(alm2);
mptr[0]=conv(map1); mptr[1]=conv(map2);
int flags=cxxjobhelper__<T>::val | (add ? SHARP_ADD : 0);
sharp_execute (SHARP_ALM2MAP,spin,aptr,mptr,ginfo,ainfo,1,flags,0,0);
alm2map_spin (reinterpret_cast<const T *>(alm1),
reinterpret_cast<const T *>(alm2), map1, map2, spin, add);
}
void alm2map_der1 (const T *alm, T *map1, T *map2, bool add) const
{
......@@ -157,12 +145,7 @@ template<typename T> class sharp_cxxjob: public sharp_base
}
void alm2map_der1 (const std::complex<T> *alm, T *map1, T *map2, bool add)
const
{
void *aptr=conv(alm), *mptr[2];
mptr[0]=conv(map1); mptr[1]=conv(map2);
int flags=cxxjobhelper__<T>::val | (add ? SHARP_ADD : 0);
sharp_execute (SHARP_ALM2MAP_DERIV1,1,&aptr,mptr,ginfo,ainfo,1,flags,0,0);
}
{ alm2map_der1(reinterpret_cast<const T *>(alm), map1, map2, add); }
void alm2map_adjoint (const T *map, T *alm, bool add) const
{
void *aptr=conv(alm), *mptr=conv(map);
......@@ -170,23 +153,30 @@ template<typename T> class sharp_cxxjob: public sharp_base
sharp_execute (SHARP_Yt,0,&aptr,&mptr,ginfo,ainfo,1,flags,0,0);
}
void alm2map_adjoint (const T *map, std::complex<T> *alm, bool add) const
{ alm2map_adjoint (map, reinterpret_cast<T *>(alm), add); }
void alm2map_spin_adjoint (const T *map1, const T *map2, T *alm1, T *alm2,
int spin, bool add) const
{
void *aptr=conv(alm), *mptr=conv(map);
void *aptr[2], *mptr[2];
aptr[0]=conv(alm1); aptr[1]=conv(alm2);
mptr[0]=conv(map1); mptr[1]=conv(map2);
int flags=cxxjobhelper__<T>::val | (add ? SHARP_ADD : 0);
sharp_execute (SHARP_Yt,0,&aptr,&mptr,ginfo,ainfo,1,flags,0,0);
sharp_execute (SHARP_Yt,spin,aptr,mptr,ginfo,ainfo,1,flags,0,0);
}
void map2alm (const T *map, T *alm, bool add) const
void alm2map_spin_adjoint (const T *map1, const T *map2,
std::complex<T> *alm1, std::complex<T> *alm2, int spin, bool add) const
{
void *aptr=conv(alm), *mptr=conv(map);
int flags=cxxjobhelper__<T>::val | (add ? SHARP_ADD : 0);
sharp_execute (SHARP_MAP2ALM,0,&aptr,&mptr,ginfo,ainfo,1,flags,0,0);
alm2map_spin_adjoint (map1, map2, reinterpret_cast<T *>(alm1),
reinterpret_cast<T *>(alm2), spin, add);
}
void map2alm (const T *map, std::complex<T> *alm, bool add) const
void map2alm (const T *map, T *alm, bool add) const
{
void *aptr=conv(alm), *mptr=conv(map);
int flags=cxxjobhelper__<T>::val | (add ? SHARP_ADD : 0);
sharp_execute (SHARP_MAP2ALM,0,&aptr,&mptr,ginfo,ainfo,1,flags,0,0);
}
void map2alm (const T *map, std::complex<T> *alm, bool add) const
{ map2alm (map, reinterpret_cast<T *>(alm),add); }
void map2alm_spin (const T *map1, const T *map2, T *alm1, T *alm2,
int spin, bool add) const
{
......@@ -199,11 +189,8 @@ template<typename T> class sharp_cxxjob: public sharp_base
void map2alm_spin (const T *map1, const T *map2, std::complex<T> *alm1,
std::complex<T> *alm2, int spin, bool add) const
{
void *aptr[2], *mptr[2];
aptr[0]=conv(alm1); aptr[1]=conv(alm2);
mptr[0]=conv(map1); mptr[1]=conv(map2);
int flags=cxxjobhelper__<T>::val | (add ? SHARP_ADD : 0);
sharp_execute (SHARP_MAP2ALM,spin,aptr,mptr,ginfo,ainfo,1,flags,0,0);
map2alm_spin (map1, map2, reinterpret_cast<T *>(alm1),
reinterpret_cast<T *>(alm2), spin, add);
}
};
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment