diff --git a/libsharp/sharp.c b/libsharp/sharp.c index 917a09e..6daada9 100644 --- a/libsharp/sharp.c +++ b/libsharp/sharp.c @@ -584,6 +584,18 @@ void sharps_build_job (sharp_job *job, sharp_jobtype type, int spin, job->fde=FLOAT; } +void sharp_build_job_ll (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) + { + if (dp) + sharpd_build_job (job, type, spin, add_output, (dcmplx **)alm, + (double **)map, geom_info, alm_info, ntrans); + else + sharps_build_job (job, type, spin, add_output, (fcmplx **)alm, + (float **)map, geom_info, alm_info, ntrans); + } + int sharp_get_nv_max (void) { return 6; } @@ -667,4 +679,6 @@ int sharp_nv_oracle (sharp_jobtype type, int spin, int ntrans) return nv_opt[ntrans-1][spin!=0][type]; } +#ifdef USE_MPI #include "sharp_mpi.c" +#endif diff --git a/libsharp/sharp.h b/libsharp/sharp.h index a984336..ab7ecf9 100644 --- a/libsharp/sharp.h +++ b/libsharp/sharp.h @@ -195,6 +195,10 @@ void sharps_build_job (sharp_job *job, sharp_jobtype type, int spin, int add_output, complex float **alm, float **map, const sharp_geom_info *geom_info, const sharp_alm_info *alm_info, int ntrans); +void sharp_build_job_ll (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); + /*! Execute the SHT job previously constructed by sharpd_build_job() or sharps_build_job(). */ void sharp_execute_job (sharp_job *job); diff --git a/libsharp/sharp_cxx.h b/libsharp/sharp_cxx.h new file mode 100644 index 0000000..5a513c0 --- /dev/null +++ b/libsharp/sharp_cxx.h @@ -0,0 +1,143 @@ +/* + * This file is part of libsharp. + * + * libsharp is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * libsharp is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with libsharp; if not, write to the Free Software + * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + */ + +/* + * libsharp is being developed at the Max-Planck-Institut fuer Astrophysik + * and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt + * (DLR). + */ + +/*! \file sharp_cxx.h + * Spherical transform library + * + * Copyright (C) 2008, 2009 Max-Planck-Society + * \author Martin Reinecke + */ + +#ifndef PLANCK_SHARP_CXX_H +#define PLANCK_SHARP_CXX_H + +#include "sharp.h" +#include "sharp_geomhelpers.h" +#include "sharp_almhelpers.h" +#include "xcomplex.h" + +class sharp_base + { + protected: + sharp_job job; + sharp_alm_info *ainfo; + sharp_geom_info *ginfo; + + public: + sharp_base() + : ainfo(0), ginfo(0) {} + ~sharp_base() + { + sharp_destroy_geom_info(ginfo); + sharp_destroy_alm_info(ainfo); + } + + void set_general_geometry (int nrings, const int *nph, const ptrdiff_t *ofs, + const int *stride, const double *phi0, const double *theta, + const double *weight) + { + sharp_make_geom_info (nrings, nph, ofs, stride, phi0, theta, weight, + &ginfo); + } + + void set_ECP_geometry (int nrings, int nphi) + { sharp_make_ecp_geom_info (nrings, nphi, 0., 1, nphi, &ginfo); } + + void set_Healpix_geometry (int nside) + { sharp_make_healpix_geom_info (nside, 1, &ginfo); } + + void set_weighted_Healpix_geometry (int nside, const double *weight) + { sharp_make_weighted_healpix_geom_info (nside, 1, weight, &ginfo); } + + void set_triangular_alm_info (int lmax, int mmax) + { sharp_make_triangular_alm_info (lmax, mmax, 1, &ainfo); } + }; + +template struct cxxjobhelper__ {}; + +template<> struct cxxjobhelper__ + { enum {val=1}; }; + +template<> struct cxxjobhelper__ + { enum {val=0}; }; + + +template class sharp_cxxjob: public sharp_base + { + private: + static void *conv (xcomplex *ptr) + { return reinterpret_cast(ptr); } + static void *conv (const xcomplex *ptr) + { return const_cast(reinterpret_cast(ptr)); } + static void *conv (T *ptr) + { return reinterpret_cast(ptr); } + static void *conv (const T *ptr) + { return const_cast(reinterpret_cast(ptr)); } + + public: + void alm2map (const xcomplex *alm, T *map, bool add) + { + void *aptr=conv(alm), *mptr=conv(map); + sharp_build_job_ll (&job, SHARP_ALM2MAP, 0, add, &aptr, &mptr, ginfo, + ainfo, 1, cxxjobhelper__::val); + sharp_execute_job(&job); + } + void alm2map_spin (const xcomplex *alm1, const xcomplex *alm2, + T *map1, T *map2, int spin, bool add) + { + void *aptr[2], *mptr[2]; + aptr[0]=conv(alm1); aptr[1]=conv(alm2); + mptr[0]=conv(map1); mptr[1]=conv(map2); + sharp_build_job_ll (&job, SHARP_ALM2MAP, spin, add, aptr, mptr, ginfo, + ainfo, 1, cxxjobhelper__::val); + sharp_execute_job(&job); + } + void alm2map_der1 (const xcomplex *alm, T *map1, T *map2, bool add) + { + void *aptr=conv(alm), *mptr[2]; + mptr[0]=conv(map1); mptr[1]=conv(map2); + sharp_build_job_ll (&job, SHARP_ALM2MAP_DERIV1, 1, add,&aptr, mptr, ginfo, + ainfo, 1, cxxjobhelper__::val); + sharp_execute_job(&job); + } + void map2alm (const T *map, xcomplex *alm, bool add) + { + void *aptr=conv(alm), *mptr=conv(map); + sharp_build_job_ll (&job, SHARP_MAP2ALM, 0, add, &aptr, &mptr, ginfo, + ainfo, 1, cxxjobhelper__::val); + sharp_execute_job(&job); + } + void map2alm_spin (const T *map1, const T *map2, xcomplex *alm1, + xcomplex *alm2, int spin, bool add) + { + void *aptr[2], *mptr[2]; + aptr[0]=conv(alm1); aptr[1]=conv(alm2); + mptr[0]=conv(map1); mptr[1]=conv(map2); + sharp_build_job_ll (&job, SHARP_MAP2ALM, spin, add, aptr, mptr, ginfo, + ainfo, 1, cxxjobhelper__::val); + sharp_execute_job(&job); + } + }; + +#endif