Merge pull request #14 from bjodah/calculate-weights-optim

Introduced finitediff::calculate_weights_optim
bjodah · Aug 24, 2016 · 3d3111c · 3d3111c
2 parents b555957 + b00b241
commit 3d3111c
Show file tree

Hide file tree

Showing 3 changed files with 48 additions and 8 deletions.
diff --git a/CHANGES.rst b/CHANGES.rst
@@ -1,7 +1,9 @@
 v0.3.1
 ======
 - Introduced ``finitediff::generate_weights_optim`` which is a pre-sorting
-  wrapper around ``finitediff::generate_weights``
+  wrapper around ``finitediff::generate_weights``  (tests indicate approx. 1 extra significant figure)
+- Introduced ``finitediff::calculate_weights_optim`` which is a pre-sorting
+  wrapper around ``finitediff::calculate_weights`` (tests indicate approx. 1 extra significant figure)
 
 v0.3.0
 ======

diff --git a/include/finitediff_templated.hpp b/include/finitediff_templated.hpp
@@ -5,13 +5,13 @@
 #if __cplusplus > 199711L
 #include <algorithm>
 #include <stdexcept>
+#include <vector>
 #endif
-
 namespace finitediff {
 
     template <typename Real_t>
-    void calculate_weights(const Real_t * const __restrict__ grid, const int len_g,
-                           const int max_deriv, Real_t * const __restrict__ weights, const Real_t around=0) {
+    void calculate_weights(const Real_t * const __restrict__ grid, const unsigned len_g,
+                           const unsigned max_deriv, Real_t * const __restrict__ weights, const Real_t around=0) {
         // Parameters
         // ----------
         // grid[len_g]: array with grid point locations
@@ -62,7 +62,6 @@ namespace finitediff {
         }
     }
 
-
     // populate_weights is deprecated due to counter-intuitive parameter "nd"
     template <typename Real_t>
     void populate_weights(const Real_t z, const Real_t * const __restrict__ x, const int nd,
@@ -147,8 +146,6 @@ namespace finitediff {
         return coeffs;
     }
 
-    template<class T> static inline T abs_(T val) { return (val < 0) ? -val : val; }
-
     template<typename Real_t, template<typename, typename...> class Cont, typename... Args>
     Cont<Real_t, Args...> generate_weights_optim(const Cont<Real_t, Args...>& grid, int maxorder=-1, const Real_t around=0){
         const unsigned n_ = grid.size();
@@ -157,7 +154,7 @@ namespace finitediff {
             index[idx] = idx;
         std::sort(index.begin(), index.end(),
              [&](const unsigned& a, const unsigned& b) {
-                 return (abs_(grid[a] - around) < abs_(grid[b] - around));
+                      return (std::abs(grid[a] - around) < std::abs(grid[b] - around));
              });
         Cont<Real_t, Args...> reordered_grid(n_);
         for (unsigned idx=0; idx < n_; ++idx){
@@ -174,6 +171,33 @@ namespace finitediff {
         }
         return coeffs;
     }
+
+    template <typename Real_t>
+    void calculate_weights_optim(const Real_t * const __restrict__ grid, const unsigned len_g,
+                                 const unsigned max_deriv, Real_t * const __restrict__ weights, const Real_t around=0) {
+        if (len_g < max_deriv + 1){
+            throw std::logic_error("size of grid insufficient");
+        }
+        std::vector<unsigned> index(len_g);
+        for (unsigned idx=0; idx < len_g; ++idx)
+            index[idx] = idx;
+        std::sort(index.begin(), index.end(),
+             [&](const unsigned& a, const unsigned& b) {
+                      return (std::abs(grid[a] - around) < std::abs(grid[b] - around));
+             });
+        Real_t * const __restrict__ reordered_grid = new Real_t[len_g];
+        for (unsigned idx=0; idx < len_g; ++idx){
+            reordered_grid[idx] = grid[index[idx]];
+        }
+        Real_t * const __restrict__ reordered_weights = new Real_t[len_g*(max_deriv + 1)];
+        calculate_weights(reordered_grid, len_g, max_deriv, reordered_weights, around);
+        for (unsigned order=0; order <= max_deriv; ++order){
+            for (unsigned idx=0; idx < len_g; ++idx)
+                weights[order*len_g + index[idx]] += reordered_weights[order*len_g + idx];
+        }
+        delete []reordered_weights;
+        delete []reordered_grid;
+    }
 #endif
 }
 

diff --git a/tests/test_finitediff_templated.cpp b/tests/test_finitediff_templated.cpp
@@ -102,4 +102,18 @@ TEST_CASE( "shifted", "finitediff::generate_weights_optim") {
         }
     }
 
+    std::vector<double> calc_naive(n*(n+3)/2);
+    std::vector<double> calc_optim(n*(n+3)/2);
+    finitediff::calculate_weights(&x_naive[0], n, (n+1)/2, &calc_naive[0], 100.0);
+    finitediff::calculate_weights_optim(&x_naive[0], n, (n+1)/2, &calc_optim[0], 100.0);
+    for (unsigned deriv_i = 0; deriv_i < ncols; deriv_i++){
+        for (unsigned idx = 0; idx < n; idx++){
+            const unsigned idx_naive = deriv_i*n + idx;
+            const unsigned idx_refer = deriv_i*n + mapping[idx];
+            const double absdiff_naive = std::abs(calc_naive[idx_naive] - c_refer[idx_refer]);
+            const double absdiff_optim = std::abs(calc_optim[idx_naive] - c_refer[idx_refer]);
+            REQUIRE( absdiff_naive*1e15 < 1 );
+            REQUIRE( absdiff_optim*1e16 < 1 );
+        }
+    }
 }