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_VECTOR_ASSIGN_ Chris@16: #define _BOOST_UBLAS_VECTOR_ASSIGN_ Chris@16: Chris@16: #include // scalar_assign 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 vector 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 vector_expression &e1, const vector_expression &e2, S epsilon, S min_norm) { Chris@101: 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 vector_expression &e1, const vector_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: // Make sparse proxies conformant 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 (V &v, const vector_expression &e) { Chris@16: BOOST_UBLAS_CHECK (v.size () == e ().size (), bad_size ()); Chris@16: typedef typename V::size_type size_type; Chris@16: typedef typename V::difference_type difference_type; Chris@16: typedef typename V::value_type value_type; Chris@16: // FIXME unbounded_array with push_back maybe better Chris@16: std::vector index; Chris@16: typename V::iterator it (v.begin ()); Chris@16: typename V::iterator it_end (v.end ()); Chris@16: typename E::const_iterator ite (e ().begin ()); Chris@16: typename E::const_iterator ite_end (e ().end ()); Chris@16: if (it != it_end && ite != ite_end) { Chris@16: size_type it_index = it.index (), ite_index = ite.index (); Chris@16: while (true) { Chris@16: difference_type compare = it_index - ite_index; Chris@16: if (compare == 0) { Chris@16: ++ it, ++ ite; Chris@16: if (it != it_end && ite != ite_end) { Chris@16: it_index = it.index (); Chris@16: ite_index = ite.index (); Chris@16: } else Chris@16: break; Chris@16: } else if (compare < 0) { Chris@16: increment (it, it_end, - compare); Chris@16: if (it != it_end) Chris@16: it_index = it.index (); Chris@16: else Chris@16: break; Chris@16: } else if (compare > 0) { Chris@16: if (*ite != value_type/*zero*/()) Chris@16: index.push_back (ite.index ()); Chris@16: ++ ite; Chris@16: if (ite != ite_end) Chris@16: ite_index = ite.index (); Chris@16: else Chris@16: break; Chris@16: } Chris@16: } Chris@16: } Chris@16: Chris@16: while (ite != ite_end) { Chris@16: if (*ite != value_type/*zero*/()) Chris@16: index.push_back (ite.index ()); Chris@16: ++ ite; Chris@16: } Chris@16: for (size_type k = 0; k < index.size (); ++ k) Chris@16: v (index [k]) = value_type/*zero*/(); Chris@16: } Chris@16: Chris@16: }//namespace detail Chris@16: Chris@16: Chris@16: // Explicitly iterating Chris@16: template