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

Subject: [Boost-commit] svn:boost r66669 - in sandbox/numeric_bindings: boost/numeric/bindings boost/numeric/bindings/traits/detail libs/numeric/bindings libs/numeric/bindings/lapack/test libs/numeric/bindings/test
From: thomas.klimpel_at_[hidden]
Date: 2010-11-21 19:42:11

Author: klimpel
Date: 2010-11-21 19:42:07 EST (Sun, 21 Nov 2010)
New Revision: 66669
URL: http://svn.boost.org/trac/boost/changeset/66669

Log:
implemented an inplace version of interlace
added a view named vector_view to turn a pointer+size into a vector
used new functionality to replace unnecessary memory allocation in the hseqr.cpp and ublas_gees.cpp tests.
Added:
   sandbox/numeric_bindings/boost/numeric/bindings/vector_view.hpp   (contents, props changed)
   sandbox/numeric_bindings/libs/numeric/bindings/test/Jamfile.v2   (contents, props changed)
   sandbox/numeric_bindings/libs/numeric/bindings/test/interlace.cpp   (contents, props changed)
Text files modified: 
   sandbox/numeric_bindings/boost/numeric/bindings/traits/detail/utils.hpp   |   119 +++++++++++++++++++++++++++++++++++---- 
   sandbox/numeric_bindings/boost/numeric/bindings/views.hpp                 |     1                                         
   sandbox/numeric_bindings/libs/numeric/bindings/Jamfile.v2                 |     1                                         
   sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/hseqr.cpp      |    23 +++++--                                 
   sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_gees.cpp |    15 ++--                                    
   5 files changed, 132 insertions(+), 27 deletions(-)

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-11-21 19:42:07 EST (Sun, 21 Nov 2010)
@@ -1,12 +1,13 @@
 /*
- * 
+ *
  * 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
  * http://www.boost.org/LICENSE_1_0.txt)
  *
- * Author acknowledges the support of the Faculty of Civil Engineering, 
+ * Author acknowledges the support of the Faculty of Civil Engineering,
  * University of Zagreb, Croatia.
  *
  */
@@ -14,7 +15,6 @@
 #ifndef BOOST_NUMERIC_BINDINGS_TRAITS_DETAIL_UTILS_HPP
 #define BOOST_NUMERIC_BINDINGS_TRAITS_DETAIL_UTILS_HPP
 
-// #include <cstring> 
 #include <iterator>
 #include <boost/numeric/bindings/traits/type_traits.hpp>
 
@@ -35,31 +35,124 @@
     void interlace (RIt r, RIt r_end, RIt ri, CIt c) {
       typedef typename std::iterator_traits<CIt>::value_type cmplx_t;
 #ifdef BOOST_NUMERIC_BINDINGS_BY_THE_BOOK
-      for (; r != r_end; ++r, ++ri, ++c) 
-        *c = cmplx_t (*r, *ri); 
+      for (; r != r_end; ++r, ++ri, ++c)
+        *c = cmplx_t (*r, *ri);
 #else
-      typedef typename type_traits<cmplx_t>::real_type real_t; 
+      typedef typename type_traits<cmplx_t>::real_type real_t;
       real_t *cp = reinterpret_cast<real_t*> (&*c);
       for (; r != r_end; ++r, ++ri) {
         *cp = *r; ++cp;
         *cp = *ri; ++cp;
       }
-#endif 
-    }    
+#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); }
-    inline std::ptrdiff_t to_int (traits::complex_f const& cf) { 
-      return static_cast<std::ptrdiff_t> (traits::real (cf)); 
+    inline std::ptrdiff_t to_int (traits::complex_f const& cf) {
+      return static_cast<std::ptrdiff_t> (traits::real (cf));
     }
-    inline std::ptrdiff_t to_int (traits::complex_d const& cd) { 
-      return static_cast<std::ptrdiff_t> (traits::real (cd)); 
+    inline std::ptrdiff_t to_int (traits::complex_d const& cd) {
+      return static_cast<std::ptrdiff_t> (traits::real (cd));
     }
 
   }
 
 }}}}
 
-#endif 
+#endif

Added: sandbox/numeric_bindings/boost/numeric/bindings/vector_view.hpp
==============================================================================
--- (empty file)
+++ sandbox/numeric_bindings/boost/numeric/bindings/vector_view.hpp	2010-11-21 19:42:07 EST (Sun, 21 Nov 2010)
@@ -0,0 +1,83 @@
+//
+// Copyright (c) 2009 Rutger ter Borg
+// 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_VECTOR_VIEW_HPP
+#define BOOST_NUMERIC_BINDINGS_VECTOR_VIEW_HPP
+
+#include <boost/numeric/bindings/detail/adaptable_type.hpp>
+
+namespace boost {
+namespace numeric {
+namespace bindings {
+namespace detail {
+
+template< typename T >
+struct vector_view_wrapper:
+        adaptable_type< vector_view_wrapper<T> > {
+    typedef T value_type;
+
+    vector_view_wrapper( T* t, std::size_t size ):
+        m_t( t ),
+        m_size( size ) {}
+    T* m_t;
+    std::size_t m_size;
+};
+
+template< typename T, typename Id, typename Enable >
+struct adaptor< vector_view_wrapper<T>, Id, Enable > {
+
+    typedef typename Id::value_type value_type;
+    typedef mpl::map<
+        mpl::pair< tag::value_type, value_type >,
+        mpl::pair< tag::entity, tag::vector >,
+        mpl::pair< tag::size_type<1>, std::ptrdiff_t >,
+        mpl::pair< tag::data_structure, tag::linear_array >,
+        mpl::pair< tag::stride_type<1>, tag::contiguous >
+    > property_map;
+
+    static std::ptrdiff_t size1( const Id& id ) {
+        return id.m_size;
+    }
+
+    static value_type* begin_value( Id& id ) {
+        return id.m_t;
+    }
+
+    static value_type* end_value( Id& id ) {
+        return id.m_t + id.m_size;
+    }
+
+};
+
+} // namespace detail
+
+namespace result_of {
+
+template< typename T >
+struct vector_view {
+    typedef detail::vector_view_wrapper<T> type;
+};
+
+} // namespace result_of
+
+template< typename T >
+detail::vector_view_wrapper<T> const vector_view( T* t, std::size_t size ) {
+    return detail::vector_view_wrapper<T>( t, size );
+}
+
+template< typename T >
+detail::vector_view_wrapper<const T> const vector_view( const T* t, std::size_t size ) {
+    return detail::vector_view_wrapper<const T>( t, size );
+}
+
+} // namespace bindings
+} // namespace numeric
+} // namespace boost
+
+#endif

Modified: sandbox/numeric_bindings/boost/numeric/bindings/views.hpp
==============================================================================
--- sandbox/numeric_bindings/boost/numeric/bindings/views.hpp	(original)
+++ sandbox/numeric_bindings/boost/numeric/bindings/views.hpp	2010-11-21 19:42:07 EST (Sun, 21 Nov 2010)
@@ -20,6 +20,7 @@
 #include <boost/numeric/bindings/unit_lower.hpp>
 #include <boost/numeric/bindings/unit_upper.hpp>
 #include <boost/numeric/bindings/upper.hpp>
+#include <boost/numeric/bindings/vector_view.hpp>
 
 #endif
 

Modified: sandbox/numeric_bindings/libs/numeric/bindings/Jamfile.v2
==============================================================================
--- sandbox/numeric_bindings/libs/numeric/bindings/Jamfile.v2	(original)
+++ sandbox/numeric_bindings/libs/numeric/bindings/Jamfile.v2	2010-11-21 19:42:07 EST (Sun, 21 Nov 2010)
@@ -11,4 +11,5 @@
 build-project blas/test ;
 build-project lapack/test ;
 build-project mumps/test ;
+build-project test ;
 build-project umfpack/test ;

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-11-21 19:42:07 EST (Sun, 21 Nov 2010)
@@ -9,6 +9,7 @@
 #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/vector_view.hpp>
 #include <boost/numeric/bindings/lapack/computational/hseqr.hpp>
 #include <boost/numeric/bindings/lapack/computational/trevc.hpp>
 
@@ -19,7 +20,7 @@
 
 namespace ublas =  boost::numeric::ublas;
 namespace lapack =  boost::numeric::bindings::lapack;
-namespace traits =  boost::numeric::bindings::traits;
+namespace bindings =  boost::numeric::bindings;
 namespace tag =  boost::numeric::bindings::tag;
 
 void hseqr(int);
@@ -43,22 +44,30 @@
     cout << "\nUpper Hessenberg matrix H:\n" << H << endl;
 
     ublas::vector<complex<double> > values(n);
-    ublas::vector<double> values_r(n);
-    ublas::vector<double> values_i(n);
     ublas::matrix<double, ublas::column_major> Z(n,n);
 
     cout << "\nHSEQR for only eigenvalues." << endl;
     ublas::matrix<double, ublas::column_major> Z_dummy(1,1);
-    lapack::hseqr('E', 'N', 1, n, H, values_r, values_i, Z_dummy);
-    traits::detail::interlace(values_r.begin(), values_r.end(), values_i.begin(), values.begin());
+    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)),
+        Z_dummy);
+    bindings::traits::detail::interlace(bindings::begin_value(values), bindings::size(values));
     cout << "\nH:\n" << H << endl;
     cout << "\nvalues: " << values << endl;
 
     cout << "\nHSEQR for eigenvalues and Schur vectors." << endl;
     Hessenberg(H);
     cout << "H:\n" << H << endl;
-    lapack::hseqr('S', 'I', 1, n, H, values_r, values_i, Z);
-    traits::detail::interlace(values_r.begin(), values_r.end(), values_i.begin(), values.begin());
+    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)),
+        Z);
+    bindings::traits::detail::interlace(bindings::begin_value(values), bindings::size(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-11-21 19:42:07 EST (Sun, 21 Nov 2010)
@@ -15,6 +15,7 @@
 #include <boost/numeric/bindings/ublas/vector.hpp>
 #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/ublas/io.hpp>
 #include <boost/type_traits/is_complex.hpp>
 #include <boost/mpl/if.hpp>
@@ -35,13 +36,13 @@
         external_fp select, MatrixA& a, fortran_int_t& sdim, VectorW& w,
         MatrixVS& vs, Workspace work ) {
     typedef typename bindings::value_type< MatrixA >::type value_type;
-    bindings::detail::array<value_type> wr(bindings::size(w));
-    bindings::detail::array<value_type> wi(bindings::size(w));
-    fortran_int_t info = lapack::gees( jobvs, sort, select, a, sdim, wr, wi, vs, work );
-    traits::detail::interlace(bindings::begin_value(wr),
-                              bindings::begin_value(wr)+bindings::size(w),
-                              bindings::begin_value(wi),
-                              bindings::begin_value(w));
+    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)),
+      vs, work );
+    traits::detail::interlace(bindings::begin_value(w), bindings::size(w));
     return info;
   }
 };

Added: sandbox/numeric_bindings/libs/numeric/bindings/test/Jamfile.v2
==============================================================================
--- (empty file)
+++ sandbox/numeric_bindings/libs/numeric/bindings/test/Jamfile.v2	2010-11-21 19:42:07 EST (Sun, 21 Nov 2010)
@@ -0,0 +1,15 @@
+# Copyright Thomas Klimpel 2008.
+# Use, modification and distribution are subject to 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)
+
+project libs/numeric/bindings/test : requirements
+        <include>$(BOOST_ROOT)
+        <include>$(B_ROOT) ;
+
+import testing ;
+
+alias bindings-tests :
+    [ run interlace.cpp ]
+;
+

Added: sandbox/numeric_bindings/libs/numeric/bindings/test/interlace.cpp
==============================================================================
--- (empty file)
+++ sandbox/numeric_bindings/libs/numeric/bindings/test/interlace.cpp	2010-11-21 19:42:07 EST (Sun, 21 Nov 2010)
@@ -0,0 +1,61 @@
+//
+// 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)
+//
+
+#include <iostream>
+#include <vector>
+#include <boost/numeric/bindings/traits/detail/utils.hpp>
+
+int test_inshuffle(std::size_t n) {
+  std::vector<std::size_t> data(2*n);
+  for (std::size_t i = 0; i < n; ++i) {
+    data[i] = 2*i+1;
+    data[i+n] = 2*i;
+  }
+  boost::numeric::bindings::traits::detail::inshuffle(&data[0],n);
+  for (std::size_t i = 1; i < 2*n; ++i)
+    if (data[i-1]+1 != data[i]) {
+      std::cout << "Test inshuffle for n=" << n << std::endl;
+      for (std::size_t j = 0; j < 2*n; ++j)
+        std::cout << " " << data[j];
+      std::cout << std::endl;
+      std::cout << "logic error" << std::endl;
+      return 255;
+    }
+  return 0;
+}
+
+int test_interlace(std::size_t n) {
+  std::vector<std::complex<double> > data(n);
+  for (std::size_t i = 0; i < n; ++i)
+    data[i] = std::complex<double>(2*i,2*i+1);
+  boost::numeric::bindings::traits::detail::interlace(&data[0],n);
+  for (std::size_t i = 0; i < n; ++i)
+    if (data[i].real() != i || data[i].imag() != i+n) {
+      std::cout << "Test interlace for n=" << n << std::endl;
+      for (std::size_t j = 0; j < n; ++j)
+        std::cout << " " << data[j];
+      std::cout << std::endl;
+      std::cout << "logic error" << std::endl;
+      return 255;
+    }
+  return 0;
+}
+
+int main() {
+  for (std::size_t n = 1; n <= 1000; ++n) {
+    int success = test_inshuffle(n);
+    if (success != 0)
+      return success;
+  }
+  for (std::size_t n = 1; n <= 1000; ++n) {
+    int success = test_interlace(n);
+    if (success != 0)
+      return success;
+  }
+  return 0;
+}