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;
+}