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