evolution

src/error_handling.cc

@@ -25,9 +25,7 @@ using namespace std;
 
 namespace mr {
 
-namespace error_handling {
-
-namespace detail {
+namespace detail_error_handling {
 
 bool abort_in_progress__ = false;
 
@@ -45,4 +43,4 @@ void killjob__()
   exit(1);
   }
 
-}}}
+}}

src/mr_util/aligned_array.h

+/*
+ *  This file is part of the MR utility library.
+ *
+ *  This code 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.
+ *
+ *  This code 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 this code; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+/* Copyright (C) 2019 Max-Planck-Society
+   Author: Martin Reinecke */
+
+#ifndef MRUTIL_ALIGNED_ARRAY_H
+#define MRUTIL_ALIGNED_ARRAY_H
+
+#include <cstdlib>
+#include <memory>
+
+namespace mr {
+
+namespace detail_aligned_array {
+
+using namespace std;
+
+/*! Simple array class guaranteeing 64-byte alignment of the data pointer.
+    Mostly useful for storing data accessed by SIMD instructions. */
+template<typename T> class aligned_array
+  {
+  private:
+    T *p;
+    size_t sz;
+
+#if (__cplusplus >= 201703L)
+    static T *ralloc(size_t num)
+      {
+      if (num==0) return nullptr;
+      void *res = aligned_alloc(64,num*sizeof(T));
+      if (!res) throw bad_alloc();
+      return reinterpret_cast<T *>(res);
+      }
+    static void dealloc(T *ptr)
+      { free(ptr); }
+#else // portable emulation
+    static T *ralloc(size_t num)
+      {
+      if (num==0) return nullptr;
+      void *ptr = malloc(num*sizeof(T)+64);
+      if (!ptr) throw bad_alloc();
+      T *res = reinterpret_cast<T *>
+        ((reinterpret_cast<size_t>(ptr) & ~(size_t(63))) + 64);
+      (reinterpret_cast<void**>(res))[-1] = ptr;
+      return res;
+      }
+    static void dealloc(T *ptr)
+      { if (ptr) free((reinterpret_cast<void**>(ptr))[-1]); }
+#endif
+
+  public:
+    aligned_array() : p(0), sz(0) {}
+    aligned_array(size_t n) : p(ralloc(n)), sz(n) {}
+    aligned_array(aligned_array &&other)
+      : p(other.p), sz(other.sz)
+      { other.p=nullptr; other.sz=0; }
+    ~aligned_array() { dealloc(p); }
+
+    void resize(size_t n)
+      {
+      if (n==sz) return;
+      dealloc(p);
+      p = ralloc(n);
+      sz = n;
+      }
+
+    T &operator[](size_t idx) { return p[idx]; }
+    const T &operator[](size_t idx) const { return p[idx]; }
+
+    T *data() { return p; }
+    const T *data() const { return p; }
+
+    size_t size() const { return sz; }
+  };
+
+}
+
+using detail_aligned_array::aligned_array;
+
+}
+
+#endif
+

src/mr_util/error_handling.h

@@ -27,23 +27,21 @@
 
 namespace mr {
 
-namespace error_handling {
-
-namespace detail {
+namespace detail_error_handling {
 
 #if defined (__GNUC__)
-#define MR_ERROR_HANDLING_LOC_ ::mr::error_handling::detail::CodeLocation(__FILE__, __LINE__, __PRETTY_FUNCTION__)
+#define MRUTIL_ERROR_HANDLING_LOC_ ::mr::detail_error_handling::CodeLocation(__FILE__, __LINE__, __PRETTY_FUNCTION__)
 #else
-#define MRERROR_HANDLING_LOC_ ::mr::error_handling::detail::CodeLocation(__FILE__, __LINE__)
+#define MRUTIL_ERROR_HANDLING_LOC_ ::mr::detail_error_handling::CodeLocation(__FILE__, __LINE__)
 #endif
 
 #define MR_fail(...) \
   do { \
-    if (!::mr::error_handling::detail::abort_in_progress__) \
+    if (!::mr::detail_error_handling::abort_in_progress__) \
       { \
-      ::mr::error_handling::detail::abort_in_progress__ = true; \
-      ::mr::error_handling::detail::streamDump__(::std::cerr, MR_ERROR_HANDLING_LOC_, "\n", ##__VA_ARGS__, "\n"); \
-      ::mr::error_handling::detail::killjob__(); \
+      ::mr::detail_error_handling::abort_in_progress__ = true; \
+      ::mr::detail_error_handling::streamDump__(::std::cerr, MRUTIL_ERROR_HANDLING_LOC_, "\n", ##__VA_ARGS__, "\n"); \
+      ::mr::detail_error_handling::killjob__(); \
       } \
     ::std::exit(1); \
     } while(0)
@@ -92,6 +90,6 @@ inline void streamDump__(::std::ostream &os, const T& value,
   }
 #endif
 
-}}}
+}}
 
 #endif

src/mr_util/gl_integrator.h

@@ -28,12 +28,9 @@
 
 namespace mr {
 
-namespace gl_integrator {
-
-namespace detail {
+namespace detail_gl_integrator {
 
 using namespace std;
-using namespace mr::threading;
 
 class GL_Integrator
   {
@@ -125,8 +122,8 @@ class GL_Integrator
 
 }
 
-using detail::GL_Integrator;
+using detail_gl_integrator::GL_Integrator;
 
-}}
+}
 
 #endif

src/mr_util/simd.h

@@ -61,7 +61,7 @@ constexpr size_t vbytes = 16;
 template<typename T, size_t len=vbytes/sizeof(T)> class vtp
   {
   protected:
-    using Tv __attribute__ ((vector_size (len*sizeof(T)))) = T;
+    using Tv [[gnu::vector_size (len*sizeof(T))]] = T;
     static_assert((len>0) && ((len&(len-1))==0), "bad vector length");
     Tv v;
 
@@ -104,7 +104,6 @@ template<typename T, size_t len=vbytes/sizeof(T)> class vtp
     template<typename I> void Set (I i, T val) { v[i]=val; }
     template<typename I> T operator[](I i) const { return v[i]; }
   };
-
 }
 
 using detail_simd::vtp;

src/mr_util/string_utils.h

@@ -27,8 +27,6 @@
 
 namespace mr {
 
-namespace string_utils {
-
 /*! \defgroup stringutilsgroup String handling helper functions */
 /*! \{ */
 
@@ -113,6 +111,6 @@ void parse_words_from_file (const std::string &filename,
 
 /*! \} */
 
-}}
+}
 
 #endif

src/mr_util/system.h

@@ -27,12 +27,10 @@
 
 namespace mr {
 
-namespace system {
-
 std::size_t getProcessInfo(const std::string &quantity);
 std::size_t getMemInfo(const std::string &quantity);
 std::size_t usable_memory();
 
-}}
+}
 
 #endif

src/mr_util/threading.h

@@ -38,9 +38,7 @@
 
 namespace mr {
 
-namespace threading {
-
-namespace detail {
+namespace detail_threading {
 
 using namespace std;
 
@@ -346,15 +344,14 @@ template<typename Func> void execGuided(size_t nwork,
   }
 } // end of namespace detail
 
-using detail::Scheduler;
-using detail::execSingle;
-using detail::execStatic;
-using detail::execDynamic;
-using detail::execGuided;
+using detail_threading::Scheduler;
+using detail_threading::execSingle;
+using detail_threading::execStatic;
+using detail_threading::execDynamic;
+using detail_threading::execGuided;
 
 // FIXME: missing execParallel(), my_thread(), num_threads()
 
-} // end of namespace threading
 } // end of namespace mr
 
 #endif

src/mr_util/timers.h

@@ -6,7 +6,6 @@
 #include <map>
 
 #include "mr_util/error_handling.h"
-#include "mr_util/string_utils.h"
 
 namespace mr {
 

src/mr_util/unity_roots.h

+/*
+ *  This file is part of the MR utility library.
+ *
+ *  This code 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.
+ *
+ *  This code 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 this code; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+/* Copyright (C) 2019 Max-Planck-Society
+   Author: Martin Reinecke */
+
+#ifndef MRUTIL_UNITY_ROOTS_H
+#define MRUTIL_UNITY_ROOTS_H
+
+#include <complex>
+#include <cmath>
+#include <type_traits>
+#include <vector>
+#include <mr_util/cmplx.h>
+
+namespace mr {
+
+namespace detail_unity_roots {
+
+using namespace std;
+
+template<typename T, typename Tc=complex<T>> class UnityRoots
+  {
+  private:
+    using Thigh = typename conditional<(sizeof(T)>sizeof(double)), T, double>::type;
+    size_t N, mask, shift;
+    vector<Cmplx<Thigh>> v1, v2;
+
+    static Cmplx<Thigh> calc(size_t x, size_t n, Thigh ang)
+      {
+      x<<=3;
+      if (x<4*n) // first half
+        {
+        if (x<2*n) // first quadrant
+          {
+          if (x<n) return Cmplx<Thigh>(cos(Thigh(x)*ang), sin(Thigh(x)*ang));
+          return Cmplx<Thigh>(sin(Thigh(2*n-x)*ang), cos(Thigh(2*n-x)*ang));
+          }
+        else // second quadrant
+          {
+          x-=2*n;
+          if (x<n) return Cmplx<Thigh>(-sin(Thigh(x)*ang), cos(Thigh(x)*ang));
+          return Cmplx<Thigh>(-cos(Thigh(2*n-x)*ang), sin(Thigh(2*n-x)*ang));
+          }
+        }
+      else
+        {
+        x=8*n-x;
+        if (x<2*n) // third quadrant
+          {
+          if (x<n) return Cmplx<Thigh>(cos(Thigh(x)*ang), -sin(Thigh(x)*ang));
+          return Cmplx<Thigh>(sin(Thigh(2*n-x)*ang), -cos(Thigh(2*n-x)*ang));
+          }
+        else // fourth quadrant
+          {
+          x-=2*n;
+          if (x<n) return Cmplx<Thigh>(-sin(Thigh(x)*ang), -cos(Thigh(x)*ang));
+          return Cmplx<Thigh>(-cos(Thigh(2*n-x)*ang), -sin(Thigh(2*n-x)*ang));
+          }
+        }
+      }
+
+  public:
+    UnityRoots(size_t n)
+      : N(n)
+      {
+      constexpr auto pi = 3.141592653589793238462643383279502884197L;
+      Thigh ang = Thigh(0.25L*pi/n);
+      size_t nval = (n+2)/2;
+      shift = 1;
+      while((size_t(1)<<shift)*(size_t(1)<<shift) < nval) ++shift;
+      mask = (size_t(1)<<shift)-1;
+      v1.resize(mask+1);
+      v1[0].Set(Thigh(1), Thigh(0));
+      for (size_t i=1; i<v1.size(); ++i)
+        v1[i]=calc(i,n,ang);
+      v2.resize((nval+mask)/(mask+1));
+      v2[0].Set(Thigh(1), Thigh(0));
+      for (size_t i=1; i<v2.size(); ++i)
+        v2[i]=calc(i*(mask+1),n,ang);
+      }
+
+    Tc operator[](size_t idx) const
+      {
+      if (2*idx<=N)
+        {
+        auto x1=v1[idx&mask], x2=v2[idx>>shift];
+        return Tc(T(x1.r*x2.r-x1.i*x2.i), T(x1.r*x2.i+x1.i*x2.r));
+        }
+      idx = N-idx;
+      auto x1=v1[idx&mask], x2=v2[idx>>shift];
+      return Tc(T(x1.r*x2.r-x1.i*x2.i), -T(x1.r*x2.i+x1.i*x2.r));
+      }
+  };
+
+}
+
+using detail_unity_roots::UnityRoots;
+
+}
+
+#endif

src/string_utils.cc

@@ -26,7 +26,7 @@
  *  This file contains the implementation of various convenience functions
  *  used by the Planck LevelS package.
  *
- *  Copyright (C) 2002-2014 Max-Planck-Society
+ *  Copyright (C) 2002-2019 Max-Planck-Society
  *  Author: Martin Reinecke
  */
 
@@ -44,8 +44,6 @@ using namespace std;
 
 namespace mr {
 
-namespace string_utils {
-
 string trim (const string &orig)
   {
   string::size_type p1=orig.find_first_not_of(" \t");
@@ -332,4 +330,4 @@ void parse_words_from_file (const string &filename, vector<string> &words)
     }
   }
 
-}}
+}

src/system.cc

+/*
+ *  This file is part of the MR utility library.
+ *
+ *  This code 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.
+ *
+ *  This code 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 this code; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+/* Copyright (C) 2019 Max-Planck-Society
+   Author: Martin Reinecke */
+
 #include <regex>
 #include <string>
 #include <fstream>
@@ -8,12 +29,9 @@
 #include "mr_util/string_utils.h"
 
 using namespace std;
-using namespace mr::string_utils;
 
 namespace mr {
 
-namespace system {
-
 namespace {
 
 string fileToString(const string &fname)
@@ -55,4 +73,4 @@ size_t usable_memory()
   return MemTotal-Committed;
   }
 
-}}
+}