diff --git a/include/boost/numeric/ublas/hermitian.hpp b/include/boost/numeric/ublas/hermitian.hpp index 88dc0d5c..600f9678 100644 --- a/include/boost/numeric/ublas/hermitian.hpp +++ b/include/boost/numeric/ublas/hermitian.hpp @@ -18,6 +18,7 @@ #define BOOST_UBLAS_HERMITIAN_H #include +#include // for resize_preserve #include // Iterators based on ideas of Jeremy Siek @@ -317,7 +318,7 @@ namespace boost { namespace numeric { namespace ublas { void resize (size_type size, bool preserve = true) { if (preserve) { self_type temporary (size, size); - detail::matrix_resize_preserve (*this, temporary); + detail::matrix_resize_preserve (*this, temporary); } else { data ().resize (triangular_type::packed_size (layout_type (), size, size)); diff --git a/include/boost/numeric/ublas/triangular.hpp b/include/boost/numeric/ublas/triangular.hpp index 3979282a..3b56dfae 100644 --- a/include/boost/numeric/ublas/triangular.hpp +++ b/include/boost/numeric/ublas/triangular.hpp @@ -25,6 +25,42 @@ namespace boost { namespace numeric { namespace ublas { + namespace detail { + using namespace boost::numeric::ublas; + + // Matrix resizing algorithm + template + BOOST_UBLAS_INLINE + void matrix_resize_preserve (M& m, M& temporary) { + typedef L layout_type; + typedef T triangular_type; + typedef typename M::size_type size_type; + const size_type msize1 (m.size1 ()); // original size + const size_type msize2 (m.size2 ()); + const size_type size1 (temporary.size1 ()); // new size is specified by temporary + const size_type size2 (temporary.size2 ()); + // Common elements to preserve + const size_type size1_min = (std::min) (size1, msize1); + const size_type size2_min = (std::min) (size2, msize2); + // Order for major and minor sizes + const size_type major_size = layout_type::size_M (size1_min, size2_min); + const size_type minor_size = layout_type::size_m (size1_min, size2_min); + // Indexing copy over major + for (size_type major = 0; major != major_size; ++major) { + for (size_type minor = 0; minor != minor_size; ++minor) { + // find indexes - use invertability of element_ functions + const size_type i1 = layout_type::index_M(major, minor); + const size_type i2 = layout_type::index_m(major, minor); + if ( triangular_type::other(i1,i2) ) { + temporary.data () [triangular_type::element (layout_type (), i1, size1, i2, size2)] = + m.data() [triangular_type::element (layout_type (), i1, msize1, i2, msize2)]; + } + } + } + m.assign_temporary (temporary); + } + } + // Array based triangular matrix class template class triangular_matrix: @@ -104,7 +140,7 @@ namespace boost { namespace numeric { namespace ublas { void resize (size_type size1, size_type size2, bool preserve = true) { if (preserve) { self_type temporary (size1, size2); - detail::matrix_resize_preserve (*this, temporary); + detail::matrix_resize_preserve (*this, temporary); } else { data ().resize (triangular_type::packed_size (layout_type (), size1, size2));