Commit 8d10272f authored by Martin Reinecke's avatar Martin Reinecke
Browse files

cleanup

parent 6e14edb3
...@@ -61,8 +61,7 @@ nobase_include_HEADERS = \ ...@@ -61,8 +61,7 @@ nobase_include_HEADERS = \
libsharp2/sharp.h \ libsharp2/sharp.h \
libsharp2/sharp_mpi.h \ libsharp2/sharp_mpi.h \
libsharp2/sharp_geomhelpers.h \ libsharp2/sharp_geomhelpers.h \
libsharp2/sharp_almhelpers.h \ libsharp2/sharp_almhelpers.h
libsharp2/sharp_cxx.h
EXTRA_DIST = \ EXTRA_DIST = \
runtest.sh fortran/sharp.f90 fortran/test_sharp.f90 libsharp2/sharp_mpi.cc runtest.sh fortran/sharp.f90 fortran/test_sharp.f90 libsharp2/sharp_mpi.cc
......
/*
* This file is part of libsharp2.
*
* libsharp2 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.
*
* libsharp2 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 libsharp2; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
/* libsharp2 is being developed at the Max-Planck-Institut fuer Astrophysik */
/*! \file sharp_cxx.h
* Spherical transform library
*
* Copyright (C) 2012-2019 Max-Planck-Society
* \author Martin Reinecke
*/
#ifndef SHARP2_CXX_H
#define SHARP2_CXX_H
#include <complex>
#include "libsharp2/sharp.h"
#include "libsharp2/sharp_geomhelpers.h"
#include "libsharp2/sharp_almhelpers.h"
class sharp_base
{
protected:
sharp_alm_info *ainfo;
sharp_geom_info *ginfo;
public:
sharp_base()
: ainfo(0), ginfo(0) {}
~sharp_base()
{
if (ginfo) sharp_destroy_geom_info(ginfo);
if (ainfo) 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 *wgt)
{
if (ginfo) sharp_destroy_geom_info(ginfo);
sharp_make_geom_info (nrings, nph, ofs, stride, phi0, theta, wgt, &ginfo);
}
void set_ECP_geometry (int nrings, int nphi)
{
if (ginfo) sharp_destroy_geom_info(ginfo);
sharp_make_ecp_geom_info (nrings, nphi, 0., 1, nphi, &ginfo);
}
void set_Gauss_geometry (int nrings, int nphi)
{
if (ginfo) sharp_destroy_geom_info(ginfo);
sharp_make_gauss_geom_info (nrings, nphi, 0., 1, nphi, &ginfo);
}
void set_Healpix_geometry (int nside)
{
if (ginfo) sharp_destroy_geom_info(ginfo);
sharp_make_healpix_geom_info (nside, 1, &ginfo);
}
void set_weighted_Healpix_geometry (int nside, const double *weight)
{
if (ginfo) sharp_destroy_geom_info(ginfo);
sharp_make_weighted_healpix_geom_info (nside, 1, weight, &ginfo);
}
void set_triangular_alm_info (int lmax, int mmax)
{
if (ainfo) sharp_destroy_alm_info(ainfo);
sharp_make_triangular_alm_info (lmax, mmax, 1, &ainfo);
}
};
template<typename T> struct cxxjobhelper__ {};
template<> struct cxxjobhelper__<double>
{ enum {val=SHARP_DP}; };
template<> struct cxxjobhelper__<float>
{ enum {val=0}; };
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
{
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, 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, flags, 0, 0);
}
void alm2map_spin (const T *alm1, const 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,flags, 0, 0);
}
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, flags,0,0);
}
void alm2map_der1 (const 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,flags,0,0);
}
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,flags,0,0);
}
void alm2map_adjoint (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_Yt,0,&aptr,&mptr,ginfo,ainfo,flags,0,0);
}
void alm2map_adjoint (const T *map, std::complex<T> *alm, bool add) const
{
void *aptr=conv(alm), *mptr=conv(map);
int flags=cxxjobhelper__<T>::val | (add ? SHARP_ADD : 0);
sharp_execute (SHARP_Yt,0,&aptr,&mptr,ginfo,ainfo,flags,0,0);
}
void alm2map_spin_adjoint (const T *map1, const T *map2, T *alm1, 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_Yt,spin,aptr,mptr,ginfo,ainfo,flags,0,0);
}
void alm2map_spin_adjoint (const T *map1, const T *map2,
std::complex<T> *alm1, std::complex<T> *alm2, int spin, bool add) const
{
alm2map_spin_adjoint (map1, map2, reinterpret_cast<T *>(alm1),
reinterpret_cast<T *>(alm2), spin, add);
}
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,flags,0,0);
}
void map2alm (const T *map, std::complex<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,flags,0,0);
}
void map2alm_spin (const T *map1, const T *map2, T *alm1, 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,flags,0,0);
}
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,flags,0,0);
}
};
#endif
...@@ -42,12 +42,9 @@ void sharp_make_subset_healpix_geom_info (int nside, int stride, int nrings, ...@@ -42,12 +42,9 @@ void sharp_make_subset_healpix_geom_info (int nside, int stride, int nrings,
ptrdiff_t npix=(ptrdiff_t)nside*nside*12; ptrdiff_t npix=(ptrdiff_t)nside*nside*12;
ptrdiff_t ncap=2*(ptrdiff_t)nside*(nside-1); ptrdiff_t ncap=2*(ptrdiff_t)nside*(nside-1);
double *theta=RALLOC(double,nrings); vector<double> theta(nrings), weight_(nrings), phi0(nrings);
double *weight_=RALLOC(double,nrings); vector<int> nph(nrings), stride_(nrings);
int *nph=RALLOC(int,nrings); vector<ptrdiff_t> ofs(nrings);
double *phi0=RALLOC(double,nrings);
ptrdiff_t *ofs=RALLOC(ptrdiff_t,nrings);
int *stride_=RALLOC(int,nrings);
ptrdiff_t curofs=0, checkofs; /* checkofs used for assertion introduced when adding rings arg */ ptrdiff_t curofs=0, checkofs; /* checkofs used for assertion introduced when adding rings arg */
for (int m=0; m<nrings; ++m) for (int m=0; m<nrings; ++m)
{ {
...@@ -88,15 +85,8 @@ void sharp_make_subset_healpix_geom_info (int nside, int stride, int nrings, ...@@ -88,15 +85,8 @@ void sharp_make_subset_healpix_geom_info (int nside, int stride, int nrings,
curofs+=nph[m]; curofs+=nph[m];
} }
sharp_make_geom_info (nrings, nph, ofs, stride_, phi0, theta, weight_, sharp_make_geom_info (nrings, nph.data(), ofs.data(), stride_.data(), phi0.data(), theta.data(), weight_.data(),
geom_info); geom_info);
DEALLOC(theta);
DEALLOC(weight_);
DEALLOC(nph);
DEALLOC(phi0);
DEALLOC(ofs);
DEALLOC(stride_);
} }
void sharp_make_weighted_healpix_geom_info (int nside, int stride, void sharp_make_weighted_healpix_geom_info (int nside, int stride,
...@@ -110,10 +100,9 @@ void sharp_make_gauss_geom_info (int nrings, int nphi, double phi0, ...@@ -110,10 +100,9 @@ void sharp_make_gauss_geom_info (int nrings, int nphi, double phi0,
{ {
const double pi=3.141592653589793238462643383279502884197; const double pi=3.141592653589793238462643383279502884197;
int *nph=RALLOC(int,nrings); vector<int> nph(nrings), stride_(nrings);
double *phi0_=RALLOC(double,nrings); vector<double> phi0_(nrings);
ptrdiff_t *ofs=RALLOC(ptrdiff_t,nrings); vector<ptrdiff_t> ofs(nrings);
int *stride_=RALLOC(int,nrings);
mr::GL_Integrator integ(nrings); mr::GL_Integrator integ(nrings);
auto theta = integ.coords(); auto theta = integ.coords();
...@@ -128,13 +117,8 @@ void sharp_make_gauss_geom_info (int nrings, int nphi, double phi0, ...@@ -128,13 +117,8 @@ void sharp_make_gauss_geom_info (int nrings, int nphi, double phi0,
weight[m]*=2*pi/nphi; weight[m]*=2*pi/nphi;
} }
sharp_make_geom_info (nrings, nph, ofs, stride_, phi0_, theta.data(), weight.data(), sharp_make_geom_info (nrings, nph.data(), ofs.data(), stride_.data(), phi0_.data(), theta.data(), weight.data(),
geom_info); geom_info);
DEALLOC(nph);
DEALLOC(phi0_);
DEALLOC(ofs);
DEALLOC(stride_);
} }
/* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */ /* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */
...@@ -143,12 +127,9 @@ void sharp_make_fejer1_geom_info (int nrings, int ppring, double phi0, ...@@ -143,12 +127,9 @@ void sharp_make_fejer1_geom_info (int nrings, int ppring, double phi0,
{ {
const double pi=3.141592653589793238462643383279502884197; const double pi=3.141592653589793238462643383279502884197;
double *theta=RALLOC(double,nrings); vector<double> theta(nrings), weight(nrings), phi0_(nrings);
double *weight=RALLOC(double,nrings); vector<int> nph(nrings), stride_(nrings);
int *nph=RALLOC(int,nrings); vector<ptrdiff_t> ofs(nrings);
double *phi0_=RALLOC(double,nrings);
ptrdiff_t *ofs=RALLOC(ptrdiff_t,nrings);
int *stride_=RALLOC(int,nrings);
weight[0]=2.; weight[0]=2.;
for (int k=1; k<=(nrings-1)/2; ++k) for (int k=1; k<=(nrings-1)/2; ++k)
...@@ -157,7 +138,7 @@ void sharp_make_fejer1_geom_info (int nrings, int ppring, double phi0, ...@@ -157,7 +138,7 @@ void sharp_make_fejer1_geom_info (int nrings, int ppring, double phi0,
weight[2*k ]=2./(1.-4.*k*k)*sin((k*pi)/nrings); weight[2*k ]=2./(1.-4.*k*k)*sin((k*pi)/nrings);
} }
if ((nrings&1)==0) weight[nrings-1]=0.; if ((nrings&1)==0) weight[nrings-1]=0.;
mr::r2r_fftpack({size_t(nrings)}, {sizeof(double)}, {sizeof(double)}, {0}, false, false, weight, weight, 1.); mr::r2r_fftpack({size_t(nrings)}, {sizeof(double)}, {sizeof(double)}, {0}, false, false, weight.data(), weight.data(), 1.);
for (int m=0; m<(nrings+1)/2; ++m) for (int m=0; m<(nrings+1)/2; ++m)
{ {
...@@ -171,15 +152,8 @@ void sharp_make_fejer1_geom_info (int nrings, int ppring, double phi0, ...@@ -171,15 +152,8 @@ void sharp_make_fejer1_geom_info (int nrings, int ppring, double phi0,
weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(nrings*nph[m]); weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(nrings*nph[m]);
} }
sharp_make_geom_info (nrings, nph, ofs, stride_, phi0_, theta, weight, sharp_make_geom_info (nrings, nph.data(), ofs.data(), stride_.data(), phi0_.data(), theta.data(), weight.data(),
geom_info); geom_info);
DEALLOC(theta);
DEALLOC(weight);
DEALLOC(nph);
DEALLOC(phi0_);
DEALLOC(ofs);
DEALLOC(stride_);
} }
/* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */ /* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */
...@@ -188,12 +162,9 @@ void sharp_make_cc_geom_info (int nrings, int ppring, double phi0, ...@@ -188,12 +162,9 @@ void sharp_make_cc_geom_info (int nrings, int ppring, double phi0,
{ {
const double pi=3.141592653589793238462643383279502884197; const double pi=3.141592653589793238462643383279502884197;
double *theta=RALLOC(double,nrings); vector<double> theta(nrings), weight(nrings,0.), phi0_(nrings);
vector<double> weight(nrings,0.); vector<int> nph(nrings), stride_(nrings);
int *nph=RALLOC(int,nrings); vector<ptrdiff_t> ofs(nrings);
double *phi0_=RALLOC(double,nrings);
ptrdiff_t *ofs=RALLOC(ptrdiff_t,nrings);
int *stride_=RALLOC(int,nrings);
int n=nrings-1; int n=nrings-1;
double dw=-1./(n*n-1.+(n&1)); double dw=-1./(n*n-1.+(n&1));
...@@ -217,14 +188,8 @@ void sharp_make_cc_geom_info (int nrings, int ppring, double phi0, ...@@ -217,14 +188,8 @@ void sharp_make_cc_geom_info (int nrings, int ppring, double phi0,
weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(n*nph[m]); weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(n*nph[m]);
} }
sharp_make_geom_info (nrings, nph, ofs, stride_, phi0_, theta, weight.data(), sharp_make_geom_info (nrings, nph.data(), ofs.data(), stride_.data(), phi0_.data(), theta.data(), weight.data(),
geom_info); geom_info);
DEALLOC(theta);
DEALLOC(nph);
DEALLOC(phi0_);
DEALLOC(ofs);
DEALLOC(stride_);
} }
/* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */ /* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */
...@@ -233,12 +198,9 @@ void sharp_make_fejer2_geom_info (int nrings, int ppring, double phi0, ...@@ -233,12 +198,9 @@ void sharp_make_fejer2_geom_info (int nrings, int ppring, double phi0,
{ {
const double pi=3.141592653589793238462643383279502884197; const double pi=3.141592653589793238462643383279502884197;
double *theta=RALLOC(double,nrings); vector<double> theta(nrings), weight(nrings+1, 0.), phi0_(nrings);
vector<double> weight(nrings+1, 0.); vector<int> nph(nrings), stride_(nrings);
int *nph=RALLOC(int,nrings); vector<ptrdiff_t> ofs(nrings);
double *phi0_=RALLOC(double,nrings);
ptrdiff_t *ofs=RALLOC(ptrdiff_t,nrings);
int *stride_=RALLOC(int,nrings);
int n=nrings+1; int n=nrings+1;
weight[0]=2.; weight[0]=2.;
...@@ -261,14 +223,8 @@ void sharp_make_fejer2_geom_info (int nrings, int ppring, double phi0, ...@@ -261,14 +223,8 @@ 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]); weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(n*nph[m]);
} }
sharp_make_geom_info (nrings, nph, ofs, stride_, phi0_, theta, weight.data(), sharp_make_geom_info (nrings, nph.data(), ofs.data(), stride_.data(), phi0_.data(), theta.data(), weight.data(),
geom_info); geom_info);
DEALLOC(theta);
DEALLOC(nph);
DEALLOC(phi0_);
DEALLOC(ofs);
DEALLOC(stride_);
} }
void sharp_make_mw_geom_info (int nrings, int ppring, double phi0, void sharp_make_mw_geom_info (int nrings, int ppring, double phi0,
...@@ -276,11 +232,9 @@ void sharp_make_mw_geom_info (int nrings, int ppring, double phi0, ...@@ -276,11 +232,9 @@ void sharp_make_mw_geom_info (int nrings, int ppring, double phi0,
{ {
const double pi=3.141592653589793238462643383279502884197; const double pi=3.141592653589793238462643383279502884197;
double *theta=RALLOC(double,nrings); vector<double> theta(nrings), phi0_(nrings);
int *nph=RALLOC(int,nrings); vector<int> nph(nrings), stride_(nrings);
double *phi0_=RALLOC(double,nrings); vector<ptrdiff_t> ofs(nrings);
ptrdiff_t *ofs=RALLOC(ptrdiff_t,nrings);
int *stride_=RALLOC(int,nrings);
for (int m=0; m<nrings; ++m) for (int m=0; m<nrings; ++m)
{ {
...@@ -292,12 +246,6 @@ void sharp_make_mw_geom_info (int nrings, int ppring, double phi0, ...@@ -292,12 +246,6 @@ void sharp_make_mw_geom_info (int nrings, int ppring, double phi0,
stride_[m]=stride_lon; stride_[m]=stride_lon;
} }
sharp_make_geom_info (nrings, nph, ofs, stride_, phi0_, theta, NULL, sharp_make_geom_info (nrings, nph.data(), ofs.data(), stride_.data(), phi0_.data(), theta.data(), NULL,
geom_info); geom_info);
DEALLOC(theta);
DEALLOC(nph);
DEALLOC(phi0_);
DEALLOC(ofs);
DEALLOC(stride_);
} }
Supports Markdown
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