[Boost-commit] svn:boost r67022 - in sandbox/numeric_bindings: boost/numeric/bindings/detail boost/numeric/bindings/traits/detail libs/numeric/bindings/lapack/test

Subject: [Boost-commit] svn:boost r67022 - in sandbox/numeric_bindings: boost/numeric/bindings/detail boost/numeric/bindings/traits/detail libs/numeric/bindings/lapack/test
From: thomas.klimpel_at_[hidden]
Date: 2010-12-05 12:30:29

Author: klimpel
Date: 2010-12-05 12:30:25 EST (Sun, 05 Dec 2010)
New Revision: 67022
URL: http://svn.boost.org/trac/boost/changeset/67022

Log:
Polished interleave functionality for complex vectors a bit, such that the usage is simplified and no reinterpret_casts are needed by the code using this new functionality.
Added:
   sandbox/numeric_bindings/boost/numeric/bindings/detail/complex_utils.hpp   (contents, props changed)
Text files modified: 
   sandbox/numeric_bindings/boost/numeric/bindings/traits/detail/utils.hpp   |    95 ----------------------------------------
   sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/hseqr.cpp      |    18 ++----                                  
   sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_gees.cpp |    10 +--                                     
   3 files changed, 10 insertions(+), 113 deletions(-)

Added: sandbox/numeric_bindings/boost/numeric/bindings/detail/complex_utils.hpp
==============================================================================
--- (empty file)
+++ sandbox/numeric_bindings/boost/numeric/bindings/detail/complex_utils.hpp	2010-12-05 12:30:25 EST (Sun, 05 Dec 2010)
@@ -0,0 +1,178 @@
+//
+// Copyright (c) 2003 Kresimir Fresl
+// Copyright (c) 2010 Thomas Klimpel
+//
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+//
+
+#ifndef BOOST_NUMERIC_BINDINGS_DETAIL_COMPLEX_UTILS_HPP
+#define BOOST_NUMERIC_BINDINGS_DETAIL_COMPLEX_UTILS_HPP
+
+#include <iterator>
+#include <boost/numeric/bindings/begin.hpp>
+#include <boost/numeric/bindings/end.hpp>
+#include <boost/numeric/bindings/is_complex.hpp>
+#include <boost/numeric/bindings/remove_imaginary.hpp>
+#include <boost/numeric/bindings/value_type.hpp>
+#include <boost/numeric/bindings/vector_view.hpp>
+#include <boost/utility/enable_if.hpp>
+
+namespace boost {
+namespace numeric {
+namespace bindings {
+
+namespace detail {
+
+#ifdef BOOST_NUMERIC_BINDINGS_BY_THE_BOOK
+template <typename It>
+void inshuffle(It it, std::size_t n) {
+  if (n==0) return;
+  for (std::size_t i = 0; 2*i < n; ++i) {
+    std::size_t k = 2*i + 1;
+    while (2*k <= n) k *= 2;
+    typename std::iterator_traits<It>::value_type tmp = it[n+i];
+    it[n+i] = it[k-1];
+    while (k % 2 == 0) {
+      it[k-1] = it[(k/2)-1];
+      k /= 2;
+    }
+    it[k-1] = tmp;
+  }
+  std::size_t kmin = 1;
+  while (2*kmin <= n) kmin *= 2;
+  for (std::size_t i = 0; 4*i+1 < n; ++i) {
+    std::size_t k = 2*i + 1;
+    while (2*k <= n) k *= 2;
+    std::size_t k1 = 2*(i+1) + 1;
+    while (2*k1 <= n) k1 *= 2;
+    if (k > k1) {
+      if (k1 < kmin) {
+        kmin = k1;
+        inshuffle(it+n, i+1);
+      }
+      else
+        inshuffle(it+n+1, i);
+    }
+  }
+  return inshuffle(it+n+(n%2), n/2);
+}
+#else
+template <typename It>
+void inshuffle(It it, std::size_t n) {
+  while (n > 0) {
+    std::size_t kmin = 1;
+    while (kmin <= n)
+      kmin *= 2;
+    {
+      std::size_t kk = kmin/2;
+      It itn = it + n;
+      for (std::size_t i = 0, s = (n+1)/2; i < s; ++i) {
+        std::size_t k = (2*i+1)*kk;
+        while (k > n) {
+          k /= 2;
+          kk /= 2;
+        }
+        // apply the cyclic permutation
+        typename std::iterator_traits<It>::value_type tmp = itn[i];
+        itn[i] = it[k-1];
+        while (k % 2 == 0) {
+          it[k-1] = it[(k/2)-1];
+          k /= 2;
+        }
+        it[k-1] = tmp;
+      }
+    }
+    // the optimized computation of k fails for n=2,
+    // so skip the 'normalization' loop when possible
+    if (n > 3) {
+      std::size_t kk = kmin/4;
+      for (std::size_t i = 1; 4*i < n+3; ++i) {
+        std::size_t k = (2*i+1)*kk;
+        if (k > n) {
+          kk /= 2;
+          if (k < kmin) {
+            kmin = k;
+            // if kmin is updated, do an in-shuffle
+            inshuffle(it+n, i);
+          }
+          else
+            // otherwise do an out-shuffle
+            inshuffle(it+n+1, i-1);
+        }
+      }
+    }
+    // implement the tail recursion as an iteration
+    it += n+(n%2);
+    n /= 2;
+  }
+}
+#endif
+
+// Reorders a real array followed by an imaginary array to a true complex array
+// where real and imaginary part of each number directly follow each other.
+template <typename VectorW>
+typename boost::enable_if< is_complex< typename bindings::value_type< VectorW >::type >, void >::type
+interlace (VectorW& w) {
+  typedef typename bindings::value_type< VectorW >::type value_type;
+  typedef typename bindings::remove_imaginary< value_type >::type real_type;
+  value_type* pw = bindings::begin_value(w);
+  std::ptrdiff_t n = bindings::end_value(w) - pw;
+  if (n < 2) return;
+  inshuffle(reinterpret_cast<real_type*> (pw)+1, n-1);
+}
+
+namespace result_of {
+
+template< typename VectorW >
+struct real_part_view {
+    typedef typename bindings::result_of::vector_view< typename
+      bindings::remove_imaginary< typename
+      bindings::value_type< VectorW >::type
+      >::type >::type type;
+};
+
+template< typename VectorW >
+struct imag_part_view {
+    typedef typename bindings::result_of::vector_view< typename
+      bindings::remove_imaginary< typename
+      bindings::value_type< VectorW >::type
+      >::type >::type type;
+};
+
+} // namespace result_of
+
+// Creates a real vector_view to the first half of the complex array,
+// which is intended to be filled by the real part
+template <typename VectorW>
+typename boost::enable_if< is_complex< typename bindings::value_type< VectorW >::type >,
+        typename result_of::real_part_view< VectorW >::type const >::type
+real_part_view (VectorW& w) {
+  typedef typename bindings::value_type< VectorW >::type value_type;
+  typedef typename bindings::remove_imaginary< value_type >::type real_type;
+  value_type* pw = bindings::begin_value(w);
+  std::ptrdiff_t n = bindings::end_value(w) - pw;
+  return bindings::vector_view(reinterpret_cast<real_type*> (pw), n);
+}
+
+// Creates a real vector_view to the second half of the complex array,
+// which is intended to be filled by the imaginary part
+template <typename VectorW>
+typename boost::enable_if< is_complex< typename bindings::value_type< VectorW >::type >,
+        typename result_of::imag_part_view< VectorW >::type const >::type
+imag_part_view (VectorW& w) {
+  typedef typename bindings::value_type< VectorW >::type value_type;
+  typedef typename bindings::remove_imaginary< value_type >::type real_type;
+  value_type* pw = bindings::begin_value(w);
+  std::ptrdiff_t n = bindings::end_value(w) - pw;
+  return bindings::vector_view(reinterpret_cast<real_type*> (pw)+n, n);
+}
+
+} // namespace detail
+
+} // namespace bindings
+} // namespace numeric
+} // namespace boost
+
+#endif

Modified: sandbox/numeric_bindings/boost/numeric/bindings/traits/detail/utils.hpp
==============================================================================
--- sandbox/numeric_bindings/boost/numeric/bindings/traits/detail/utils.hpp	(original)
+++ sandbox/numeric_bindings/boost/numeric/bindings/traits/detail/utils.hpp	2010-12-05 12:30:25 EST (Sun, 05 Dec 2010)
@@ -1,7 +1,6 @@
 /*
  *
  * Copyright (c) Kresimir Fresl 2003
- * Copyright (c) 2010 Thomas Klimpel
  *
  * Distributed under the Boost Software License, Version 1.0.
  * (See accompanying file LICENSE_1_0.txt or copy at
@@ -47,100 +46,6 @@
 #endif
     }
 
-#ifdef BOOST_NUMERIC_BINDINGS_BY_THE_BOOK
-    template <typename It>
-    void inshuffle(It it, std::size_t n) {
-      if (n==0) return;
-      for (std::size_t i = 0; 2*i < n; ++i) {
-        std::size_t k = 2*i + 1;
-        while (2*k <= n) k *= 2;
-        typename std::iterator_traits<It>::value_type tmp = it[n+i];
-        it[n+i] = it[k-1];
-        while (k % 2 == 0) {
-          it[k-1] = it[(k/2)-1];
-          k /= 2;
-        }
-        it[k-1] = tmp;
-      }
-      std::size_t kmin = 1;
-      while (2*kmin <= n) kmin *= 2;
-      for (std::size_t i = 0; 4*i+1 < n; ++i) {
-        std::size_t k = 2*i + 1;
-        while (2*k <= n) k *= 2;
-        std::size_t k1 = 2*(i+1) + 1;
-        while (2*k1 <= n) k1 *= 2;
-        if (k > k1) {
-          if (k1 < kmin) {
-            kmin = k1;
-            inshuffle(it+n, i+1);
-          }
-          else
-            inshuffle(it+n+1, i);
-        }
-      }
-      return inshuffle(it+n+(n%2), n/2);
-    }
-#else
-    template <typename It>
-    void inshuffle(It it, std::size_t n) {
-      while (n > 0) {
-        std::size_t kmin = 1;
-        while (kmin <= n)
-          kmin *= 2;
-        {
-          std::size_t kk = kmin/2;
-          It itn = it + n;
-          for (std::size_t i = 0, s = (n+1)/2; i < s; ++i) {
-            std::size_t k = (2*i+1)*kk;
-            while (k > n) {
-              k /= 2;
-              kk /= 2;
-            }
-            // apply the cyclic permutation
-            typename std::iterator_traits<It>::value_type tmp = itn[i];
-            itn[i] = it[k-1];
-            while (k % 2 == 0) {
-              it[k-1] = it[(k/2)-1];
-              k /= 2;
-            }
-            it[k-1] = tmp;
-          }
-        }
-        // the optimized computation of k fails for n=2,
-        // so skip the 'normalization' loop when possible
-        if (n > 3) {
-          std::size_t kk = kmin/4;
-          for (std::size_t i = 1; 4*i < n+3; ++i) {
-            std::size_t k = (2*i+1)*kk;
-            if (k > n) {
-              kk /= 2;
-              if (k < kmin) {
-                kmin = k;
-                // if kmin is updated, do an in-shuffle
-                inshuffle(it+n, i);
-              }
-              else
-                // otherwise do an out-shuffle
-                inshuffle(it+n+1, i-1);
-            }
-          }
-        }
-        // implement the tail recursion as an iteration
-        it += n+(n%2);
-        n /= 2;
-      }
-    }
-#endif
-
-    // real & imaginary arrays => complex array
-    template <typename CIt>
-    void interlace (CIt c, std::ptrdiff_t n) {
-      typedef typename std::iterator_traits<CIt>::value_type cmplx_t;
-      typedef typename type_traits<cmplx_t>::real_type real_t;
-      if (n < 2) return;
-      inshuffle(reinterpret_cast<real_t*> (&*c)+1, n-1);
-    }
-
     // converts real/complex to std::ptrdiff_t
     inline std::ptrdiff_t to_int (float f) { return static_cast<std::ptrdiff_t> (f); }
     inline std::ptrdiff_t to_int (double d) { return static_cast<std::ptrdiff_t> (d); }

Modified: sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/hseqr.cpp
==============================================================================
--- sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/hseqr.cpp	(original)
+++ sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/hseqr.cpp	2010-12-05 12:30:25 EST (Sun, 05 Dec 2010)
@@ -8,7 +8,7 @@
 #include <boost/numeric/ublas/io.hpp>
 #include <boost/numeric/bindings/ublas/matrix.hpp>
 #include <boost/numeric/bindings/ublas/vector.hpp>
-#include <boost/numeric/bindings/traits/detail/utils.hpp>
+#include <boost/numeric/bindings/detail/complex_utils.hpp>
 #include <boost/numeric/bindings/vector_view.hpp>
 #include <boost/numeric/bindings/lapack/computational/hseqr.hpp>
 #include <boost/numeric/bindings/lapack/computational/trevc.hpp>
@@ -49,12 +49,10 @@
     cout << "\nHSEQR for only eigenvalues." << endl;
     ublas::matrix<double, ublas::column_major> Z_dummy(1,1);
     lapack::hseqr('E', 'N', 1, n, H,
-        bindings::vector_view(reinterpret_cast<double*>(
-            &*bindings::begin_value(values)), bindings::size(values)),
-        bindings::vector_view(reinterpret_cast<double*>(
-            &*bindings::begin_value(values))+bindings::size(values), bindings::size(values)),
+        bindings::detail::real_part_view(values),
+        bindings::detail::imag_part_view(values),
         Z_dummy);
-    bindings::traits::detail::interlace(bindings::begin_value(values), bindings::size(values));
+    bindings::detail::interlace(values);
     cout << "\nH:\n" << H << endl;
     cout << "\nvalues: " << values << endl;
 
@@ -62,12 +60,10 @@
     Hessenberg(H);
     cout << "H:\n" << H << endl;
     lapack::hseqr('S', 'I', 1, n, H,
-        bindings::vector_view(reinterpret_cast<double*>(
-            &*bindings::begin_value(values)), bindings::size(values)),
-        bindings::vector_view(reinterpret_cast<double*>(
-            &*bindings::begin_value(values))+bindings::size(values), bindings::size(values)),
+        bindings::detail::real_part_view(values),
+        bindings::detail::imag_part_view(values),
         Z);
-    bindings::traits::detail::interlace(bindings::begin_value(values), bindings::size(values));
+    bindings::detail::interlace(values);
     cout << "\nH: " << H << endl;
     cout << "\nvalues: " << values << endl;
     cout << "\nZ: " << Z << endl;

Modified: sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_gees.cpp
==============================================================================
--- sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_gees.cpp	(original)
+++ sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_gees.cpp	2010-12-05 12:30:25 EST (Sun, 05 Dec 2010)
@@ -16,6 +16,7 @@
 #include <boost/numeric/bindings/lapack/driver/gees.hpp>
 #include <boost/numeric/bindings/detail/array.hpp>
 #include <boost/numeric/bindings/vector_view.hpp>
+#include <boost/numeric/bindings/detail/complex_utils.hpp>
 #include <boost/numeric/ublas/io.hpp>
 #include <boost/type_traits/is_complex.hpp>
 #include <boost/mpl/if.hpp>
@@ -26,7 +27,6 @@
 
 namespace ublas = boost::numeric::ublas;
 namespace lapack = boost::numeric::bindings::lapack;
-namespace traits = boost::numeric::bindings::traits;
 namespace bindings = boost::numeric::bindings;
 
 struct apply_real {
@@ -35,14 +35,10 @@
   static inline std::ptrdiff_t gees( const char jobvs, const char sort,
         external_fp select, MatrixA& a, fortran_int_t& sdim, VectorW& w,
         MatrixVS& vs, Workspace work ) {
-    typedef typename bindings::value_type< MatrixA >::type value_type;
     fortran_int_t info = lapack::gees( jobvs, sort, select, a, sdim,
-      bindings::vector_view(reinterpret_cast<value_type*>(
-        &*bindings::begin_value(w)), bindings::size(w)),
-      bindings::vector_view(reinterpret_cast<value_type*>(
-        &*bindings::begin_value(w))+bindings::size(w), bindings::size(w)),
+      bindings::detail::real_part_view(w), bindings::detail::imag_part_view(w),
       vs, work );
-    traits::detail::interlace(bindings::begin_value(w), bindings::size(w));
+    bindings::detail::interlace(w);
     return info;
   }
 };