Chris@16: // Chris@16: // Copyright (c) 2000-2002 Chris@16: // Joerg Walter, Mathias Koch Chris@16: // Chris@16: // Distributed under the Boost Software License, Version 1.0. (See Chris@16: // accompanying file LICENSE_1_0.txt or copy at Chris@16: // http://www.boost.org/LICENSE_1_0.txt) Chris@16: // Chris@16: // The authors gratefully acknowledge the support of Chris@16: // GeNeSys mbH & Co. KG in producing this work. Chris@16: // Chris@16: Chris@16: #ifndef _BOOST_UBLAS_MATRIX_ASSIGN_ Chris@16: #define _BOOST_UBLAS_MATRIX_ASSIGN_ Chris@16: Chris@16: #include Chris@16: // Required for make_conformant storage Chris@16: #include Chris@16: Chris@16: // Iterators based on ideas of Jeremy Siek Chris@16: Chris@16: namespace boost { namespace numeric { namespace ublas { Chris@16: namespace detail { Chris@16: Chris@16: // Weak equality check - useful to compare equality two arbitary matrix expression results. Chris@16: // Since the actual expressions are unknown, we check for and arbitary error bound Chris@16: // on the relative error. Chris@16: // For a linear expression the infinity norm makes sense as we do not know how the elements will be Chris@16: // combined in the expression. False positive results are inevitable for arbirary expressions! Chris@16: template Chris@16: BOOST_UBLAS_INLINE Chris@16: bool equals (const matrix_expression &e1, const matrix_expression &e2, S epsilon, S min_norm) { Chris@16: return norm_inf (e1 - e2) < epsilon * Chris@16: std::max (std::max (norm_inf (e1), norm_inf (e2)), min_norm); Chris@16: } Chris@16: Chris@16: template Chris@16: BOOST_UBLAS_INLINE Chris@16: bool expression_type_check (const matrix_expression &e1, const matrix_expression &e2) { Chris@16: typedef typename type_traits::promote_type>::real_type real_type; Chris@16: return equals (e1, e2, BOOST_UBLAS_TYPE_CHECK_EPSILON, BOOST_UBLAS_TYPE_CHECK_MIN); Chris@16: } Chris@16: Chris@16: Chris@16: template Chris@16: // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. Chris@16: void make_conformant (M &m, const matrix_expression &e, row_major_tag, R) { Chris@16: BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ()); Chris@16: BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ()); Chris@16: typedef R conformant_restrict_type; Chris@16: typedef typename M::size_type size_type; Chris@16: typedef typename M::difference_type difference_type; Chris@16: typedef typename M::value_type value_type; Chris@16: // FIXME unbounded_array with push_back maybe better Chris@16: std::vector > index; Chris@16: typename M::iterator1 it1 (m.begin1 ()); Chris@16: typename M::iterator1 it1_end (m.end1 ()); Chris@16: typename E::const_iterator1 it1e (e ().begin1 ()); Chris@16: typename E::const_iterator1 it1e_end (e ().end1 ()); Chris@16: while (it1 != it1_end && it1e != it1e_end) { Chris@16: difference_type compare = it1.index1 () - it1e.index1 (); Chris@16: if (compare == 0) { Chris@16: #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION Chris@16: typename M::iterator2 it2 (it1.begin ()); Chris@16: typename M::iterator2 it2_end (it1.end ()); Chris@16: typename E::const_iterator2 it2e (it1e.begin ()); Chris@16: typename E::const_iterator2 it2e_end (it1e.end ()); Chris@16: #else Chris@16: typename M::iterator2 it2 (begin (it1, iterator1_tag ())); Chris@16: typename M::iterator2 it2_end (end (it1, iterator1_tag ())); Chris@16: typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ())); Chris@16: typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ())); Chris@16: #endif Chris@16: if (it2 != it2_end && it2e != it2e_end) { Chris@16: size_type it2_index = it2.index2 (), it2e_index = it2e.index2 (); Chris@16: while (true) { Chris@101: difference_type compare2 = it2_index - it2e_index; Chris@101: if (compare2 == 0) { Chris@16: ++ it2, ++ it2e; Chris@16: if (it2 != it2_end && it2e != it2e_end) { Chris@16: it2_index = it2.index2 (); Chris@16: it2e_index = it2e.index2 (); Chris@16: } else Chris@16: break; Chris@101: } else if (compare2 < 0) { Chris@101: increment (it2, it2_end, - compare2); Chris@16: if (it2 != it2_end) Chris@16: it2_index = it2.index2 (); Chris@16: else Chris@16: break; Chris@101: } else if (compare2 > 0) { Chris@16: if (conformant_restrict_type::other (it2e.index1 (), it2e.index2 ())) Chris@16: if (static_cast(*it2e) != value_type/*zero*/()) Chris@16: index.push_back (std::pair (it2e.index1 (), it2e.index2 ())); Chris@16: ++ it2e; Chris@16: if (it2e != it2e_end) Chris@16: it2e_index = it2e.index2 (); Chris@16: else Chris@16: break; Chris@16: } Chris@16: } Chris@16: } Chris@16: while (it2e != it2e_end) { Chris@16: if (conformant_restrict_type::other (it2e.index1 (), it2e.index2 ())) Chris@16: if (static_cast(*it2e) != value_type/*zero*/()) Chris@16: index.push_back (std::pair (it2e.index1 (), it2e.index2 ())); Chris@16: ++ it2e; Chris@16: } Chris@16: ++ it1, ++ it1e; Chris@16: } else if (compare < 0) { Chris@16: increment (it1, it1_end, - compare); Chris@16: } else if (compare > 0) { Chris@16: #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION Chris@16: typename E::const_iterator2 it2e (it1e.begin ()); Chris@16: typename E::const_iterator2 it2e_end (it1e.end ()); Chris@16: #else Chris@16: typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ())); Chris@16: typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ())); Chris@16: #endif Chris@16: while (it2e != it2e_end) { Chris@16: if (conformant_restrict_type::other (it2e.index1 (), it2e.index2 ())) Chris@16: if (static_cast(*it2e) != value_type/*zero*/()) Chris@16: index.push_back (std::pair (it2e.index1 (), it2e.index2 ())); Chris@16: ++ it2e; Chris@16: } Chris@16: ++ it1e; Chris@16: } Chris@16: } Chris@16: while (it1e != it1e_end) { Chris@16: #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION Chris@16: typename E::const_iterator2 it2e (it1e.begin ()); Chris@16: typename E::const_iterator2 it2e_end (it1e.end ()); Chris@16: #else Chris@16: typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ())); Chris@16: typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ())); Chris@16: #endif Chris@16: while (it2e != it2e_end) { Chris@16: if (conformant_restrict_type::other (it2e.index1 (), it2e.index2 ())) Chris@16: if (static_cast(*it2e) != value_type/*zero*/()) Chris@16: index.push_back (std::pair (it2e.index1 (), it2e.index2 ())); Chris@16: ++ it2e; Chris@16: } Chris@16: ++ it1e; Chris@16: } Chris@16: // ISSUE proxies require insert_element Chris@16: for (size_type k = 0; k < index.size (); ++ k) Chris@16: m (index [k].first, index [k].second) = value_type/*zero*/(); Chris@16: } Chris@16: template Chris@16: // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. Chris@16: void make_conformant (M &m, const matrix_expression &e, column_major_tag, R) { Chris@16: BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ()); Chris@16: BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ()); Chris@16: typedef R conformant_restrict_type; Chris@16: typedef typename M::size_type size_type; Chris@16: typedef typename M::difference_type difference_type; Chris@16: typedef typename M::value_type value_type; Chris@16: std::vector > index; Chris@16: typename M::iterator2 it2 (m.begin2 ()); Chris@16: typename M::iterator2 it2_end (m.end2 ()); Chris@16: typename E::const_iterator2 it2e (e ().begin2 ()); Chris@16: typename E::const_iterator2 it2e_end (e ().end2 ()); Chris@16: while (it2 != it2_end && it2e != it2e_end) { Chris@16: difference_type compare = it2.index2 () - it2e.index2 (); Chris@16: if (compare == 0) { Chris@16: #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION Chris@16: typename M::iterator1 it1 (it2.begin ()); Chris@16: typename M::iterator1 it1_end (it2.end ()); Chris@16: typename E::const_iterator1 it1e (it2e.begin ()); Chris@16: typename E::const_iterator1 it1e_end (it2e.end ()); Chris@16: #else Chris@16: typename M::iterator1 it1 (begin (it2, iterator2_tag ())); Chris@16: typename M::iterator1 it1_end (end (it2, iterator2_tag ())); Chris@16: typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ())); Chris@16: typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ())); Chris@16: #endif Chris@16: if (it1 != it1_end && it1e != it1e_end) { Chris@16: size_type it1_index = it1.index1 (), it1e_index = it1e.index1 (); Chris@16: while (true) { Chris@101: difference_type compare2 = it1_index - it1e_index; Chris@101: if (compare2 == 0) { Chris@16: ++ it1, ++ it1e; Chris@16: if (it1 != it1_end && it1e != it1e_end) { Chris@16: it1_index = it1.index1 (); Chris@16: it1e_index = it1e.index1 (); Chris@16: } else Chris@16: break; Chris@101: } else if (compare2 < 0) { Chris@101: increment (it1, it1_end, - compare2); Chris@16: if (it1 != it1_end) Chris@16: it1_index = it1.index1 (); Chris@16: else Chris@16: break; Chris@101: } else if (compare2 > 0) { Chris@16: if (conformant_restrict_type::other (it1e.index1 (), it1e.index2 ())) Chris@16: if (static_cast(*it1e) != value_type/*zero*/()) Chris@16: index.push_back (std::pair (it1e.index1 (), it1e.index2 ())); Chris@16: ++ it1e; Chris@16: if (it1e != it1e_end) Chris@16: it1e_index = it1e.index1 (); Chris@16: else Chris@16: break; Chris@16: } Chris@16: } Chris@16: } Chris@16: while (it1e != it1e_end) { Chris@16: if (conformant_restrict_type::other (it1e.index1 (), it1e.index2 ())) Chris@16: if (static_cast(*it1e) != value_type/*zero*/()) Chris@16: index.push_back (std::pair (it1e.index1 (), it1e.index2 ())); Chris@16: ++ it1e; Chris@16: } Chris@16: ++ it2, ++ it2e; Chris@16: } else if (compare < 0) { Chris@16: increment (it2, it2_end, - compare); Chris@16: } else if (compare > 0) { Chris@16: #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION Chris@16: typename E::const_iterator1 it1e (it2e.begin ()); Chris@16: typename E::const_iterator1 it1e_end (it2e.end ()); Chris@16: #else Chris@16: typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ())); Chris@16: typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ())); Chris@16: #endif Chris@16: while (it1e != it1e_end) { Chris@16: if (conformant_restrict_type::other (it1e.index1 (), it1e.index2 ())) Chris@16: if (static_cast(*it1e) != value_type/*zero*/()) Chris@16: index.push_back (std::pair (it1e.index1 (), it1e.index2 ())); Chris@16: ++ it1e; Chris@16: } Chris@16: ++ it2e; Chris@16: } Chris@16: } Chris@16: while (it2e != it2e_end) { Chris@16: #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION Chris@16: typename E::const_iterator1 it1e (it2e.begin ()); Chris@16: typename E::const_iterator1 it1e_end (it2e.end ()); Chris@16: #else Chris@16: typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ())); Chris@16: typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ())); Chris@16: #endif Chris@16: while (it1e != it1e_end) { Chris@16: if (conformant_restrict_type::other (it1e.index1 (), it1e.index2 ())) Chris@16: if (static_cast(*it1e) != value_type/*zero*/()) Chris@16: index.push_back (std::pair (it1e.index1 (), it1e.index2 ())); Chris@16: ++ it1e; Chris@16: } Chris@16: ++ it2e; Chris@16: } Chris@16: // ISSUE proxies require insert_element Chris@16: for (size_type k = 0; k < index.size (); ++ k) Chris@16: m (index [k].first, index [k].second) = value_type/*zero*/(); Chris@16: } Chris@16: Chris@16: }//namespace detail Chris@16: Chris@16: Chris@16: // Explicitly iterating row major Chris@16: template