try to provide SIMD data types even if no SIMD support was detected in compiler

libsharp2/sharp_core_inc.cc

@@ -52,7 +52,7 @@ using namespace std;
 using Tv=native_simd<double>;
 static constexpr size_t VLEN=Tv::size();
 
-#if (defined(__AVX__) && (!defined(__AVX512F__)))
+#if ((!defined(MRUTIL_NO_SIMD)) && defined(__AVX__) && (!defined(__AVX512F__)))
 static inline void vhsum_cmplx_special (Tv a, Tv b, Tv c, Tv d,
   complex<double> * MRUTIL_RESTRICT cc)
   {

mr_util/simd.h

@@ -16,7 +16,7 @@
  *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  */
 
-/* Copyright (C) 2019 Max-Planck-Society
+/* Copyright (C) 2019-2020 Max-Planck-Society
    Author: Martin Reinecke */
 
 #ifndef MRUTIL_SIMD_H
@@ -41,10 +41,11 @@
 #endif
 #endif
 
-#ifndef MRUTIL_NO_SIMD
 #include <cstdlib>
 #include <cmath>
+#ifndef MRUTIL_NO_SIMD
 #include <x86intrin.h>
+#endif
 
 namespace mr {
 
@@ -227,6 +228,8 @@ template<typename T> class helper_<T,1>
     static size_t maskbits(Tm v) { return v; }
   };
 
+#ifndef MRUTIL_NO_SIMD
+
 #if defined(__AVX512F__)
 template<> class helper_<double,8>
   {
@@ -379,7 +382,11 @@ template<typename T> using native_simd = vtp<T,16/sizeof(T)>;
 #else
 template<typename T> using native_simd = vtp<T,1>;
 #endif
+#else
+template<typename T> using native_simd = vtp<T,1>;
+#endif
 }
+
 using detail_simd::native_simd;
 using detail_simd::reduce;
 using detail_simd::max;
@@ -391,5 +398,3 @@ using detail_simd::all_of;
 }
 
 #endif
-
-#endif

pysharp/pysharp.cc

+/*
+ *  This file is part of pyHealpix.
+ *
+ *  pyHealpix 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.
+ *
+ *  pyHealpix 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 pyHealpix; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ *
+ *  For more information about HEALPix, see http://healpix.sourceforge.net
+ */
+
+/*
+ *  pyHealpix is being developed at the Max-Planck-Institut fuer Astrophysik
+ */
+
+/*
+ *  Copyright (C) 2017-2020 Max-Planck-Society
+ *  Author: Martin Reinecke
+ */
+
+#include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
+#include <vector>
+#include <complex>
+
+#include "libsharp2/sharp.h"
+#include "libsharp2/sharp_geomhelpers.h"
+#include "libsharp2/sharp_almhelpers.h"
+
+using namespace std;
+using namespace mr;
+
+namespace py = pybind11;
+
+namespace {
+
+using a_d_c = py::array_t<double, py::array::c_style | py::array::forcecast>;
+using a_c_c = py::array_t<complex<double>,
+  py::array::c_style | py::array::forcecast>;
+
+template<typename T> class py_sharpjob
+  {
+  private:
+    unique_ptr<sharp_geom_info> ginfo;
+    unique_ptr<sharp_alm_info> ainfo;
+    int64_t lmax_, mmax_, npix_;
+
+  public:
+    py_sharpjob () : lmax_(0), mmax_(0), npix_(0) {}
+
+    string repr() const
+      {
+      return "<sharpjob_d: lmax=" + dataToString(lmax_) +
+        ", mmax=" + dataToString(mmax_) + ", npix=", dataToString(npix_) +".>";
+      }
+
+    void set_Gauss_geometry(int64_t nrings, int64_t nphi)
+      {
+      MR_assert((nrings>0)&&(nphi>0),"bad grid dimensions");
+      npix_=nrings*nphi;
+      ginfo = sharp_make_gauss_geom_info (nrings, nphi, 0., 1, nphi);
+      }
+    void set_Healpix_geometry(int64_t nside)
+      {
+      MR_assert(nside>0,"bad Nside value");
+      npix_=12*nside*nside;
+      ginfo = sharp_make_healpix_geom_info (nside, 1);
+      }
+    void set_ECP_geometry(int64_t nrings, int64_t nphi)
+      {
+      MR_assert(nrings>0,"bad nrings value");
+      MR_assert(nphi>0,"bad nphi value");
+      npix_=nrings*nphi;
+      ginfo = sharp_make_ecp_geom_info (nrings, nphi, 0., 1, nphi);
+      }
+    void set_triangular_alm_info (int64_t lmax, int64_t mmax)
+      {
+      MR_assert(mmax>=0,"negative mmax");
+      MR_assert(mmax<=lmax,"mmax must not be larger than lmax");
+      lmax_=lmax; mmax_=mmax;
+      ainfo = sharp_make_triangular_alm_info(lmax,mmax,1);
+      }
+
+    int64_t n_alm() const
+      { return ((mmax_+1)*(mmax_+2))/2 + (mmax_+1)*(lmax_-mmax_); }
+
+    a_d_c alm2map (const a_c_c &alm) const
+      {
+      MR_assert(npix_>0,"no map geometry specified");
+      MR_assert (alm.size()==n_alm(),
+        "incorrect size of a_lm array");
+      a_d_c map(npix_);
+      auto mr=map.mutable_unchecked<1>();
+      auto ar=alm.unchecked<1>();
+      sharp_alm2map(&ar[0], &mr[0], *ginfo, *ainfo, 0, nullptr, nullptr);
+      return map;
+      }
+    a_c_c alm2map_adjoint (const a_d_c &map) const
+      {
+      MR_assert(npix_>0,"no map geometry specified");
+      MR_assert (map.size()==npix_,"incorrect size of map array");
+      a_c_c alm(n_alm());
+      auto mr=map.unchecked<1>();
+      auto ar=alm.mutable_unchecked<1>();
+      sharp_map2alm(&ar[0], &mr[0], *ginfo, *ainfo, 0, nullptr, nullptr);
+      return alm;
+      }
+    a_c_c map2alm (const a_d_c &map) const
+      {
+      MR_assert(npix_>0,"no map geometry specified");
+      MR_assert (map.size()==npix_,"incorrect size of map array");
+      a_c_c alm(n_alm());
+      auto mr=map.unchecked<1>();
+      auto ar=alm.mutable_unchecked<1>();
+      sharp_map2alm(&ar[0], &mr[0], *ginfo, *ainfo, SHARP_USE_WEIGHTS, nullptr, nullptr);
+      return alm;
+      }
+    a_d_c alm2map_spin (const a_c_c &alm, int64_t spin) const
+      {
+      MR_assert(npix_>0,"no map geometry specified");
+      auto ar=alm.unchecked<2>();
+      MR_assert((ar.shape(0)==2)&&(ar.shape(1)==n_alm()),
+        "incorrect size of a_lm array");
+      a_d_c map(vector<size_t>{2,size_t(npix_)});
+      auto mr=map.mutable_unchecked<2>();
+      sharp_alm2map_spin(spin, &ar(0,0), &ar(1,0), &mr(0,0), &mr(1,0), *ginfo, *ainfo, 0, nullptr, nullptr);
+      return map;
+      }
+    a_c_c map2alm_spin (const a_d_c &map, int64_t spin) const
+      {
+      MR_assert(npix_>0,"no map geometry specified");
+      auto mr=map.unchecked<2>();
+      MR_assert ((mr.shape(0)==2)&&(mr.shape(1)==npix_),
+        "incorrect size of map array");
+      a_c_c alm(vector<size_t>{2,size_t(n_alm())});
+      auto ar=alm.mutable_unchecked<2>();
+      sharp_map2alm_spin(spin, &ar(0,0), &ar(1,0), &mr(0,0), &mr(1,0), *ginfo, *ainfo, SHARP_USE_WEIGHTS, nullptr, nullptr);
+      return alm;
+      }
+  };
+
+a_d_c GL_weights(int64_t nlat, int64_t nlon)
+  {
+  a_d_c res(nlat);
+  auto rr=res.mutable_unchecked<1>();
+  GL_Integrator integ(nlat);
+  auto wgt = integ.weights();
+  for (size_t i=0; i<size_t(rr.shape(0)); ++i)
+    rr[i]=wgt[i]*twopi/nlon;
+  return res;
+  }
+
+a_d_c GL_thetas(int64_t nlat)
+  {
+  a_d_c res(nlat);
+  auto rr=res.mutable_unchecked<1>();
+  GL_Integrator integ(nlat);
+  auto coord = integ.coords();
+  for (size_t i=0; i<size_t(rr.shape(0)); ++i)
+    rr[i]=acos(-coord[i]);
+  return res;
+  }
+
+const char *pyHealpix_DS = R"DELIM(
+Python interface for some of the libsharp functionality
+
+Error conditions are reported by raising exceptions.
+)DELIM";
+
+} // unnamed namespace
+
+PYBIND11_MODULE(pysharp, m)
+  {
+  using namespace pybind11::literals;
+
+  m.doc() = pysharp_DS;
+
+  py::class_<py_sharpjob<double>> (m, "sharpjob_d")
+    .def(py::init<>())
+    .def("set_Gauss_geometry", &py_sharpjob<double>::set_Gauss_geometry,
+      "nrings"_a,"nphi"_a)
+    .def("set_Healpix_geometry", &py_sharpjob<double>::set_Healpix_geometry,
+      "nside"_a)
+    .def("set_ECP_geometry", &py_sharpjob<double>::set_ECP_geometry,
+      "nrings"_a, "nphi"_a)
+    .def("set_triangular_alm_info",
+      &py_sharpjob<double>::set_triangular_alm_info, "lmax"_a, "mmax"_a)
+    .def("n_alm", &py_sharpjob<double>::n_alm)
+    .def("alm2map", &py_sharpjob<double>::alm2map,"alm"_a)
+    .def("alm2map_adjoint", &py_sharpjob<double>::alm2map_adjoint,"map"_a)
+    .def("map2alm", &py_sharpjob<double>::map2alm,"map"_a)
+    .def("alm2map_spin", &py_sharpjob<double>::alm2map_spin,"alm"_a,"spin"_a)
+    .def("map2alm_spin", &py_sharpjob<double>::map2alm_spin,"map"_a,"spin"_a)
+    .def("__repr__", &py_sharpjob<double>::repr)
+    ;
+
+  m.def("GL_weights",&GL_weights, "nlat"_a, "nlon"_a);
+  m.def("GL_thetas",&GL_thetas, "nlat"_a);
+  }

pysharp/setup.py

+from setuptools import setup, Extension
+import sys
+
+
+class _deferred_pybind11_include(object):
+    def __init__(self, user=False):
+        self.user = user
+
+    def __str__(self):
+        import pybind11
+        return pybind11.get_include(self.user)
+
+
+include_dirs = ['..', _deferred_pybind11_include(True),
+                _deferred_pybind11_include()]
+extra_compile_args = ['--std=c++17', '-march=native', '-ffast-math', '-O3']
+python_module_link_args = []
+define_macros = []
+
+if sys.platform == 'darwin':
+    import distutils.sysconfig
+    extra_compile_args += ['-mmacosx-version-min=10.9']
+    python_module_link_args += ['-mmacosx-version-min=10.9', '-bundle']
+    vars = distutils.sysconfig.get_config_vars()
+    vars['LDSHARED'] = vars['LDSHARED'].replace('-bundle', '')
+elif sys.platform == 'win32':
+    extra_compile_args = ['/Ox', '/EHsc']
+else:
+    extra_compile_args += ['-Wfatal-errors', '-Wfloat-conversion', '-Wsign-conversion', '-Wconversion' ,'-W', '-Wall', '-Wstrict-aliasing=2', '-Wwrite-strings', '-Wredundant-decls', '-Woverloaded-virtual', '-Wcast-qual', '-Wcast-align', '-Wpointer-arith']
+    python_module_link_args += ['-march=native', '-ffast-math', '-Wl,-rpath,$ORIGIN']
+
+# if you don't want debugging info, add "-s" to python_module_link_args
+
+def get_extension_modules():
+    return [Extension('pysharp',
+                      language='c++',
+                      sources=['pysharp.cc','../mr_util/threading.cc',
+                               '../libsharp2/sharp.cc', '../libsharp2/sharp_core.cc', '../libsharp2/sharp_geomhelpers.cc',
+                               '../libsharp2/sharp_almhelpers.cc','../libsharp2/sharp_ylmgen.cc'],
+                      depends=['../mr_util/fft.h', '../mr_util/mav.h', '../mr_util/threading.h',
+                               '../mr_util/aligned_array.h', '../mr_util/simd.h',
+                               '../mr_util/cmplx.h', '../mr_util/unity_roots.h', '../mr_util/error_handling.h',
+                               'setup.py'],
+
+                      include_dirs=include_dirs,
+                      define_macros=define_macros,
+                      extra_compile_args=extra_compile_args,
+                      extra_link_args=python_module_link_args)]
+
+
+setup(name='pysharp',
+      version='0.0.1',
+      description='Python bindings for some libsharp functionality',
+      include_package_data=True,
+      author='Martin Reinecke',
+      author_email='martin@mpa-garching.mpg.de',
+      packages=[],
+      setup_requires=['numpy>=1.15.0', 'pybind11>=2.2.4'],
+      ext_modules=get_extension_modules(),
+      install_requires=['numpy>=1.15.0', 'pybind11>=2.2.4']
+      )