Squashed 'third_party/boostorg/ublas/' content from commit e8607b3
Change-Id: Ia06afd642157a24e17fa9ddea28fb8601810b78e
git-subtree-dir: third_party/boostorg/ublas
git-subtree-split: e8607b3eea238e590eca93bfe498c21f470155c1
diff --git a/include/boost/numeric/ublas/detail/concepts.hpp b/include/boost/numeric/ublas/detail/concepts.hpp
new file mode 100644
index 0000000..674c610
--- /dev/null
+++ b/include/boost/numeric/ublas/detail/concepts.hpp
@@ -0,0 +1,1575 @@
+//
+// Copyright (c) 2000-2002
+// Joerg Walter, Mathias Koch
+//
+// 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)
+//
+// The authors gratefully acknowledge the support of
+// GeNeSys mbH & Co. KG in producing this work.
+//
+
+#ifndef _BOOST_UBLAS_CONCEPTS_
+#define _BOOST_UBLAS_CONCEPTS_
+
+#include <boost/concept_check.hpp>
+
+// Concept checks based on ideas of Jeremy Siek
+
+namespace boost { namespace numeric { namespace ublas {
+
+
+ template<class I>
+ struct Indexed1DIteratorConcept {
+ typedef I iterator_type;
+
+ void constraints () {
+ iterator_type it = iterator_type ();
+ // Index
+ it.index ();
+ }
+ };
+
+ template<class I>
+ struct IndexedBidirectional1DIteratorConcept {
+ typedef I iterator_type;
+
+ void constraints () {
+ function_requires< BidirectionalIteratorConcept<iterator_type> >();
+ function_requires< Indexed1DIteratorConcept<iterator_type> >();
+ }
+ };
+
+ template<class I>
+ struct Mutable_IndexedBidirectional1DIteratorConcept {
+ typedef I iterator_type;
+
+ void constraints () {
+ function_requires< Mutable_BidirectionalIteratorConcept<iterator_type> >();
+ function_requires< Indexed1DIteratorConcept<iterator_type> >();
+ }
+ };
+
+ template<class I>
+ struct IndexedRandomAccess1DIteratorConcept {
+ typedef I iterator_type;
+
+ void constraints () {
+ function_requires< RandomAccessIteratorConcept<iterator_type> >();
+ function_requires< Indexed1DIteratorConcept<iterator_type> >();
+ }
+ };
+
+ template<class I>
+ struct Mutable_IndexedRandomAccess1DIteratorConcept {
+ typedef I iterator_type;
+
+ void constraints () {
+ function_requires< Mutable_RandomAccessIteratorConcept<iterator_type> >();
+ function_requires< Indexed1DIteratorConcept<iterator_type> >();
+ }
+ };
+
+ template<class I>
+ struct Indexed2DIteratorConcept {
+ typedef I iterator_type;
+ typedef typename I::dual_iterator_type dual_iterator_type;
+ typedef typename I::dual_reverse_iterator_type dual_reverse_iterator_type;
+
+ void constraints () {
+ iterator_type it = iterator_type ();
+ // Indices
+ it.index1 ();
+ it.index2 ();
+ // Iterator begin/end
+ dual_iterator_type it_begin (it.begin ());
+ dual_iterator_type it_end (it.end ());
+ // Reverse iterator begin/end
+ dual_reverse_iterator_type it_rbegin (it.rbegin ());
+ dual_reverse_iterator_type it_rend (it.rend ());
+ ignore_unused_variable_warning (it_begin);
+ ignore_unused_variable_warning (it_end);
+ ignore_unused_variable_warning (it_rbegin);
+ ignore_unused_variable_warning (it_rend);
+ }
+ };
+
+ template<class I1, class I2>
+ struct IndexedBidirectional2DIteratorConcept {
+ typedef I1 subiterator1_type;
+ typedef I2 subiterator2_type;
+
+ void constraints () {
+ function_requires< BidirectionalIteratorConcept<subiterator1_type> >();
+ function_requires< BidirectionalIteratorConcept<subiterator2_type> >();
+ function_requires< Indexed2DIteratorConcept<subiterator1_type> >();
+ function_requires< Indexed2DIteratorConcept<subiterator2_type> >();
+ }
+ };
+
+ template<class I1, class I2>
+ struct Mutable_IndexedBidirectional2DIteratorConcept {
+ typedef I1 subiterator1_type;
+ typedef I2 subiterator2_type;
+
+ void constraints () {
+ function_requires< Mutable_BidirectionalIteratorConcept<subiterator1_type> >();
+ function_requires< Mutable_BidirectionalIteratorConcept<subiterator2_type> >();
+ function_requires< Indexed2DIteratorConcept<subiterator1_type> >();
+ function_requires< Indexed2DIteratorConcept<subiterator2_type> >();
+ }
+ };
+
+ template<class I1, class I2>
+ struct IndexedRandomAccess2DIteratorConcept {
+ typedef I1 subiterator1_type;
+ typedef I2 subiterator2_type;
+
+ void constraints () {
+ function_requires< RandomAccessIteratorConcept<subiterator1_type> >();
+ function_requires< RandomAccessIteratorConcept<subiterator2_type> >();
+ function_requires< Indexed2DIteratorConcept<subiterator1_type> >();
+ function_requires< Indexed2DIteratorConcept<subiterator2_type> >();
+ }
+ };
+
+ template<class I1, class I2>
+ struct Mutable_IndexedRandomAccess2DIteratorConcept {
+ typedef I1 subiterator1_type;
+ typedef I2 subiterator2_type;
+
+ void constraints () {
+ function_requires< Mutable_RandomAccessIteratorConcept<subiterator1_type> >();
+ function_requires< Mutable_RandomAccessIteratorConcept<subiterator2_type> >();
+ function_requires< Indexed2DIteratorConcept<subiterator1_type> >();
+ function_requires< Indexed2DIteratorConcept<subiterator2_type> >();
+ }
+ };
+
+ template<class C>
+ struct StorageArrayConcept {
+ typedef C container_type;
+ typedef typename C::size_type size_type;
+ typedef typename C::value_type value_type;
+
+ void constraints () {
+ function_requires< RandomAccessContainerConcept<container_type> >();
+ size_type n (0);
+ // Sizing constructor
+ container_type c = container_type (n);
+ // Initialised sizing constructor
+ container_type (n, value_type (5));
+ ignore_unused_variable_warning (c);
+ }
+ };
+
+ template<class C>
+ struct Mutable_StorageArrayConcept {
+ typedef C container_type;
+ typedef typename C::size_type size_type;
+ typedef typename C::value_type value_type;
+ typedef typename C::iterator iterator_type;
+
+ void constraints () {
+ function_requires< Mutable_RandomAccessContainerConcept<container_type> > ();
+ size_type n (0);
+ // Sizing constructor
+ container_type c = container_type (n);
+ // Initialised sizing constructor
+ c = container_type (n, value_type (3));
+ // Resize
+ c.resize (n, value_type (5));
+ // Resize - none preserving
+ c.resize (n);
+ }
+ };
+
+ template<class C>
+ struct StorageSparseConcept {
+ typedef C container_type;
+ typedef typename C::size_type size_type;
+
+ void constraints () {
+ function_requires< ReversibleContainerConcept<container_type> > ();
+ }
+ };
+
+ template<class C>
+ struct Mutable_StorageSparseConcept {
+ typedef C container_type;
+ typedef typename C::size_type size_type;
+ typedef typename C::value_type value_type;
+ typedef typename C::iterator iterator_type;
+
+ void constraints () {
+ // NOTE - Not Mutable_ReversibleContainerConcept
+ function_requires< ReversibleContainerConcept<container_type> >();
+ container_type c = container_type ();
+ value_type t = value_type ();
+ iterator_type it = iterator_type (), it1 = iterator_type (), it2 = iterator_type ();
+ // Insert
+ c.insert (it, t);
+ // Erase
+ c.erase (it);
+ // Range erase
+ c.erase (it1, it2);
+ // Clear
+ c.clear ();
+ }
+ };
+
+ template<class G>
+ struct IndexSetConcept {
+ typedef G generator_type;
+ typedef typename G::size_type size_type;
+ typedef typename G::value_type value_type;
+
+ void constraints () {
+ function_requires< AssignableConcept<generator_type> >();
+ function_requires< ReversibleContainerConcept<generator_type> >();
+ generator_type g = generator_type ();
+ size_type n (0);
+ value_type t;
+ // Element access
+ t = g (n);
+ ignore_unused_variable_warning (t);
+ }
+ };
+
+ /** \brief Scalar expression concept.
+ *
+ * requirements
+ * \li \c SE::value_type is the type of the scalar expression
+ * \li \c SE must be convertable to \c SE::value_type
+ * \li the constant \c SE::complexity must exist
+ *
+ * \param SE the type of the scalar expression
+ */
+ template<class SE>
+ struct ScalarExpressionConcept {
+ typedef SE scalar_expression_type;
+ typedef typename SE::value_type value_type;
+
+ static const unsigned complexity = SE::complexity;
+
+ void constraints () {
+ scalar_expression_type *sp;
+ scalar_expression_type s = *sp;
+ value_type t;
+ // Conversion
+ t = s;
+ ignore_unused_variable_warning (t);
+ }
+ };
+
+ /** \brief Vector expression concept.
+ *
+ * requirements
+ * \li \c VE::value_type is the type of the elements
+ * \li \c VE::const_reference The return type when accessing an element of a constant vector
+ * expression. Must be convertable to a \c value_type.
+ * \li \c VE::size_type is the (unsigned) type of the indices
+ * \li \c VE::difference_type is the (signed) type of distances between indices
+ * \li \c VE::category
+ *
+ * \li the constant \c SE::complexity must exist
+ *
+ * \param SE the type of the scalar expression
+ */
+ template<class VE>
+ struct VectorExpressionConcept {
+ typedef VE vector_expression_type;
+ typedef typename VE::type_category type_category;
+ typedef typename VE::size_type size_type;
+ typedef typename VE::difference_type difference_type;
+ typedef typename VE::value_type value_type;
+ typedef typename VE::const_reference const_reference;
+ typedef typename VE::const_iterator const_iterator_type;
+ typedef typename VE::const_reverse_iterator const_reverse_iterator_type;
+
+ void constraints () {
+ vector_expression_type *vp;
+ const vector_expression_type *cvp;
+ vector_expression_type v = *vp;
+ const vector_expression_type cv = *cvp;
+ size_type n (0), i (0);
+ value_type t;
+ // Find (internal?)
+ const_iterator_type cit (v.find (i));
+ // Beginning of range
+ const_iterator_type cit_begin (v.begin ());
+ // End of range
+ const_iterator_type cit_end (v.end ());
+ // Size
+ n = v.size ();
+ // Beginning of reverse range
+ const_reverse_iterator_type crit_begin (cv.rbegin ());
+ // End of reverse range
+ const_reverse_iterator_type crit_end (cv.rend ());
+ // Element access
+ t = v (i);
+ ignore_unused_variable_warning (n);
+ ignore_unused_variable_warning (cit);
+ ignore_unused_variable_warning (cit_begin);
+ ignore_unused_variable_warning (cit_end);
+ ignore_unused_variable_warning (crit_begin);
+ ignore_unused_variable_warning (crit_end);
+ ignore_unused_variable_warning (t);
+ }
+ };
+
+ template<class VE>
+ struct Mutable_VectorExpressionConcept {
+ typedef VE vector_expression_type;
+ typedef typename VE::size_type size_type;
+ typedef typename VE::value_type value_type;
+ typedef typename VE::iterator iterator_type;
+ typedef typename VE::reverse_iterator reverse_iterator_type;
+
+ void constraints () {
+ function_requires< AssignableConcept<vector_expression_type> >();
+ function_requires< VectorExpressionConcept<vector_expression_type> >();
+ vector_expression_type *vp;
+ vector_expression_type v = *vp, v1 = *vp, v2 = *vp;
+ size_type i (0);
+ value_type t = value_type ();
+ // Find (internal?)
+ iterator_type it (v.find (i));
+ // Beginning of range
+ iterator_type it_begin (v.begin ());
+ // End of range
+ iterator_type it_end (v.end ());
+ // Swap
+ v1.swap (v2);
+ // Beginning of reverse range
+ reverse_iterator_type rit_begin (v.rbegin ());
+ // End of reverse range
+ reverse_iterator_type rit_end (v.rend ());
+ // Assignments
+ v2 = v1;
+ v2.assign (v1);
+ v2 += v1;
+ v2.plus_assign (v1);
+ v2 -= v1;
+ v2.minus_assign (v1);
+ v *= t;
+ ignore_unused_variable_warning (it);
+ ignore_unused_variable_warning (it_begin);
+ ignore_unused_variable_warning (it_end);
+ ignore_unused_variable_warning (rit_begin);
+ ignore_unused_variable_warning (rit_end);
+ }
+ };
+
+ template<class ME>
+ struct MatrixExpressionConcept {
+ typedef ME matrix_expression_type;
+ typedef typename ME::type_category type_category;
+ typedef typename ME::size_type size_type;
+ typedef typename ME::value_type value_type;
+ typedef typename ME::const_iterator1 const_subiterator1_type;
+ typedef typename ME::const_iterator2 const_subiterator2_type;
+ typedef typename ME::const_reverse_iterator1 const_reverse_subiterator1_type;
+ typedef typename ME::const_reverse_iterator2 const_reverse_subiterator2_type;
+
+ void constraints () {
+ matrix_expression_type *mp;
+ const matrix_expression_type *cmp;
+ matrix_expression_type m = *mp;
+ const matrix_expression_type cm = *cmp;
+ size_type n (0), i (0), j (0);
+ value_type t;
+ // Find (internal?)
+ const_subiterator1_type cit1 (m.find1 (0, i, j));
+ const_subiterator2_type cit2 (m.find2 (0, i, j));
+ // Beginning of range
+ const_subiterator1_type cit1_begin (m.begin1 ());
+ const_subiterator2_type cit2_begin (m.begin2 ());
+ // End of range
+ const_subiterator1_type cit1_end (m.end1 ());
+ const_subiterator2_type cit2_end (m.end2 ());
+ // Size
+ n = m.size1 ();
+ n = m.size2 ();
+ // Beginning of reverse range
+ const_reverse_subiterator1_type crit1_begin (cm.rbegin1 ());
+ const_reverse_subiterator2_type crit2_begin (cm.rbegin2 ());
+ // End of reverse range
+ const_reverse_subiterator1_type crit1_end (cm.rend1 ());
+ const_reverse_subiterator2_type crit2_end (cm.rend2 ());
+ // Element access
+ t = m (i, j);
+ ignore_unused_variable_warning (n);
+ ignore_unused_variable_warning (cit1);
+ ignore_unused_variable_warning (cit2);
+ ignore_unused_variable_warning (cit1_begin);
+ ignore_unused_variable_warning (cit2_begin);
+ ignore_unused_variable_warning (cit1_end);
+ ignore_unused_variable_warning (cit2_end);
+ ignore_unused_variable_warning (crit1_begin);
+ ignore_unused_variable_warning (crit2_begin);
+ ignore_unused_variable_warning (crit1_end);
+ ignore_unused_variable_warning (crit2_end);
+ ignore_unused_variable_warning (t);
+ }
+ };
+
+ template<class ME>
+ struct Mutable_MatrixExpressionConcept {
+ typedef ME matrix_expression_type;
+ typedef typename ME::size_type size_type;
+ typedef typename ME::value_type value_type;
+ typedef typename ME::iterator1 subiterator1_type;
+ typedef typename ME::iterator2 subiterator2_type;
+ typedef typename ME::reverse_iterator1 reverse_subiterator1_type;
+ typedef typename ME::reverse_iterator2 reverse_subiterator2_type;
+
+ void constraints () {
+ function_requires< AssignableConcept<matrix_expression_type> >();
+ function_requires< MatrixExpressionConcept<matrix_expression_type> >();
+ matrix_expression_type *mp;
+ matrix_expression_type m = *mp, m1 = *mp, m2 = *mp;
+ size_type i (0), j (0);
+ value_type t = value_type ();
+ // Find (internal?)
+ subiterator1_type it1 (m.find1 (0, i, j));
+ subiterator2_type it2 (m.find2 (0, i, j));
+ // Beginning of range
+ subiterator1_type it1_begin (m.begin1 ());
+ subiterator2_type it2_begin (m.begin2 ());
+ // End of range
+ subiterator1_type it1_end (m.end1 ());
+ subiterator2_type it2_end (m.end2 ());
+ // Swap
+ m1.swap (m2);
+ // Beginning of reverse range
+ reverse_subiterator1_type rit1_begin (m.rbegin1 ());
+ reverse_subiterator2_type rit2_begin (m.rbegin2 ());
+ // End of reverse range
+ reverse_subiterator1_type rit1_end (m.rend1 ());
+ reverse_subiterator2_type rit2_end (m.rend2 ());
+ // Assignments
+ m2 = m1;
+ m2.assign (m1);
+ m2 += m1;
+ m2.plus_assign (m1);
+ m2 -= m1;
+ m2.minus_assign (m1);
+ m *= t;
+ ignore_unused_variable_warning (it1);
+ ignore_unused_variable_warning (it2);
+ ignore_unused_variable_warning (it1_begin);
+ ignore_unused_variable_warning (it2_begin);
+ ignore_unused_variable_warning (it1_end);
+ ignore_unused_variable_warning (it2_end);
+ ignore_unused_variable_warning (rit1_begin);
+ ignore_unused_variable_warning (rit2_begin);
+ ignore_unused_variable_warning (rit1_end);
+ ignore_unused_variable_warning (rit2_end);
+ }
+ };
+
+ template<class V>
+ struct VectorConcept {
+ typedef V vector_type;
+ typedef typename V::size_type size_type;
+ typedef typename V::value_type value_type;
+ typedef const value_type *const_pointer;
+
+ void constraints () {
+ function_requires< VectorExpressionConcept<vector_type> >();
+ size_type n (0);
+ size_type i (0);
+ // Sizing constructor
+ vector_type v (n);
+ // Element support
+ const_pointer p = v.find_element (i);
+
+ ignore_unused_variable_warning (p);
+ }
+ };
+
+ template<class V>
+ struct Mutable_VectorConcept {
+ typedef V vector_type;
+ typedef typename V::size_type size_type;
+ typedef typename V::value_type value_type;
+ typedef value_type *pointer;
+
+ void constraints () {
+ function_requires< VectorConcept<vector_type> >();
+ function_requires< DefaultConstructible<vector_type> >();
+ function_requires< Mutable_VectorExpressionConcept<vector_type> >();
+ size_type n (0);
+ value_type t = value_type ();
+ size_type i (0);
+ vector_type v;
+ // Element support
+ pointer p = v.find_element (i);
+ // Element assignment
+ value_type r = v.insert_element (i, t);
+ v.insert_element (i, t) = r;
+ // Zeroing
+ v.clear ();
+ // Resize
+ v.resize (n);
+
+ ignore_unused_variable_warning (p);
+ ignore_unused_variable_warning (r);
+ }
+ };
+
+ template<class V>
+ struct SparseVectorConcept {
+ typedef V vector_type;
+ typedef typename V::size_type size_type;
+
+ void constraints () {
+ function_requires< VectorConcept<vector_type> >();
+ }
+ };
+
+ template<class V>
+ struct Mutable_SparseVectorConcept {
+ typedef V vector_type;
+ typedef typename V::size_type size_type;
+ typedef typename V::value_type value_type;
+
+ void constraints () {
+ function_requires< SparseVectorConcept<vector_type> >();
+ function_requires< Mutable_VectorConcept<vector_type> >();
+ size_type i (0);
+ vector_type v;
+ // Element erasure
+ v.erase_element (i);
+ }
+ };
+
+ template<class M>
+ struct MatrixConcept {
+ typedef M matrix_type;
+ typedef typename M::size_type size_type;
+ typedef typename M::value_type value_type;
+ typedef const value_type *const_pointer;
+
+ void constraints () {
+ function_requires< MatrixExpressionConcept<matrix_type> >();
+ size_type n (0);
+ size_type i (0), j (0);
+ // Sizing constructor
+ matrix_type m (n, n);
+ // Element support
+#ifndef SKIP_BAD
+ const_pointer p = m.find_element (i, j);
+#else
+ const_pointer p;
+ ignore_unused_variable_warning (i);
+ ignore_unused_variable_warning (j);
+#endif
+ ignore_unused_variable_warning (p);
+ }
+ };
+
+ template<class M>
+ struct Mutable_MatrixConcept {
+ typedef M matrix_type;
+ typedef typename M::size_type size_type;
+ typedef typename M::value_type value_type;
+ typedef value_type *pointer;
+
+ void constraints () {
+ function_requires< MatrixConcept<matrix_type> >();
+ function_requires< DefaultConstructible<matrix_type> >();
+ function_requires< Mutable_MatrixExpressionConcept<matrix_type> >();
+ size_type n (0);
+ value_type t = value_type ();
+ size_type i (0), j (0);
+ matrix_type m;
+ // Element support
+#ifndef SKIP_BAD
+ pointer p = m.find_element (i, j);
+ ignore_unused_variable_warning (i);
+ ignore_unused_variable_warning (j);
+#else
+ pointer p;
+#endif
+ // Element assigment
+ value_type r = m.insert_element (i, j, t);
+ m.insert_element (i, j, t) = r;
+ // Zeroing
+ m.clear ();
+ // Resize
+ m.resize (n, n);
+ m.resize (n, n, false);
+
+ ignore_unused_variable_warning (p);
+ ignore_unused_variable_warning (r);
+ }
+ };
+
+ template<class M>
+ struct SparseMatrixConcept {
+ typedef M matrix_type;
+ typedef typename M::size_type size_type;
+
+ void constraints () {
+ function_requires< MatrixConcept<matrix_type> >();
+ }
+ };
+
+ template<class M>
+ struct Mutable_SparseMatrixConcept {
+ typedef M matrix_type;
+ typedef typename M::size_type size_type;
+ typedef typename M::value_type value_type;
+
+ void constraints () {
+ function_requires< SparseMatrixConcept<matrix_type> >();
+ function_requires< Mutable_MatrixConcept<matrix_type> >();
+ size_type i (0), j (0);
+ matrix_type m;
+ // Elemnent erasure
+ m.erase_element (i, j);
+ }
+ };
+
+ /** introduce anonymous namespace to make following functions
+ * local to the current compilation unit.
+ */
+ namespace {
+
+ // Replaced the ZeroElement and OneElement functions with the templated versions
+ // because the former where giving warnings with clang
+ template<class T>
+ T
+ ZeroElement (T) {
+ return T(0.0);
+ }
+
+ template<class T>
+ vector<T>
+ ZeroElement (vector<T>) {
+ return zero_vector<T> ();
+ }
+
+ template<class T>
+ matrix<T>
+ ZeroElement (matrix<T>) {
+ return zero_matrix<T> ();
+ }
+
+ template<class T>
+ T
+ OneElement (T) {
+ return T(0.0);
+ }
+
+ template<class T>
+ vector<T>
+ OneElement (vector<T>) {
+ return zero_vector<T> ();
+ }
+
+ template<class T>
+ matrix<T>
+ OneElement (matrix<T>) {
+ return identity_matrix<T> ();
+ }
+
+// template<>
+// float
+// ZeroElement (float) {
+// return 0.f;
+// }
+// template<>
+// double
+// ZeroElement (double) {
+// return 0.;
+// }
+// template<>
+// vector<float>
+// ZeroElement (vector<float>) {
+// return zero_vector<float> ();
+// }
+// template<>
+// vector<double>
+// ZeroElement (vector<double>) {
+// return zero_vector<double> ();
+// }
+// template<>
+// matrix<float>
+// ZeroElement (matrix<float>) {
+// return zero_matrix<float> ();
+// }
+// template<>
+// matrix<double>
+// ZeroElement (matrix<double>) {
+// return zero_matrix<double> ();
+// }
+// template<>
+// std::complex<float>
+// ZeroElement (std::complex<float>) {
+// return std::complex<float> (0.f);
+// }
+// template<>
+// std::complex<double>
+// ZeroElement (std::complex<double>) {
+// return std::complex<double> (0.);
+// }
+// template<>
+// vector<std::complex<float> >
+// ZeroElement (vector<std::complex<float> >) {
+// return zero_vector<std::complex<float> > ();
+// }
+// template<>
+// vector<std::complex<double> >
+// ZeroElement (vector<std::complex<double> >) {
+// return zero_vector<std::complex<double> > ();
+// }
+// template<>
+// matrix<std::complex<float> >
+// ZeroElement (matrix<std::complex<float> >) {
+// return zero_matrix<std::complex<float> > ();
+// }
+// template<>
+// matrix<std::complex<double> >
+// ZeroElement (matrix<std::complex<double> >) {
+// return zero_matrix<std::complex<double> > ();
+// }
+
+// template<class T>
+// T
+// OneElement (T);
+// template<>
+// float
+// OneElement (float) {
+// return 1.f;
+// }
+// template<>
+// double
+// OneElement (double) {
+// return 1.;
+// }
+// template<>
+// matrix<float>
+// OneElement (matrix<float>) {
+// return identity_matrix<float> ();
+// }
+// template<>
+// matrix<double>
+// OneElement (matrix<double>) {
+// return identity_matrix<double> ();
+// }
+// template<>
+// std::complex<float>
+// OneElement (std::complex<float>) {
+// return std::complex<float> (1.f);
+// }
+// template<>
+// std::complex<double>
+// OneElement (std::complex<double>) {
+// return std::complex<double> (1.);
+// }
+// template<>
+// matrix<std::complex<float> >
+// OneElement (matrix<std::complex<float> >) {
+// return identity_matrix<std::complex<float> > ();
+// }
+// template<>
+// matrix<std::complex<double> >
+// OneElement (matrix<std::complex<double> >) {
+// return identity_matrix<std::complex<double> > ();
+// }
+
+ template<class E1, class E2>
+ bool
+ operator == (const vector_expression<E1> &e1, const vector_expression<E2> &e2) {
+ typedef typename promote_traits<typename E1::value_type,
+ typename E2::value_type>::promote_type value_type;
+ typedef typename type_traits<value_type>::real_type real_type;
+ return norm_inf (e1 - e2) == real_type/*zero*/();
+ }
+ template<class E1, class E2>
+ bool
+ operator == (const matrix_expression<E1> &e1, const matrix_expression<E2> &e2) {
+ typedef typename promote_traits<typename E1::value_type,
+ typename E2::value_type>::promote_type value_type;
+ typedef typename type_traits<value_type>::real_type real_type;
+ return norm_inf (e1 - e2) == real_type/*zero*/();
+ }
+
+ template<class T>
+ struct AdditiveAbelianGroupConcept {
+ typedef T value_type;
+
+ void constraints () {
+ bool r;
+ value_type a = value_type (), b = value_type (), c = value_type ();
+ r = (a + b) + c == a + (b + c);
+ r = ZeroElement (value_type ()) + a == a;
+ r = a + ZeroElement (value_type ()) == a;
+ r = a + (- a) == ZeroElement (value_type ());
+ r = (- a) + a == ZeroElement (value_type ());
+ r = a + b == b + a;
+ ignore_unused_variable_warning (r);
+ }
+ };
+
+ template<class T>
+ struct MultiplicativeAbelianGroupConcept {
+ typedef T value_type;
+
+ void constraints () {
+ bool r;
+ value_type a = value_type (), b = value_type (), c = value_type ();
+ r = (a * b) * c == a * (b * c);
+ r = OneElement (value_type ()) * a == a;
+ r = a * OneElement (value_type ()) == a;
+ r = a * (OneElement (value_type ()) / a) == a;
+ r = (OneElement (value_type ()) / a) * a == a;
+ r = a * b == b * a;
+ ignore_unused_variable_warning (r);
+ }
+ };
+
+ template<class T>
+ struct RingWithIdentityConcept {
+ typedef T value_type;
+
+ void constraints () {
+ function_requires< AdditiveAbelianGroupConcept<value_type> >();
+ bool r;
+ value_type a = value_type (), b = value_type (), c = value_type ();
+ r = (a * b) * c == a * (b * c);
+ r = (a + b) * c == a * c + b * c;
+ r = OneElement (value_type ()) * a == a;
+ r = a * OneElement (value_type ()) == a;
+ ignore_unused_variable_warning (r);
+ }
+ };
+
+ template<class T>
+ struct Prod_RingWithIdentityConcept {
+ typedef T value_type;
+
+ void constraints () {
+ function_requires< AdditiveAbelianGroupConcept<value_type> >();
+ bool r;
+ value_type a = value_type (), b = value_type (), c = value_type ();
+ r = prod (T (prod (a, b)), c) == prod (a, T (prod (b, c)));
+ r = prod (a + b, c) == prod (a, c) + prod (b, c);
+ r = prod (OneElement (value_type ()), a) == a;
+ r = prod (a, OneElement (value_type ())) == a;
+ ignore_unused_variable_warning (r);
+ }
+ };
+
+ template<class T>
+ struct CommutativeRingWithIdentityConcept {
+ typedef T value_type;
+
+ void constraints () {
+ function_requires< RingWithIdentityConcept<value_type> >();
+ bool r;
+ value_type a = value_type (), b = value_type ();
+ r = a * b == b * a;
+ ignore_unused_variable_warning (r);
+ }
+ };
+
+ template<class T>
+ struct FieldConcept {
+ typedef T value_type;
+
+ void constraints () {
+ function_requires< CommutativeRingWithIdentityConcept<value_type> >();
+ bool r;
+ value_type a = value_type ();
+ r = a == ZeroElement (value_type ()) || a * (OneElement (value_type ()) / a) == a;
+ r = a == ZeroElement (value_type ()) || (OneElement (value_type ()) / a) * a == a;
+ ignore_unused_variable_warning (r);
+ }
+ };
+
+ template<class T, class V>
+ struct VectorSpaceConcept {
+ typedef T value_type;
+ typedef V vector_type;
+
+ void constraints () {
+ function_requires< FieldConcept<value_type> >();
+ function_requires< AdditiveAbelianGroupConcept<vector_type> >();
+ bool r;
+ value_type alpha = value_type (), beta = value_type ();
+ vector_type a = vector_type (), b = vector_type ();
+ r = alpha * (a + b) == alpha * a + alpha * b;
+ r = (alpha + beta) * a == alpha * a + beta * a;
+ r = (alpha * beta) * a == alpha * (beta * a);
+ r = OneElement (value_type ()) * a == a;
+ ignore_unused_variable_warning (r);
+ }
+ };
+
+ template<class T, class V, class M>
+ struct LinearOperatorConcept {
+ typedef T value_type;
+ typedef V vector_type;
+ typedef M matrix_type;
+
+ void constraints () {
+ function_requires< VectorSpaceConcept<value_type, vector_type> >();
+ bool r;
+ value_type alpha = value_type (), beta = value_type ();
+ vector_type a = vector_type (), b = vector_type ();
+ matrix_type A = matrix_type ();
+ r = prod (A, alpha * a + beta * b) == alpha * prod (A, a) + beta * prod (A, b);
+ ignore_unused_variable_warning (r);
+ }
+ };
+
+inline void concept_checks () {
+
+ // Allow tests to be group to keep down compiler storage requirement
+#ifdef INTERAL
+#define INTERNAL_STORAGE
+#define INTERNAL_VECTOR
+#define INTERNAL_MATRIX
+#define INTERNAL_SPECIAL
+#define INTERNAL_SPARSE
+#define INTERNAL_EXPRESSION
+#endif
+
+ // TODO enable this for development
+ // #define VIEW_CONCEPTS
+
+ // Element value type for tests
+ typedef float T;
+
+ // Storage Array
+#if defined (INTERNAL_STORAGE) || defined (INTERNAL_STORAGE_DENSE)
+ {
+ typedef std::vector<T> container_model;
+ function_requires< Mutable_StorageArrayConcept<container_model> >();
+ function_requires< RandomAccessIteratorConcept<container_model::const_iterator> >();
+ function_requires< Mutable_RandomAccessIteratorConcept<container_model::iterator> >();
+ }
+
+ {
+ typedef bounded_array<T, 1> container_model;
+ function_requires< Mutable_StorageArrayConcept<container_model> >();
+ function_requires< RandomAccessIteratorConcept<container_model::const_iterator> >();
+ function_requires< Mutable_RandomAccessIteratorConcept<container_model::iterator> >();
+ }
+
+ {
+ typedef unbounded_array<T> container_model;
+ function_requires< Mutable_StorageArrayConcept<container_model> >();
+ function_requires< RandomAccessIteratorConcept<container_model::const_iterator> >();
+ function_requires< Mutable_RandomAccessIteratorConcept<container_model::iterator> >();
+ }
+
+/* FIXME array_adaptors are in progress
+ {
+ typedef array_adaptor<T> container_model;
+ function_requires< Mutable_StorageArrayConcept<container_model> >();
+ function_requires< RandomAccessIteratorConcept<container_model::const_iterator> >();
+ function_requires< Mutable_RandomAccessIteratorConcept<container_model::iterator> >();
+ }
+*/
+
+ {
+ typedef range container_model;
+ function_requires< IndexSetConcept<range> >();
+ function_requires< RandomAccessIteratorConcept<range::const_iterator> >();
+ }
+
+ {
+ typedef slice container_model;
+ function_requires< IndexSetConcept<range> >();
+ function_requires< RandomAccessIteratorConcept<range::const_iterator> >();
+ }
+
+ {
+ typedef indirect_array<> container_model;
+ function_requires< IndexSetConcept<range> >();
+ function_requires< RandomAccessIteratorConcept<range::const_iterator> >();
+ }
+#endif
+
+ // Storage Sparse
+#if defined (INTERNAL_STORAGE) || defined (INTERNAL_STORAGE_SPARSE)
+ {
+ typedef map_array<std::size_t, T> container_model;
+ function_requires< Mutable_StorageSparseConcept<container_model> >();
+ function_requires< RandomAccessIteratorConcept<container_model::const_iterator> >();
+ function_requires< RandomAccessIteratorConcept<container_model::iterator> >();
+ }
+
+ {
+ typedef std::map<std::size_t, T> container_model;
+ function_requires< Mutable_StorageSparseConcept<container_model > >();
+ function_requires< BidirectionalIteratorConcept<container_model::const_iterator> >();
+ function_requires< BidirectionalIteratorConcept<container_model::iterator> >();
+ }
+#endif
+
+#ifdef VIEW_CONCEPTS
+ // read only vectors
+ {
+ typedef vector_view<T> container_model;
+ function_requires< RandomAccessContainerConcept<container_model> >();
+ function_requires< VectorConcept<container_model> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
+ }
+#endif
+
+ // Vector
+#if defined (INTERNAL_VECTOR) || defined (INTERNAL_VECTOR_DENSE)
+ {
+ typedef vector<T> container_model;
+ function_requires< RandomAccessContainerConcept<container_model> >();
+ function_requires< Mutable_VectorConcept<container_model> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
+ function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::iterator> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
+ function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::reverse_iterator> >();
+ }
+
+ {
+ typedef zero_vector<T> container_model;
+ function_requires< VectorConcept<container_model> >();
+ function_requires< IndexedBidirectional1DIteratorConcept<container_model::const_iterator> >();
+ function_requires< IndexedBidirectional1DIteratorConcept<container_model::const_reverse_iterator> >();
+ }
+
+ {
+ typedef unit_vector<T> container_model;
+ function_requires< VectorConcept<container_model> >();
+ function_requires< IndexedBidirectional1DIteratorConcept<container_model::const_iterator> >();
+ function_requires< IndexedBidirectional1DIteratorConcept<container_model::const_reverse_iterator> >();
+ }
+
+ {
+ typedef scalar_vector<T> container_model;
+ function_requires< VectorConcept<container_model> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
+ }
+
+ {
+ typedef c_vector<T, 1> container_model;
+ function_requires< Mutable_VectorConcept<container_model> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
+ function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::iterator> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
+ function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::reverse_iterator> >();
+ }
+#endif
+
+ // Vector Proxies
+#if defined (INTERNAL_VECTOR) || defined (INTERNAL_VECTOR_PROXY)
+ {
+ typedef vector_range<vector<T> > container_model;
+ function_requires< Mutable_VectorExpressionConcept<container_model> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
+ function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::iterator> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
+ function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::reverse_iterator> >();
+ }
+
+ {
+ typedef vector_slice<vector<T> > container_model;
+ function_requires< Mutable_VectorExpressionConcept<container_model> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
+ function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::iterator> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
+ function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::reverse_iterator> >();
+ }
+
+ {
+ typedef vector_indirect<vector<T> > container_model;
+ function_requires< Mutable_VectorExpressionConcept<container_model> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
+ function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::iterator> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
+ function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::reverse_iterator> >();
+ }
+#endif
+
+ // Sparse Vector
+#if defined (INTERNAL_SPARSE) || defined (INTERNAL_VECTOR_SPARSE)
+ {
+ typedef mapped_vector<T> container_model;
+ function_requires< Mutable_SparseVectorConcept<container_model> >();
+ function_requires< IndexedBidirectional1DIteratorConcept<container_model::const_iterator> >();
+ function_requires< Mutable_IndexedBidirectional1DIteratorConcept<container_model::iterator> >();
+ function_requires< IndexedBidirectional1DIteratorConcept<container_model::const_reverse_iterator> >();
+ function_requires< Mutable_IndexedBidirectional1DIteratorConcept<container_model::reverse_iterator> >();
+ }
+
+ {
+ typedef compressed_vector<T> container_model;
+ function_requires< Mutable_SparseVectorConcept<container_model> >();
+ function_requires< IndexedBidirectional1DIteratorConcept<container_model::const_iterator> >();
+ function_requires< Mutable_IndexedBidirectional1DIteratorConcept<container_model::iterator> >();
+ function_requires< IndexedBidirectional1DIteratorConcept<container_model::const_reverse_iterator> >();
+ function_requires< Mutable_IndexedBidirectional1DIteratorConcept<container_model::reverse_iterator> >();
+ }
+
+ {
+ typedef coordinate_vector<T> container_model;
+ function_requires< Mutable_SparseVectorConcept<container_model> >();
+ function_requires< IndexedBidirectional1DIteratorConcept<container_model::const_iterator> >();
+ function_requires< Mutable_IndexedBidirectional1DIteratorConcept<container_model::iterator> >();
+ function_requires< IndexedBidirectional1DIteratorConcept<container_model::const_reverse_iterator> >();
+ function_requires< Mutable_IndexedBidirectional1DIteratorConcept<container_model::reverse_iterator> >();
+ }
+#endif
+
+ // Matrix
+#if defined (INTERNAL_MATRIX) || defined (INTERNAL_MATRIX_DENSE)
+ {
+ typedef matrix<T> container_model;
+ function_requires< Mutable_MatrixConcept<matrix<T> > >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
+ }
+
+ {
+ typedef vector_of_vector<T> container_model;
+ function_requires< Mutable_MatrixConcept<matrix<T> > >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
+ }
+
+ {
+ typedef zero_matrix<T> container_model;
+ function_requires< Mutable_MatrixConcept<matrix<T> > >();
+ function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
+ function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
+ }
+
+ {
+ typedef identity_matrix<T> container_model;
+ function_requires< Mutable_MatrixConcept<matrix<T> > >();
+ function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
+ function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
+ }
+
+ {
+ typedef scalar_matrix<T> container_model;
+ function_requires< Mutable_MatrixConcept<matrix<T> > >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
+ }
+
+ {
+ typedef c_matrix<T, 1, 1> container_model;
+ function_requires< Mutable_MatrixConcept<matrix<T> > >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
+ }
+#endif
+
+ // Matrix Proxies
+#if defined (INTERNAL_MATRIX) || defined (INTERNAL_MATRIX_PROXY)
+ {
+ typedef matrix_row<matrix<T> > container_model;
+ function_requires< Mutable_VectorExpressionConcept<container_model> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
+ function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::iterator> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
+ function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::reverse_iterator> >();
+ }
+
+ {
+ typedef matrix_column<matrix<T> > container_model;
+ function_requires< Mutable_VectorExpressionConcept<container_model> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
+ function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::iterator> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
+ function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::reverse_iterator> >();
+ }
+
+ {
+ typedef matrix_vector_range<matrix<T> > container_model;
+ function_requires< Mutable_VectorExpressionConcept<container_model> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
+ function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::iterator> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
+ function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::reverse_iterator> >();
+ }
+
+ {
+ typedef matrix_vector_slice<matrix<T> > container_model;
+ function_requires< Mutable_VectorExpressionConcept<container_model> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
+ function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::iterator> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
+ function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::reverse_iterator> >();
+ }
+
+ {
+ typedef matrix_vector_indirect<matrix<T> > container_model;
+ function_requires< Mutable_VectorExpressionConcept<container_model> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
+ function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::iterator> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
+ function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::reverse_iterator> >();
+ }
+
+ {
+ typedef matrix_range<matrix<T> > container_model;
+ function_requires< Mutable_MatrixExpressionConcept<container_model> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
+ }
+
+ {
+ typedef matrix_slice<matrix<T> > container_model;
+ function_requires< Mutable_MatrixExpressionConcept<container_model> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
+ }
+
+ {
+ typedef matrix_indirect<matrix<T> > container_model;
+ function_requires< Mutable_MatrixExpressionConcept<container_model> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
+ }
+#endif
+
+ // Banded Matrix
+#if defined (INTERNAL_SPECIAL) || defined (INTERNAL_BANDED)
+ {
+ typedef banded_matrix<T> container_model;
+ function_requires< Mutable_MatrixConcept<container_model> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
+ }
+
+ {
+ typedef banded_adaptor<matrix<T> > container_model;
+ function_requires< Mutable_MatrixExpressionConcept<container_model> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
+ }
+#endif
+
+ // Triangular Matrix
+#if defined (INTERNAL_SPECIAL) || defined (INTERNAL_TRIANGULAR)
+ {
+ typedef triangular_matrix<T> container_model;
+ function_requires< Mutable_MatrixConcept<container_model> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
+ }
+
+ {
+ typedef triangular_adaptor<matrix<T> > container_model;
+ function_requires< Mutable_MatrixExpressionConcept<container_model> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
+ }
+#endif
+
+ // Symmetric Matrix
+#if defined (INTERNA_SPECIAL) || defined (INTERNAL_SYMMETRIC)
+ {
+ typedef symmetric_matrix<T> container_model;
+ function_requires< Mutable_MatrixConcept<container_model> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
+ }
+
+ {
+ typedef banded_adaptor<matrix<T> > container_model;
+#ifndef SKIP_BAD
+ // const_iterator (iterator) constructor is bad
+ function_requires< Mutable_MatrixExpressionConcept<container_model> >();
+#endif
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
+ }
+#endif
+
+ // Hermitian Matrix
+#if defined (INTERNAL_SPECIAL) || defined (INTERNAL_HERMITIAN)
+ {
+ typedef hermitian_matrix<T> container_model;
+ function_requires< Mutable_MatrixConcept<container_model> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
+ }
+
+ {
+ typedef hermitian_adaptor<matrix<T> > container_model;
+#ifndef SKIP_BAD
+ // const_iterator (iterator) constructor is bad
+ function_requires< Mutable_MatrixExpressionConcept<container_model> >();
+#endif
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
+ }
+#endif
+
+ // Sparse Matrix
+#if defined (INTERNAL_SPARSE) || defined (INTERNAL_MATRIX_SPARSE)
+ {
+ typedef mapped_matrix<T> container_model;
+ function_requires< Mutable_SparseMatrixConcept<container_model> >();
+ function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
+ function_requires< Mutable_IndexedBidirectional2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
+ function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
+ function_requires< Mutable_IndexedBidirectional2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
+ }
+ {
+ typedef mapped_vector_of_mapped_vector<T> container_model;
+ function_requires< Mutable_SparseMatrixConcept<container_model> >();
+ function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
+ function_requires< Mutable_IndexedBidirectional2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
+ function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
+ function_requires< Mutable_IndexedBidirectional2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
+ }
+ {
+ typedef compressed_matrix<T> container_model;
+ function_requires< Mutable_SparseMatrixConcept<container_model> >();
+ function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
+ function_requires< Mutable_IndexedBidirectional2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
+ function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
+ function_requires< Mutable_IndexedBidirectional2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
+ }
+ {
+ typedef coordinate_matrix<T> container_model;
+ function_requires< Mutable_SparseMatrixConcept<container_model> >();
+ function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
+ function_requires< Mutable_IndexedBidirectional2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
+ function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
+ function_requires< Mutable_IndexedBidirectional2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
+ }
+ {
+ typedef generalized_vector_of_vector<T, row_major, vector< coordinate_vector<T> > > container_model;
+ function_requires< Mutable_SparseMatrixConcept<container_model> >();
+ function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
+ function_requires< Mutable_IndexedBidirectional2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
+ function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
+ function_requires< Mutable_IndexedBidirectional2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
+ }
+
+#endif
+
+ // Scalar Expressions
+#if defined (INTERNAL_EXPRESSION) || defined (INTERNAL_VECTOR_EXPRESSION)
+ function_requires< ScalarExpressionConcept<scalar_value<T> > >();
+ function_requires< ScalarExpressionConcept<scalar_reference<T> > >();
+
+ // Vector Expressions
+ {
+ typedef vector_reference<vector<T> > expression_model;
+ function_requires< VectorExpressionConcept<expression_model> >();
+ function_requires< Mutable_VectorExpressionConcept<expression_model> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_iterator> >();
+ function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<expression_model::iterator> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_reverse_iterator> >();
+ function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<expression_model::reverse_iterator> >();
+ }
+
+ {
+ typedef vector_unary<vector<T>, scalar_identity<T> > expression_model;
+ function_requires< VectorExpressionConcept<expression_model> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_iterator> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_reverse_iterator> >();
+ }
+
+ {
+ typedef vector_binary<vector<T>, vector<T>, scalar_plus<T, T> > expression_model;
+ function_requires< VectorExpressionConcept<expression_model> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_iterator> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_reverse_iterator> >();
+ }
+
+ {
+ typedef vector_binary_scalar1<T, vector<T>, scalar_multiplies<T, T> > expression_model;
+ function_requires< VectorExpressionConcept<expression_model> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_iterator> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_reverse_iterator> >();
+ }
+
+ {
+ typedef vector_binary_scalar2<vector<T>, scalar_value<T>, scalar_multiplies<T, T> > expression_model;
+ function_requires< VectorExpressionConcept<expression_model> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_iterator> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_reverse_iterator> >();
+ }
+
+ {
+ typedef vector_binary_scalar1<scalar_value<T>, vector<T>, scalar_multiplies<T, T> > expression_model;
+ function_requires< VectorExpressionConcept<expression_model> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_iterator> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_reverse_iterator> >();
+ }
+
+ {
+ typedef vector_binary_scalar2<vector<T>, scalar_value<T>, scalar_multiplies<T, T> > expression_model;
+ function_requires< VectorExpressionConcept<expression_model> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_iterator> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_reverse_iterator> >();
+ }
+
+ function_requires< ScalarExpressionConcept<vector_scalar_unary<vector<T>, vector_sum<vector<T> > > > >();
+ function_requires< ScalarExpressionConcept<vector_scalar_unary<vector<T>, vector_norm_1<vector<T> > > > >();
+ function_requires< ScalarExpressionConcept<vector_scalar_unary<vector<T>, vector_norm_2<vector<T> > > > >();
+ function_requires< ScalarExpressionConcept<vector_scalar_unary<vector<T>, vector_norm_inf<vector<T> > > > >();
+
+ function_requires< ScalarExpressionConcept<vector_scalar_binary<vector<T>, vector<T>, vector_inner_prod<vector<T>, vector<T>, T> > > >();
+#endif
+
+ // Matrix Expressions
+#if defined (INTERNAL_EXPRESSION) || defined (INTERNAL_MATRIX_EXPRESSION)
+ {
+ typedef matrix_reference<matrix<T> > expression_model;
+ function_requires< MatrixExpressionConcept<expression_model> >();
+ function_requires< Mutable_MatrixExpressionConcept<expression_model> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_iterator1, expression_model::const_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<expression_model::iterator1, expression_model::iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_reverse_iterator1, expression_model::const_reverse_iterator2> >();
+ function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<expression_model::reverse_iterator1, expression_model::reverse_iterator2> >();
+ }
+
+ {
+ typedef vector_matrix_binary<vector<T>, vector<T>, scalar_multiplies<T, T> > expression_model;
+ function_requires< MatrixExpressionConcept<expression_model> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_iterator1, expression_model::const_iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_reverse_iterator1, expression_model::const_reverse_iterator2> >();
+ }
+
+ {
+ typedef matrix_unary1<matrix<T>, scalar_identity<T> > expression_model;
+ function_requires< MatrixExpressionConcept<expression_model> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_iterator1, expression_model::const_iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_reverse_iterator1, expression_model::const_reverse_iterator2> >();
+ }
+
+ {
+ typedef matrix_unary2<matrix<T>, scalar_identity<T> > expression_model;
+ function_requires< MatrixExpressionConcept<expression_model> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_iterator1, expression_model::const_iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_reverse_iterator1, expression_model::const_reverse_iterator2> >();
+ }
+
+ {
+ typedef matrix_binary<matrix<T>, matrix<T>, scalar_plus<T, T> > expression_model;
+ function_requires< MatrixExpressionConcept<expression_model> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_iterator1, expression_model::const_iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_reverse_iterator1, expression_model::const_reverse_iterator2> >();
+ }
+
+ {
+ typedef matrix_binary_scalar1<T, matrix<T>, scalar_multiplies<T, T> > expression_model;
+ function_requires< MatrixExpressionConcept<expression_model> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_iterator1, expression_model::const_iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_reverse_iterator1, expression_model::const_reverse_iterator2> >();
+ }
+
+ {
+ typedef matrix_binary_scalar2<matrix<T>, T, scalar_multiplies<T, T> > expression_model;
+ function_requires< MatrixExpressionConcept<expression_model> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_iterator1, expression_model::const_iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_reverse_iterator1, expression_model::const_reverse_iterator2> >();
+ }
+
+ {
+ typedef matrix_binary_scalar1<scalar_value<T>, matrix<T>, scalar_multiplies<T, T> > expression_model;
+ function_requires< MatrixExpressionConcept<expression_model> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_iterator1, expression_model::const_iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_reverse_iterator1, expression_model::const_reverse_iterator2> >();
+ }
+
+ {
+ typedef matrix_binary_scalar2<matrix<T>, scalar_value<T>, scalar_multiplies<T, T> > expression_model;
+ function_requires< MatrixExpressionConcept<expression_model> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_iterator1, expression_model::const_iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_reverse_iterator1, expression_model::const_reverse_iterator2> >();
+ }
+
+ {
+ typedef matrix_vector_binary1<matrix<T>, vector<T>, matrix_vector_prod1<matrix<T>, vector<T>, T> > expression_model;
+ function_requires< VectorExpressionConcept<expression_model> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_iterator> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_reverse_iterator> >();
+ }
+
+ {
+ typedef matrix_vector_binary2<vector<T>, matrix<T>, matrix_vector_prod2<matrix<T>, vector<T>, T > > expression_model;
+ function_requires< VectorExpressionConcept<expression_model> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_iterator> >();
+ function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_reverse_iterator> >();
+ }
+
+ {
+ typedef matrix_matrix_binary<matrix<T>, matrix<T>, matrix_matrix_prod<matrix<T>, matrix<T>, T > > expression_model;
+ function_requires< MatrixExpressionConcept<expression_model> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_iterator1, expression_model::const_iterator2> >();
+ function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_reverse_iterator1, expression_model::const_reverse_iterator2> >();
+ }
+
+ function_requires< ScalarExpressionConcept<matrix_scalar_unary<matrix<T>, matrix_norm_1<vector<T> > > > >();
+ function_requires< ScalarExpressionConcept<matrix_scalar_unary<matrix<T>, matrix_norm_frobenius<vector<T> > > > >();
+ function_requires< ScalarExpressionConcept<matrix_scalar_unary<matrix<T>, matrix_norm_inf<vector<T> > > > >();
+#endif
+
+#ifdef EXTERNAL
+ function_requires< AdditiveAbelianGroupConcept<T> >();
+ function_requires< CommutativeRingWithIdentityConcept<T> >();
+ function_requires< FieldConcept<T> >();
+ function_requires< VectorSpaceConcept<T, vector<T> > >();
+ function_requires< Prod_RingWithIdentityConcept<matrix<T> > >();
+ function_requires< VectorSpaceConcept<T, matrix<T> > >();
+ function_requires< LinearOperatorConcept<T, vector<T>, matrix<T> > >();
+
+ function_requires< AdditiveAbelianGroupConcept<std::complex<T> > >();
+ function_requires< CommutativeRingWithIdentityConcept<std::complex<T> > >();
+ function_requires< FieldConcept<std::complex<T> > >();
+ function_requires< VectorSpaceConcept<std::complex<T>, vector<std::complex<T> > > >();
+ function_requires< Prod_RingWithIdentityConcept<matrix<std::complex<T> > > >();
+ function_requires< VectorSpaceConcept<std::complex<T>, matrix<std::complex<T> > > >();
+ function_requires< LinearOperatorConcept<std::complex<T>, vector<std::complex<T> >, matrix<std::complex<T> > > >();
+#endif
+ }
+
+ } // end of anonymous namespace
+
+}}}
+
+#endif
diff --git a/include/boost/numeric/ublas/detail/config.hpp b/include/boost/numeric/ublas/detail/config.hpp
new file mode 100644
index 0000000..9e67410
--- /dev/null
+++ b/include/boost/numeric/ublas/detail/config.hpp
@@ -0,0 +1,304 @@
+//
+// Copyright (c) 2000-2002
+// Joerg Walter, Mathias Koch
+//
+// 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)
+//
+// The authors gratefully acknowledge the support of
+// GeNeSys mbH & Co. KG in producing this work.
+//
+
+#ifndef _BOOST_UBLAS_CONFIG_
+#define _BOOST_UBLAS_CONFIG_
+
+#include <cassert>
+#include <cstddef>
+#include <algorithm>
+#include <limits>
+
+#include <boost/config.hpp>
+#include <boost/static_assert.hpp>
+#include <boost/noncopyable.hpp>
+#include <boost/mpl/if.hpp>
+#include <boost/mpl/and.hpp>
+#include <boost/type_traits/is_same.hpp>
+#include <boost/type_traits/is_convertible.hpp>
+#include <boost/type_traits/is_const.hpp>
+#include <boost/type_traits/remove_reference.hpp>
+
+// C++11
+#if defined(__cplusplus) && __cplusplus >= 201103L
+
+#define BOOST_UBLAS_CPP_GE_2011
+
+#elif BOOST_MSVC >= 1800
+
+#define BOOST_UBLAS_CPP_GE_2011
+
+#else
+
+#undef BOOST_UBLAS_CPP_GE_2011 // Make sure no one defined it
+
+#endif
+
+// Microsoft Visual C++
+#if defined (BOOST_MSVC) && ! defined (BOOST_STRICT_CONFIG)
+
+// Version 7.1
+#if BOOST_MSVC == 1310
+// One of these workarounds is needed for MSVC 7.1 AFAIK
+// (thanks to John Maddock and Martin Lauer).
+#if !(defined(BOOST_UBLAS_NO_NESTED_CLASS_RELATION) || defined(BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION))
+#define BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+#endif
+
+#endif
+
+#endif
+
+
+// GNU Compiler Collection
+#if defined (__GNUC__) && ! defined (BOOST_STRICT_CONFIG)
+
+#if __GNUC__ >= 4 || (__GNUC__ >= 3 && __GNUC_MINOR__ >= 4)
+// Specified by ABI definition see GCC bug id 9982
+#define BOOST_UBLAS_USEFUL_ARRAY_PLACEMENT_NEW
+#endif
+
+#if __GNUC__ < 3
+#define BOOST_UBLAS_UNSUPPORTED_COMPILER 1
+#endif
+
+#endif
+
+
+// Intel Compiler
+#if defined (BOOST_INTEL) && ! defined (BOOST_STRICT_CONFIG)
+
+#if defined (BOOST_INTEL_LINUX) && (BOOST_INTEL_LINUX >= 800)
+// By inspection of compiler results
+#define BOOST_UBLAS_USEFUL_ARRAY_PLACEMENT_NEW
+#endif
+
+#if (BOOST_INTEL < 700)
+#define BOOST_UBLAS_UNSUPPORTED_COMPILER 1
+#endif
+
+// Define swap for index_pair and triple.
+#if (BOOST_INTEL <= 800)
+namespace boost { namespace numeric { namespace ublas {
+ template<class C, class IC>
+ class indexed_iterator;
+
+ template<class V>
+ class index_pair;
+ template<class M>
+ class index_triple;
+}}}
+
+namespace std {
+ template<class V>
+ inline
+ void swap (boost::numeric::ublas::index_pair<V> i1, boost::numeric::ublas::index_pair<V> i2) {
+ i1.swap (i2);
+ }
+ template<class M>
+ inline
+ void swap (boost::numeric::ublas::index_triple<M> i1, boost::numeric::ublas::index_triple<M> i2) {
+ i1.swap (i2);
+ }
+ // iter_swap also needed for ICC on Itanium?
+ template<class C, class IC>
+ inline
+ void iter_swap (boost::numeric::ublas::indexed_iterator<C, IC> it1,
+ boost::numeric::ublas::indexed_iterator<C, IC> it2) {
+ swap (*it1, *it2);
+ }
+}
+#endif
+
+#endif
+
+
+// Comeau compiler - thanks to Kresimir Fresl
+#if defined (__COMO__) && ! defined (BOOST_STRICT_CONFIG)
+
+// Missing std::abs overloads for float types in <cmath> are in <cstdlib>
+#if defined(__LIBCOMO__) && (__LIBCOMO_VERSION__ <= 31)
+#include <cstdlib>
+#endif
+
+#endif
+
+// PGI compiler
+#ifdef __PGIC__
+#define BOOST_UBLAS_UNSUPPORTED_COMPILER 0
+#endif
+
+// HP aCC C++ compiler
+#if defined (__HP_aCC) && ! defined (BOOST_STRICT_CONFIG)
+# if (__HP_aCC >= 60000 )
+# define BOOST_UBLAS_USEFUL_ARRAY_PLACEMENT_NEW
+#endif
+#endif
+
+
+// SGI MIPSpro C++ compiler
+#if defined (__sgi) && ! defined (BOOST_STRICT_CONFIG)
+
+// Missing std::abs overloads for float types in <cmath> are in <cstdlib>
+// This test should be library version specific.
+#include <cstdlib>
+
+#if __COMPILER_VERSION >=650
+// By inspection of compiler results - thanks to Peter Schmitteckert
+#define BOOST_UBLAS_USEFUL_ARRAY_PLACEMENT_NEW
+#endif
+
+#endif
+
+
+// Metrowerks Codewarrior
+#if defined (__MWERKS__) && ! defined (BOOST_STRICT_CONFIG)
+
+// 8.x
+#if __MWERKS__ <= 0x3003
+#define BOOST_UBLAS_UNSUPPORTED_COMPILER 1
+#endif
+
+#endif
+
+
+// Detect other compilers with serious defects - override by defineing BOOST_UBLAS_UNSUPPORTED_COMPILER=0
+#ifndef BOOST_UBLAS_UNSUPPORTED_COMPILER
+#if defined(BOOST_NO_FUNCTION_TEMPLATE_ORDERING) || defined(BOOST_NO_SFINAE) || defined(BOOST_NO_STDC_NAMESPACE)
+#define BOOST_UBLAS_UNSUPPORTED_COMPILER 1
+#endif
+#endif
+
+// Cannot continue with an unsupported compiler
+#if defined(BOOST_UBLAS_UNSUPPORTED_COMPILER) && (BOOST_UBLAS_UNSUPPORTED_COMPILER != 0)
+#error Your compiler and/or configuration is unsupported by this verions of uBLAS. Define BOOST_UBLAS_UNSUPPORTED_COMPILER=0 to override this message. Boost 1.32.0 includes uBLAS with support for many older compilers.
+#endif
+
+
+
+// Enable performance options in RELEASE mode
+#if defined (NDEBUG) || defined (BOOST_UBLAS_NDEBUG)
+
+#ifndef BOOST_UBLAS_INLINE
+#define BOOST_UBLAS_INLINE inline
+#endif
+
+// Do not check sizes!
+#define BOOST_UBLAS_USE_FAST_SAME
+
+// NO runtime error checks with BOOST_UBLAS_CHECK macro
+#ifndef BOOST_UBLAS_CHECK_ENABLE
+#define BOOST_UBLAS_CHECK_ENABLE 0
+#endif
+
+// NO type compatibility numeric checks
+#ifndef BOOST_UBLAS_TYPE_CHECK
+#define BOOST_UBLAS_TYPE_CHECK 0
+#endif
+
+
+// Disable performance options in DEBUG mode
+#else
+
+#ifndef BOOST_UBLAS_INLINE
+#define BOOST_UBLAS_INLINE
+#endif
+
+// Enable runtime error checks with BOOST_UBLAS_CHECK macro. Check bounds etc
+#ifndef BOOST_UBLAS_CHECK_ENABLE
+#define BOOST_UBLAS_CHECK_ENABLE 1
+#endif
+
+// Type compatibiltity numeric checks
+#ifndef BOOST_UBLAS_TYPE_CHECK
+#define BOOST_UBLAS_TYPE_CHECK 1
+#endif
+
+#endif
+
+
+/*
+ * Type compatibility checks
+ * Control type compatibility numeric runtime checks for non dense matrices.
+ * Require additional storage and complexity
+ */
+#if BOOST_UBLAS_TYPE_CHECK
+template <class Dummy>
+struct disable_type_check
+{
+ static bool value;
+};
+template <class Dummy>
+bool disable_type_check<Dummy>::value = false;
+#endif
+#ifndef BOOST_UBLAS_TYPE_CHECK_EPSILON
+#define BOOST_UBLAS_TYPE_CHECK_EPSILON (type_traits<real_type>::type_sqrt (std::numeric_limits<real_type>::epsilon ()))
+#endif
+#ifndef BOOST_UBLAS_TYPE_CHECK_MIN
+#define BOOST_UBLAS_TYPE_CHECK_MIN (type_traits<real_type>::type_sqrt ( (std::numeric_limits<real_type>::min) ()))
+#endif
+
+
+/*
+ * General Configuration
+ */
+
+// Proxy shortcuts overload the alreadly heavily over used operator ()
+//#define BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
+
+// In order to simplify debugging is is possible to simplify expression template
+// so they are restricted to a single operation
+// #define BOOST_UBLAS_SIMPLE_ET_DEBUG
+
+// Use invariant hoisting.
+// #define BOOST_UBLAS_USE_INVARIANT_HOISTING
+
+// Use Duff's device in element access loops
+// #define BOOST_UBLAS_USE_DUFF_DEVICE
+
+// Choose evaluation method for dense vectors and matrices
+#if !(defined(BOOST_UBLAS_USE_INDEXING) || defined(BOOST_UBLAS_USE_ITERATING))
+#define BOOST_UBLAS_USE_INDEXING
+#endif
+// #define BOOST_UBLAS_ITERATOR_THRESHOLD 0
+
+// Use indexed iterators - unsupported implementation experiment
+// #define BOOST_UBLAS_USE_INDEXED_ITERATOR
+
+// Alignment of bounded_array type
+#ifndef BOOST_UBLAS_BOUNDED_ARRAY_ALIGN
+#define BOOST_UBLAS_BOUNDED_ARRAY_ALIGN
+#endif
+
+// Enable different sparse element proxies
+#ifndef BOOST_UBLAS_NO_ELEMENT_PROXIES
+// Sparse proxies prevent reference invalidation problems in expressions such as:
+// a [1] = a [0] = 1 Thanks to Marc Duflot for spotting this.
+// #define BOOST_UBLAS_STRICT_MAP_ARRAY
+#define BOOST_UBLAS_STRICT_VECTOR_SPARSE
+#define BOOST_UBLAS_STRICT_MATRIX_SPARSE
+// Hermitian matrices use element proxies to allow assignment to conjugate triangle
+#define BOOST_UBLAS_STRICT_HERMITIAN
+#endif
+
+// Define to configure special settings for reference returning members
+// #define BOOST_UBLAS_REFERENCE_CONST_MEMBER
+// #define BOOST_UBLAS_PROXY_CONST_MEMBER
+
+
+// Include type declerations and functions
+#include <boost/numeric/ublas/fwd.hpp>
+#include <boost/numeric/ublas/detail/definitions.hpp>
+
+
+#endif
+
diff --git a/include/boost/numeric/ublas/detail/definitions.hpp b/include/boost/numeric/ublas/detail/definitions.hpp
new file mode 100644
index 0000000..c5e1cfc
--- /dev/null
+++ b/include/boost/numeric/ublas/detail/definitions.hpp
@@ -0,0 +1,212 @@
+//
+// Copyright (c) 2000-2002
+// Joerg Walter, Mathias Koch
+//
+// 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)
+//
+// The authors gratefully acknowledge the support of
+// GeNeSys mbH & Co. KG in producing this work.
+//
+
+#ifndef _BOOST_UBLAS_DEFINITIONS_
+#define _BOOST_UBLAS_DEFINITIONS_
+
+
+namespace boost { namespace numeric { namespace ublas {
+
+ namespace detail {
+ /* Borrowed from boost/concept_checks.hpp
+ "inline" is used for ignore_unused_variable_warning()
+ to make sure there is no overhead with g++.
+ */
+ template <class T> inline
+ void ignore_unused_variable_warning(const T&) {}
+ } // namespace detail
+
+ // Borrowed from Dave Abraham's noncopyable.
+ // I believe this should be part of utility.hpp one day...
+ namespace nonassignable_ // protection from unintended ADL
+ {
+ class nonassignable {
+ protected:
+ nonassignable () {}
+ ~nonassignable () {}
+ private: // emphasize the following members are private
+ const nonassignable& operator= (const nonassignable &);
+ }; // nonassignable
+ }
+ typedef nonassignable_::nonassignable nonassignable;
+
+
+ // Assignment proxy.
+ // Provides temporary free assigment when LHS has no alias on RHS
+ template<class C>
+ class noalias_proxy:
+ private nonassignable {
+ public:
+ typedef typename C::closure_type closure_type;
+
+ BOOST_UBLAS_INLINE
+ noalias_proxy (C& lval):
+ nonassignable (), lval_ (lval) {}
+ BOOST_UBLAS_INLINE
+ noalias_proxy (const noalias_proxy& p):
+ nonassignable (), lval_ (p.lval_) {}
+
+ template <class E>
+ BOOST_UBLAS_INLINE
+ closure_type &operator= (const E& e) {
+ lval_.assign (e);
+ return lval_;
+ }
+
+ template <class E>
+ BOOST_UBLAS_INLINE
+ closure_type &operator+= (const E& e) {
+ lval_.plus_assign (e);
+ return lval_;
+ }
+
+ template <class E>
+ BOOST_UBLAS_INLINE
+ closure_type &operator-= (const E& e) {
+ lval_.minus_assign (e);
+ return lval_;
+ }
+
+ private:
+ closure_type lval_;
+ };
+
+ // Improve syntax of efficient assignment where no aliases of LHS appear on the RHS
+ // noalias(lhs) = rhs_expression
+ template <class C>
+ BOOST_UBLAS_INLINE
+ noalias_proxy<C> noalias (C& lvalue) {
+ return noalias_proxy<C> (lvalue);
+ }
+ template <class C>
+ BOOST_UBLAS_INLINE
+ noalias_proxy<const C> noalias (const C& lvalue) {
+ return noalias_proxy<const C> (lvalue);
+ }
+
+ // Possible future compatible syntax where lvalue possible has an unsafe alias on the RHS
+ // safe(lhs) = rhs_expression
+ template <class C>
+ BOOST_UBLAS_INLINE
+ C& safe (C& lvalue) {
+ return lvalue;
+ }
+ template <class C>
+ BOOST_UBLAS_INLINE
+ const C& safe (const C& lvalue) {
+ return lvalue;
+ }
+
+
+ // Dimension accessors
+ namespace dimension {
+
+ // Generic accessors
+ template<unsigned dimension>
+ struct dimension_properties {};
+
+ template<>
+ struct dimension_properties<1> {
+ template <class E>
+ BOOST_UBLAS_INLINE static
+ typename E::size_type size (const vector_expression<E> &e) {
+ return e ().size ();
+ }
+ template <class E>
+ BOOST_UBLAS_INLINE static
+ typename E::size_type size (const matrix_expression<E> &e) {
+ return e ().size1 ();
+ }
+ // Note: Index functions cannot deduce dependant template parameter V or M from i
+ template <class V>
+ BOOST_UBLAS_INLINE static
+ typename V::size_type index (const typename V::iterator &i) {
+ return i.index ();
+ }
+ template <class M>
+ BOOST_UBLAS_INLINE static
+ typename M::size_type index (const typename M::iterator1 &i) {
+ return i.index1 ();
+ }
+ template <class M>
+ BOOST_UBLAS_INLINE static
+ typename M::size_type index (const typename M::iterator2 &i) {
+ return i.index1 ();
+ }
+ };
+ template<>
+ struct dimension_properties<2> {
+ template <class E>
+ BOOST_UBLAS_INLINE static
+ typename E::size_type size (const vector_expression<E> &) {
+ return 1;
+ }
+ template <class E>
+ BOOST_UBLAS_INLINE static
+ typename E::size_type size (const matrix_expression<E> &e) {
+ return e ().size2 ();
+ }
+ template <class V>
+ BOOST_UBLAS_INLINE static
+ typename V::size_type index (const typename V::iterator &) {
+ return 1;
+ }
+ template <class M>
+ BOOST_UBLAS_INLINE static
+ typename M::size_type index (const typename M::iterator1 &i) {
+ return i.index2 ();
+ }
+ template <class M>
+ BOOST_UBLAS_INLINE static
+ typename M::size_type index (const typename M::iterator2 &i) {
+ return i.index2 ();
+ }
+ };
+
+ template<unsigned dimension, class E>
+ BOOST_UBLAS_INLINE
+ typename E::size_type size (const E& e) {
+ return dimension_properties<dimension>::size (e);
+ }
+
+ template<unsigned dimension, class I>
+ BOOST_UBLAS_INLINE
+ typename I::container_type::size_type
+ index (const I& i) {
+ typedef typename I::container_type container_type;
+ return dimension_properties<dimension>::template index<container_type> (i);
+ }
+
+
+ // Named accessors - just syntactic sugar
+ template<class V>
+ typename V::size_type num_elements (const V &v) {
+ return v.size ();
+ }
+ template<class M>
+ typename M::size_type num_rows (const M &m) {
+ return m.size1 ();
+ }
+ template<class M>
+ typename M::size_type num_columns (const M &m) {
+ return m.size2 ();
+ }
+ template<class MV>
+ typename MV::size_type num_non_zeros (const MV &mv) {
+ return mv.non_zeros ();
+ }
+ }
+
+
+}}}
+
+#endif
diff --git a/include/boost/numeric/ublas/detail/documentation.hpp b/include/boost/numeric/ublas/detail/documentation.hpp
new file mode 100644
index 0000000..4b2bcf0
--- /dev/null
+++ b/include/boost/numeric/ublas/detail/documentation.hpp
@@ -0,0 +1,33 @@
+//
+// Copyright (c) 2000-2004
+// Joerg Walter, Mathias Koch
+//
+// 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)
+//
+// The authors gratefully acknowledge the support of
+// GeNeSys mbH & Co. KG in producing this work.
+//
+
+// this file should not contain any code, but the documentation
+// global to all files
+
+/** \namespace boost::numeric::ublas
+ \brief contains all important classes and functions of uBLAS
+
+ all ublas definitions ...
+ \todo expand this section
+ */
+
+/** \defgroup blas1 Level 1 BLAS
+ \brief level 1 basic linear algebra subroutines
+*/
+
+/** \defgroup blas2 Level 2 BLAS
+ \brief level 2 basic linear algebra subroutines
+*/
+
+/** \defgroup blas3 Level 3 BLAS
+ \brief level 3 basic linear algebra subroutines
+*/
diff --git a/include/boost/numeric/ublas/detail/duff.hpp b/include/boost/numeric/ublas/detail/duff.hpp
new file mode 100644
index 0000000..b0ec08c
--- /dev/null
+++ b/include/boost/numeric/ublas/detail/duff.hpp
@@ -0,0 +1,56 @@
+//
+// Copyright (c) 2000-2002
+// Joerg Walter, Mathias Koch
+//
+// 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)
+//
+// The authors gratefully acknowledge the support of
+// GeNeSys mbH & Co. KG in producing this work.
+//
+
+#ifndef _BOOST_UBLAS_DUFF_
+#define _BOOST_UBLAS_DUFF_
+
+#define DD_SWITCH(n, d, r, expr) \
+ { \
+ unsigned r = ((n) + (d) - 1) / (d); \
+ switch ((n) % (d)) { \
+ case 0: do { expr;
+#define DD_CASE_I(i, expr) \
+ case (i): expr;
+#define DD_WHILE(r) \
+ } while (-- (r) > 0); \
+ } \
+ }
+
+#define DD_1T(n, d, r, expr) \
+ DD_WHILE(r)
+#define DD_2T(n, d, r, expr) \
+ DD_CASE_I(1, expr) \
+ DD_1T(n, d, r, expr)
+#define DD_3T(n, d, r, expr) \
+ DD_CASE_I(2, expr) \
+ DD_2T(n, d, r, expr)
+#define DD_4T(n, d, r, expr) \
+ DD_CASE_I(3, expr) \
+ DD_3T(n, d, r, expr)
+#define DD_5T(n, d, r, expr) \
+ DD_CASE_I(4, expr) \
+ DD_4T(n, d, r, expr)
+#define DD_6T(n, d, r, expr) \
+ DD_CASE_I(5, expr) \
+ DD_5T(n, d, r, expr)
+#define DD_7T(n, d, r, expr) \
+ DD_CASE_I(6, expr) \
+ DD_6T(n, d, r, expr)
+#define DD_8T(n, d, r, expr) \
+ DD_CASE_I(7, expr) \
+ DD_7T(n, d, r, expr)
+
+#define DD(n, d, r, expr) \
+ DD_SWITCH(n, d, r, expr) \
+ DD_##d##T(n, d, r, expr)
+
+#endif
diff --git a/include/boost/numeric/ublas/detail/iterator.hpp b/include/boost/numeric/ublas/detail/iterator.hpp
new file mode 100644
index 0000000..1723a30
--- /dev/null
+++ b/include/boost/numeric/ublas/detail/iterator.hpp
@@ -0,0 +1,1436 @@
+//
+// Copyright (c) 2000-2002
+// Joerg Walter, Mathias Koch
+//
+// 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)
+//
+// The authors gratefully acknowledge the support of
+// GeNeSys mbH & Co. KG in producing this work.
+//
+
+#ifndef _BOOST_UBLAS_ITERATOR_
+#define _BOOST_UBLAS_ITERATOR_
+
+#include <boost/numeric/ublas/exception.hpp>
+#include <iterator>
+
+
+namespace boost { namespace numeric { namespace ublas {
+
+ /** \brief Base class of all proxy classes that contain
+ * a (redirectable) reference to an immutable object.
+ *
+ * \param C the type of the container referred to
+ */
+ template<class C>
+ class container_const_reference:
+ private nonassignable {
+ public:
+ typedef C container_type;
+
+ BOOST_UBLAS_INLINE
+ container_const_reference ():
+ c_ (0) {}
+ BOOST_UBLAS_INLINE
+ container_const_reference (const container_type &c):
+ c_ (&c) {}
+
+ BOOST_UBLAS_INLINE
+ const container_type &operator () () const {
+ return *c_;
+ }
+
+ BOOST_UBLAS_INLINE
+ container_const_reference &assign (const container_type *c) {
+ c_ = c;
+ return *this;
+ }
+
+ // Closure comparison
+ BOOST_UBLAS_INLINE
+ bool same_closure (const container_const_reference &cr) const {
+ return c_ == cr.c_;
+ }
+
+ private:
+ const container_type *c_;
+ };
+
+ /** \brief Base class of all proxy classes that contain
+ * a (redirectable) reference to a mutable object.
+ *
+ * \param C the type of the container referred to
+ */
+ template<class C>
+ class container_reference:
+ private nonassignable {
+ public:
+ typedef C container_type;
+
+ BOOST_UBLAS_INLINE
+ container_reference ():
+ c_ (0) {}
+ BOOST_UBLAS_INLINE
+ container_reference (container_type &c):
+ c_ (&c) {}
+
+ BOOST_UBLAS_INLINE
+ container_type &operator () () const {
+ return *c_;
+ }
+
+ BOOST_UBLAS_INLINE
+ container_reference &assign (container_type *c) {
+ c_ = c;
+ return *this;
+ }
+
+ // Closure comparison
+ BOOST_UBLAS_INLINE
+ bool same_closure (const container_reference &cr) const {
+ return c_ == cr.c_;
+ }
+
+ private:
+ container_type *c_;
+ };
+
+ /** \brief Base class of all forward iterators.
+ *
+ * \param IC the iterator category
+ * \param I the derived iterator type
+ * \param T the value type
+ *
+ * The forward iterator can only proceed in one direction
+ * via the post increment operator.
+ */
+ template<class IC, class I, class T>
+ struct forward_iterator_base:
+ public std::iterator<IC, T> {
+ typedef I derived_iterator_type;
+ typedef T derived_value_type;
+
+ // Arithmetic
+ BOOST_UBLAS_INLINE
+ derived_iterator_type operator ++ (int) {
+ derived_iterator_type &d (*static_cast<const derived_iterator_type *> (this));
+ derived_iterator_type tmp (d);
+ ++ d;
+ return tmp;
+ }
+ BOOST_UBLAS_INLINE
+ friend derived_iterator_type operator ++ (derived_iterator_type &d, int) {
+ derived_iterator_type tmp (d);
+ ++ d;
+ return tmp;
+ }
+
+ // Comparison
+ BOOST_UBLAS_INLINE
+ bool operator != (const derived_iterator_type &it) const {
+ const derived_iterator_type *d = static_cast<const derived_iterator_type *> (this);
+ return ! (*d == it);
+ }
+ };
+
+ /** \brief Base class of all bidirectional iterators.
+ *
+ * \param IC the iterator category
+ * \param I the derived iterator type
+ * \param T the value type
+ *
+ * The bidirectional iterator can proceed in both directions
+ * via the post increment and post decrement operator.
+ */
+ template<class IC, class I, class T>
+ struct bidirectional_iterator_base:
+ public std::iterator<IC, T> {
+ typedef I derived_iterator_type;
+ typedef T derived_value_type;
+
+ // Arithmetic
+ BOOST_UBLAS_INLINE
+ derived_iterator_type operator ++ (int) {
+ derived_iterator_type &d (*static_cast<const derived_iterator_type *> (this));
+ derived_iterator_type tmp (d);
+ ++ d;
+ return tmp;
+ }
+ BOOST_UBLAS_INLINE
+ friend derived_iterator_type operator ++ (derived_iterator_type &d, int) {
+ derived_iterator_type tmp (d);
+ ++ d;
+ return tmp;
+ }
+ BOOST_UBLAS_INLINE
+ derived_iterator_type operator -- (int) {
+ derived_iterator_type &d (*static_cast<const derived_iterator_type *> (this));
+ derived_iterator_type tmp (d);
+ -- d;
+ return tmp;
+ }
+ BOOST_UBLAS_INLINE
+ friend derived_iterator_type operator -- (derived_iterator_type &d, int) {
+ derived_iterator_type tmp (d);
+ -- d;
+ return tmp;
+ }
+
+ // Comparison
+ BOOST_UBLAS_INLINE
+ bool operator != (const derived_iterator_type &it) const {
+ const derived_iterator_type *d = static_cast<const derived_iterator_type *> (this);
+ return ! (*d == it);
+ }
+ };
+
+ /** \brief Base class of all random access iterators.
+ *
+ * \param IC the iterator category
+ * \param I the derived iterator type
+ * \param T the value type
+ * \param D the difference type, default: std::ptrdiff_t
+ *
+ * The random access iterator can proceed in both directions
+ * via the post increment/decrement operator or in larger steps
+ * via the +, - and +=, -= operators. The random access iterator
+ * is LessThan Comparable.
+ */
+ template<class IC, class I, class T, class D = std::ptrdiff_t>
+ // ISSUE the default for D seems rather dangerous as it can easily be (silently) incorrect
+ struct random_access_iterator_base:
+ public std::iterator<IC, T> {
+ typedef I derived_iterator_type;
+ typedef T derived_value_type;
+ typedef D derived_difference_type;
+
+ /* FIXME Need to explicitly pass derived_reference_type as otherwise I undefined type or forward declared
+ typedef typename derived_iterator_type::reference derived_reference_type;
+ // Indexed element
+ BOOST_UBLAS_INLINE
+ derived_reference_type operator [] (derived_difference_type n) {
+ return *(*this + n);
+ }
+ */
+
+ // Arithmetic
+ BOOST_UBLAS_INLINE
+ derived_iterator_type operator ++ (int) {
+ derived_iterator_type &d (*static_cast<derived_iterator_type *> (this));
+ derived_iterator_type tmp (d);
+ ++ d;
+ return tmp;
+ }
+ BOOST_UBLAS_INLINE
+ friend derived_iterator_type operator ++ (derived_iterator_type &d, int) {
+ derived_iterator_type tmp (d);
+ ++ d;
+ return tmp;
+ }
+ BOOST_UBLAS_INLINE
+ derived_iterator_type operator -- (int) {
+ derived_iterator_type &d (*static_cast<derived_iterator_type *> (this));
+ derived_iterator_type tmp (d);
+ -- d;
+ return tmp;
+ }
+ BOOST_UBLAS_INLINE
+ friend derived_iterator_type operator -- (derived_iterator_type &d, int) {
+ derived_iterator_type tmp (d);
+ -- d;
+ return tmp;
+ }
+ BOOST_UBLAS_INLINE
+ derived_iterator_type operator + (derived_difference_type n) const {
+ derived_iterator_type tmp (*static_cast<const derived_iterator_type *> (this));
+ return tmp += n;
+ }
+ BOOST_UBLAS_INLINE
+ friend derived_iterator_type operator + (const derived_iterator_type &d, derived_difference_type n) {
+ derived_iterator_type tmp (d);
+ return tmp += n;
+ }
+ BOOST_UBLAS_INLINE
+ friend derived_iterator_type operator + (derived_difference_type n, const derived_iterator_type &d) {
+ derived_iterator_type tmp (d);
+ return tmp += n;
+ }
+ BOOST_UBLAS_INLINE
+ derived_iterator_type operator - (derived_difference_type n) const {
+ derived_iterator_type tmp (*static_cast<const derived_iterator_type *> (this));
+ return tmp -= n;
+ }
+ BOOST_UBLAS_INLINE
+ friend derived_iterator_type operator - (const derived_iterator_type &d, derived_difference_type n) {
+ derived_iterator_type tmp (d);
+ return tmp -= n;
+ }
+
+ // Comparison
+ BOOST_UBLAS_INLINE
+ bool operator != (const derived_iterator_type &it) const {
+ const derived_iterator_type *d = static_cast<const derived_iterator_type *> (this);
+ return ! (*d == it);
+ }
+ BOOST_UBLAS_INLINE
+ bool operator <= (const derived_iterator_type &it) const {
+ const derived_iterator_type *d = static_cast<const derived_iterator_type *> (this);
+ return ! (it < *d);
+ }
+ BOOST_UBLAS_INLINE
+ bool operator >= (const derived_iterator_type &it) const {
+ const derived_iterator_type *d = static_cast<const derived_iterator_type *> (this);
+ return ! (*d < it);
+ }
+ BOOST_UBLAS_INLINE
+ bool operator > (const derived_iterator_type &it) const {
+ const derived_iterator_type *d = static_cast<const derived_iterator_type *> (this);
+ return it < *d;
+ }
+ };
+
+ /** \brief Base class of all reverse iterators. (non-MSVC version)
+ *
+ * \param I the derived iterator type
+ * \param T the value type
+ * \param R the reference type
+ *
+ * The reverse iterator implements a bidirectional iterator
+ * reversing the elements of the underlying iterator. It
+ * implements most operators of a random access iterator.
+ *
+ * uBLAS extension: it.index()
+ */
+
+ // Renamed this class from reverse_iterator to get
+ // typedef reverse_iterator<...> reverse_iterator
+ // working. Thanks to Gabriel Dos Reis for explaining this.
+ template <class I>
+ class reverse_iterator_base:
+ public std::reverse_iterator<I> {
+ public:
+ typedef typename I::container_type container_type;
+ typedef typename container_type::size_type size_type;
+ typedef typename I::difference_type difference_type;
+ typedef I iterator_type;
+
+ // Construction and destruction
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base ():
+ std::reverse_iterator<iterator_type> () {}
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base (const iterator_type &it):
+ std::reverse_iterator<iterator_type> (it) {}
+
+ // Arithmetic
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base &operator ++ () {
+ return *this = -- this->base ();
+ }
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base operator ++ (int) {
+ reverse_iterator_base tmp (*this);
+ *this = -- this->base ();
+ return tmp;
+ }
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base &operator -- () {
+ return *this = ++ this->base ();
+ }
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base operator -- (int) {
+ reverse_iterator_base tmp (*this);
+ *this = ++ this->base ();
+ return tmp;
+ }
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base &operator += (difference_type n) {
+ return *this = this->base () - n;
+ }
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base &operator -= (difference_type n) {
+ return *this = this->base () + n;
+ }
+
+ BOOST_UBLAS_INLINE
+ friend reverse_iterator_base operator + (const reverse_iterator_base &it, difference_type n) {
+ reverse_iterator_base tmp (it);
+ return tmp += n;
+ }
+ BOOST_UBLAS_INLINE
+ friend reverse_iterator_base operator + (difference_type n, const reverse_iterator_base &it) {
+ reverse_iterator_base tmp (it);
+ return tmp += n;
+ }
+ BOOST_UBLAS_INLINE
+ friend reverse_iterator_base operator - (const reverse_iterator_base &it, difference_type n) {
+ reverse_iterator_base tmp (it);
+ return tmp -= n;
+ }
+ BOOST_UBLAS_INLINE
+ friend difference_type operator - (const reverse_iterator_base &it1, const reverse_iterator_base &it2) {
+ return it2.base () - it1.base ();
+ }
+
+ BOOST_UBLAS_INLINE
+ const container_type &operator () () const {
+ return this->base () ();
+ }
+
+ BOOST_UBLAS_INLINE
+ size_type index () const {
+ iterator_type tmp (this->base ());
+ return (-- tmp).index ();
+ }
+ };
+
+ /** \brief 1st base class of all matrix reverse iterators. (non-MSVC version)
+ *
+ * \param I the derived iterator type
+ *
+ * The reverse iterator implements a bidirectional iterator
+ * reversing the elements of the underlying iterator. It
+ * implements most operators of a random access iterator.
+ *
+ * uBLAS extension: it.index1(), it.index2() and access to
+ * the dual iterator via begin(), end(), rbegin(), rend()
+ */
+
+ // Renamed this class from reverse_iterator1 to get
+ // typedef reverse_iterator1<...> reverse_iterator1
+ // working. Thanks to Gabriel Dos Reis for explaining this.
+ template <class I>
+ class reverse_iterator_base1:
+ public std::reverse_iterator<I> {
+ public:
+ typedef typename I::container_type container_type;
+ typedef typename container_type::size_type size_type;
+ typedef typename I::difference_type difference_type;
+ typedef I iterator_type;
+ typedef typename I::dual_iterator_type dual_iterator_type;
+ typedef typename I::dual_reverse_iterator_type dual_reverse_iterator_type;
+
+ // Construction and destruction
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base1 ():
+ std::reverse_iterator<iterator_type> () {}
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base1 (const iterator_type &it):
+ std::reverse_iterator<iterator_type> (it) {}
+
+ // Arithmetic
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base1 &operator ++ () {
+ return *this = -- this->base ();
+ }
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base1 operator ++ (int) {
+ reverse_iterator_base1 tmp (*this);
+ *this = -- this->base ();
+ return tmp;
+ }
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base1 &operator -- () {
+ return *this = ++ this->base ();
+ }
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base1 operator -- (int) {
+ reverse_iterator_base1 tmp (*this);
+ *this = ++ this->base ();
+ return tmp;
+ }
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base1 &operator += (difference_type n) {
+ return *this = this->base () - n;
+ }
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base1 &operator -= (difference_type n) {
+ return *this = this->base () + n;
+ }
+
+ BOOST_UBLAS_INLINE
+ friend reverse_iterator_base1 operator + (const reverse_iterator_base1 &it, difference_type n) {
+ reverse_iterator_base1 tmp (it);
+ return tmp += n;
+ }
+ BOOST_UBLAS_INLINE
+ friend reverse_iterator_base1 operator + (difference_type n, const reverse_iterator_base1 &it) {
+ reverse_iterator_base1 tmp (it);
+ return tmp += n;
+ }
+ BOOST_UBLAS_INLINE
+ friend reverse_iterator_base1 operator - (const reverse_iterator_base1 &it, difference_type n) {
+ reverse_iterator_base1 tmp (it);
+ return tmp -= n;
+ }
+ BOOST_UBLAS_INLINE
+ friend difference_type operator - (const reverse_iterator_base1 &it1, const reverse_iterator_base1 &it2) {
+ return it2.base () - it1.base ();
+ }
+
+ BOOST_UBLAS_INLINE
+ const container_type &operator () () const {
+ return this->base () ();
+ }
+
+ BOOST_UBLAS_INLINE
+ size_type index1 () const {
+ iterator_type tmp (this->base ());
+ return (-- tmp).index1 ();
+ }
+ BOOST_UBLAS_INLINE
+ size_type index2 () const {
+ iterator_type tmp (this->base ());
+ return (-- tmp).index2 ();
+ }
+
+ BOOST_UBLAS_INLINE
+ dual_iterator_type begin () const {
+ iterator_type tmp (this->base ());
+ return (-- tmp).begin ();
+ }
+ BOOST_UBLAS_INLINE
+ dual_iterator_type end () const {
+ iterator_type tmp (this->base ());
+ return (-- tmp).end ();
+ }
+ BOOST_UBLAS_INLINE
+ dual_reverse_iterator_type rbegin () const {
+ return dual_reverse_iterator_type (end ());
+ }
+ BOOST_UBLAS_INLINE
+ dual_reverse_iterator_type rend () const {
+ return dual_reverse_iterator_type (begin ());
+ }
+ };
+
+ /** \brief 2nd base class of all matrix reverse iterators. (non-MSVC version)
+ *
+ * \param I the derived iterator type
+ *
+ * The reverse iterator implements a bidirectional iterator
+ * reversing the elements of the underlying iterator. It
+ * implements most operators of a random access iterator.
+ *
+ * uBLAS extension: it.index1(), it.index2() and access to
+ * the dual iterator via begin(), end(), rbegin(), rend()
+ *
+ * Note: this type is _identical_ to reverse_iterator_base1
+ */
+
+ // Renamed this class from reverse_iterator2 to get
+ // typedef reverse_iterator2<...> reverse_iterator2
+ // working. Thanks to Gabriel Dos Reis for explaining this.
+ template <class I>
+ class reverse_iterator_base2:
+ public std::reverse_iterator<I> {
+ public:
+ typedef typename I::container_type container_type;
+ typedef typename container_type::size_type size_type;
+ typedef typename I::difference_type difference_type;
+ typedef I iterator_type;
+ typedef typename I::dual_iterator_type dual_iterator_type;
+ typedef typename I::dual_reverse_iterator_type dual_reverse_iterator_type;
+
+ // Construction and destruction
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base2 ():
+ std::reverse_iterator<iterator_type> () {}
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base2 (const iterator_type &it):
+ std::reverse_iterator<iterator_type> (it) {}
+
+ // Arithmetic
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base2 &operator ++ () {
+ return *this = -- this->base ();
+ }
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base2 operator ++ (int) {
+ reverse_iterator_base2 tmp (*this);
+ *this = -- this->base ();
+ return tmp;
+ }
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base2 &operator -- () {
+ return *this = ++ this->base ();
+ }
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base2 operator -- (int) {
+ reverse_iterator_base2 tmp (*this);
+ *this = ++ this->base ();
+ return tmp;
+ }
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base2 &operator += (difference_type n) {
+ return *this = this->base () - n;
+ }
+ BOOST_UBLAS_INLINE
+ reverse_iterator_base2 &operator -= (difference_type n) {
+ return *this = this->base () + n;
+ }
+
+ BOOST_UBLAS_INLINE
+ friend reverse_iterator_base2 operator + (const reverse_iterator_base2 &it, difference_type n) {
+ reverse_iterator_base2 tmp (it);
+ return tmp += n;
+ }
+ BOOST_UBLAS_INLINE
+ friend reverse_iterator_base2 operator + (difference_type n, const reverse_iterator_base2 &it) {
+ reverse_iterator_base2 tmp (it);
+ return tmp += n;
+ }
+ BOOST_UBLAS_INLINE
+ friend reverse_iterator_base2 operator - (const reverse_iterator_base2 &it, difference_type n) {
+ reverse_iterator_base2 tmp (it);
+ return tmp -= n;
+ }
+ BOOST_UBLAS_INLINE
+ friend difference_type operator - (const reverse_iterator_base2 &it1, const reverse_iterator_base2 &it2) {
+ return it2.base () - it1.base ();
+ }
+
+ BOOST_UBLAS_INLINE
+ const container_type &operator () () const {
+ return this->base () ();
+ }
+
+ BOOST_UBLAS_INLINE
+ size_type index1 () const {
+ iterator_type tmp (this->base ());
+ return (-- tmp).index1 ();
+ }
+ BOOST_UBLAS_INLINE
+ size_type index2 () const {
+ iterator_type tmp (this->base ());
+ return (-- tmp).index2 ();
+ }
+
+ BOOST_UBLAS_INLINE
+ dual_iterator_type begin () const {
+ iterator_type tmp (this->base ());
+ return (-- tmp).begin ();
+ }
+ BOOST_UBLAS_INLINE
+ dual_iterator_type end () const {
+ iterator_type tmp (this->base ());
+ return (-- tmp).end ();
+ }
+ BOOST_UBLAS_INLINE
+ dual_reverse_iterator_type rbegin () const {
+ return dual_reverse_iterator_type (end ());
+ }
+ BOOST_UBLAS_INLINE
+ dual_reverse_iterator_type rend () const {
+ return dual_reverse_iterator_type (begin ());
+ }
+ };
+
+ /** \brief A class implementing an indexed random access iterator.
+ *
+ * \param C the (mutable) container type
+ * \param IC the iterator category
+ *
+ * This class implements a random access iterator. The current
+ * position is stored as the unsigned integer it_ and the
+ * values are accessed via operator()(it_) of the container.
+ *
+ * uBLAS extension: index()
+ */
+
+ template<class C, class IC>
+ class indexed_iterator:
+ public container_reference<C>,
+ public random_access_iterator_base<IC,
+ indexed_iterator<C, IC>,
+ typename C::value_type,
+ typename C::difference_type> {
+ public:
+ typedef C container_type;
+ typedef IC iterator_category;
+ typedef typename container_type::size_type size_type;
+ typedef typename container_type::difference_type difference_type;
+ typedef typename container_type::value_type value_type;
+ typedef typename container_type::reference reference;
+
+ // Construction and destruction
+ BOOST_UBLAS_INLINE
+ indexed_iterator ():
+ container_reference<container_type> (), it_ () {}
+ BOOST_UBLAS_INLINE
+ indexed_iterator (container_type &c, size_type it):
+ container_reference<container_type> (c), it_ (it) {}
+
+ // Arithmetic
+ BOOST_UBLAS_INLINE
+ indexed_iterator &operator ++ () {
+ ++ it_;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ indexed_iterator &operator -- () {
+ -- it_;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ indexed_iterator &operator += (difference_type n) {
+ it_ += n;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ indexed_iterator &operator -= (difference_type n) {
+ it_ -= n;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ difference_type operator - (const indexed_iterator &it) const {
+ BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
+ return it_ - it.it_;
+ }
+
+ // Dereference
+ BOOST_UBLAS_INLINE
+ reference operator * () const {
+ BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
+ return (*this) () (it_);
+ }
+ BOOST_UBLAS_INLINE
+ reference operator [] (difference_type n) const {
+ return *((*this) + n);
+ }
+
+ // Index
+ BOOST_UBLAS_INLINE
+ size_type index () const {
+ return it_;
+ }
+
+ // Assignment
+ BOOST_UBLAS_INLINE
+ indexed_iterator &operator = (const indexed_iterator &it) {
+ // FIX: ICC needs full qualification?!
+ // assign (&it ());
+ container_reference<C>::assign (&it ());
+ it_ = it.it_;
+ return *this;
+ }
+
+ // Comparison
+ BOOST_UBLAS_INLINE
+ bool operator == (const indexed_iterator &it) const {
+ BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
+ return it_ == it.it_;
+ }
+ BOOST_UBLAS_INLINE
+ bool operator < (const indexed_iterator &it) const {
+ BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
+ return it_ < it.it_;
+ }
+
+ private:
+ size_type it_;
+ };
+
+ /** \brief A class implementing an indexed random access iterator.
+ *
+ * \param C the (immutable) container type
+ * \param IC the iterator category
+ *
+ * This class implements a random access iterator. The current
+ * position is stored as the unsigned integer \c it_ and the
+ * values are accessed via \c operator()(it_) of the container.
+ *
+ * uBLAS extension: \c index()
+ *
+ * Note: there is an automatic conversion from
+ * \c indexed_iterator to \c indexed_const_iterator
+ */
+
+ template<class C, class IC>
+ class indexed_const_iterator:
+ public container_const_reference<C>,
+ public random_access_iterator_base<IC,
+ indexed_const_iterator<C, IC>,
+ typename C::value_type,
+ typename C::difference_type> {
+ public:
+ typedef C container_type;
+ typedef IC iterator_category;
+ typedef typename container_type::size_type size_type;
+ typedef typename container_type::difference_type difference_type;
+ typedef typename container_type::value_type value_type;
+ typedef typename container_type::const_reference reference;
+ typedef indexed_iterator<container_type, iterator_category> iterator_type;
+
+ // Construction and destruction
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator ():
+ container_const_reference<container_type> (), it_ () {}
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator (const container_type &c, size_type it):
+ container_const_reference<container_type> (c), it_ (it) {}
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator (const iterator_type &it):
+ container_const_reference<container_type> (it ()), it_ (it.index ()) {}
+
+ // Arithmetic
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator &operator ++ () {
+ ++ it_;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator &operator -- () {
+ -- it_;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator &operator += (difference_type n) {
+ it_ += n;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator &operator -= (difference_type n) {
+ it_ -= n;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ difference_type operator - (const indexed_const_iterator &it) const {
+ BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
+ return it_ - it.it_;
+ }
+
+ // Dereference
+ BOOST_UBLAS_INLINE
+ reference operator * () const {
+ BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
+ return (*this) () (it_);
+ }
+ BOOST_UBLAS_INLINE
+ reference operator [] (difference_type n) const {
+ return *((*this) + n);
+ }
+
+ // Index
+ BOOST_UBLAS_INLINE
+ size_type index () const {
+ return it_;
+ }
+
+ // Assignment
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator &operator = (const indexed_const_iterator &it) {
+ // FIX: ICC needs full qualification?!
+ // assign (&it ());
+ container_const_reference<C>::assign (&it ());
+ it_ = it.it_;
+ return *this;
+ }
+
+ // Comparison
+ BOOST_UBLAS_INLINE
+ bool operator == (const indexed_const_iterator &it) const {
+ BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
+ return it_ == it.it_;
+ }
+ BOOST_UBLAS_INLINE
+ bool operator < (const indexed_const_iterator &it) const {
+ BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
+ return it_ < it.it_;
+ }
+
+ private:
+ size_type it_;
+
+ friend class indexed_iterator<container_type, iterator_category>;
+ };
+
+ template<class C, class IC>
+ class indexed_iterator2;
+
+ /** \brief A class implementing an indexed random access iterator
+ * of a matrix.
+ *
+ * \param C the (mutable) container type
+ * \param IC the iterator category
+ *
+ * This class implements a random access iterator. The current
+ * position is stored as two unsigned integers \c it1_ and \c it2_
+ * and the values are accessed via \c operator()(it1_, it2_) of the
+ * container. The iterator changes the first index.
+ *
+ * uBLAS extension: \c index1(), \c index2() and access to the
+ * dual iterator via \c begin(), \c end(), \c rbegin() and \c rend()
+ *
+ * Note: The container has to support the \code find2(rank, i, j) \endcode
+ * method
+ */
+
+ template<class C, class IC>
+ class indexed_iterator1:
+ public container_reference<C>,
+ public random_access_iterator_base<IC,
+ indexed_iterator1<C, IC>,
+ typename C::value_type,
+ typename C::difference_type> {
+ public:
+ typedef C container_type;
+ typedef IC iterator_category;
+ typedef typename container_type::size_type size_type;
+ typedef typename container_type::difference_type difference_type;
+ typedef typename container_type::value_type value_type;
+ typedef typename container_type::reference reference;
+
+ typedef indexed_iterator2<container_type, iterator_category> dual_iterator_type;
+ typedef reverse_iterator_base2<dual_iterator_type> dual_reverse_iterator_type;
+
+ // Construction and destruction
+ BOOST_UBLAS_INLINE
+ indexed_iterator1 ():
+ container_reference<container_type> (), it1_ (), it2_ () {}
+ BOOST_UBLAS_INLINE
+ indexed_iterator1 (container_type &c, size_type it1, size_type it2):
+ container_reference<container_type> (c), it1_ (it1), it2_ (it2) {}
+
+ // Arithmetic
+ BOOST_UBLAS_INLINE
+ indexed_iterator1 &operator ++ () {
+ ++ it1_;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ indexed_iterator1 &operator -- () {
+ -- it1_;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ indexed_iterator1 &operator += (difference_type n) {
+ it1_ += n;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ indexed_iterator1 &operator -= (difference_type n) {
+ it1_ -= n;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ difference_type operator - (const indexed_iterator1 &it) const {
+ BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
+ BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
+ return it1_ - it.it1_;
+ }
+
+ // Dereference
+ BOOST_UBLAS_INLINE
+ reference operator * () const {
+ BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
+ BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
+ return (*this) () (it1_, it2_);
+ }
+ BOOST_UBLAS_INLINE
+ reference operator [] (difference_type n) const {
+ return *((*this) + n);
+ }
+
+ // Index
+ BOOST_UBLAS_INLINE
+ size_type index1 () const {
+ return it1_;
+ }
+ BOOST_UBLAS_INLINE
+ size_type index2 () const {
+ return it2_;
+ }
+
+ BOOST_UBLAS_INLINE
+ dual_iterator_type begin () const {
+ return (*this) ().find2 (1, index1 (), 0);
+ }
+ BOOST_UBLAS_INLINE
+ dual_iterator_type end () const {
+ return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
+ }
+ BOOST_UBLAS_INLINE
+ dual_reverse_iterator_type rbegin () const {
+ return dual_reverse_iterator_type (end ());
+ }
+ BOOST_UBLAS_INLINE
+ dual_reverse_iterator_type rend () const {
+ return dual_reverse_iterator_type (begin ());
+ }
+
+ // Assignment
+ BOOST_UBLAS_INLINE
+ indexed_iterator1 &operator = (const indexed_iterator1 &it) {
+ // FIX: ICC needs full qualification?!
+ // assign (&it ());
+ container_reference<C>::assign (&it ());
+ it1_ = it.it1_;
+ it2_ = it.it2_;
+ return *this;
+ }
+
+ // Comparison
+ BOOST_UBLAS_INLINE
+ bool operator == (const indexed_iterator1 &it) const {
+ BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
+ BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
+ return it1_ == it.it1_;
+ }
+ BOOST_UBLAS_INLINE
+ bool operator < (const indexed_iterator1 &it) const {
+ BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
+ BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
+ return it1_ < it.it1_;
+ }
+
+ private:
+ size_type it1_;
+ size_type it2_;
+ };
+
+ template<class C, class IC>
+ class indexed_const_iterator2;
+
+ /** \brief A class implementing an indexed random access iterator
+ * of a matrix.
+ *
+ * \param C the (immutable) container type
+ * \param IC the iterator category
+ *
+ * This class implements a random access iterator. The current
+ * position is stored as two unsigned integers \c it1_ and \c it2_
+ * and the values are accessed via \c operator()(it1_, it2_) of the
+ * container. The iterator changes the first index.
+ *
+ * uBLAS extension: \c index1(), \c index2() and access to the
+ * dual iterator via \c begin(), \c end(), \c rbegin() and \c rend()
+ *
+ * Note 1: The container has to support the find2(rank, i, j) method
+ *
+ * Note 2: there is an automatic conversion from
+ * \c indexed_iterator1 to \c indexed_const_iterator1
+ */
+
+ template<class C, class IC>
+ class indexed_const_iterator1:
+ public container_const_reference<C>,
+ public random_access_iterator_base<IC,
+ indexed_const_iterator1<C, IC>,
+ typename C::value_type,
+ typename C::difference_type> {
+ public:
+ typedef C container_type;
+ typedef IC iterator_category;
+ typedef typename container_type::size_type size_type;
+ typedef typename container_type::difference_type difference_type;
+ typedef typename container_type::value_type value_type;
+ typedef typename container_type::const_reference reference;
+
+ typedef indexed_iterator1<container_type, iterator_category> iterator_type;
+ typedef indexed_const_iterator2<container_type, iterator_category> dual_iterator_type;
+ typedef reverse_iterator_base2<dual_iterator_type> dual_reverse_iterator_type;
+
+ // Construction and destruction
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator1 ():
+ container_const_reference<container_type> (), it1_ (), it2_ () {}
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator1 (const container_type &c, size_type it1, size_type it2):
+ container_const_reference<container_type> (c), it1_ (it1), it2_ (it2) {}
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator1 (const iterator_type &it):
+ container_const_reference<container_type> (it ()), it1_ (it.index1 ()), it2_ (it.index2 ()) {}
+
+ // Arithmetic
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator1 &operator ++ () {
+ ++ it1_;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator1 &operator -- () {
+ -- it1_;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator1 &operator += (difference_type n) {
+ it1_ += n;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator1 &operator -= (difference_type n) {
+ it1_ -= n;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ difference_type operator - (const indexed_const_iterator1 &it) const {
+ BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
+ BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
+ return it1_ - it.it1_;
+ }
+
+ // Dereference
+ BOOST_UBLAS_INLINE
+ reference operator * () const {
+ BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
+ BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
+ return (*this) () (it1_, it2_);
+ }
+ BOOST_UBLAS_INLINE
+ reference operator [] (difference_type n) const {
+ return *((*this) + n);
+ }
+
+ // Index
+ BOOST_UBLAS_INLINE
+ size_type index1 () const {
+ return it1_;
+ }
+ BOOST_UBLAS_INLINE
+ size_type index2 () const {
+ return it2_;
+ }
+
+ BOOST_UBLAS_INLINE
+ dual_iterator_type begin () const {
+ return (*this) ().find2 (1, index1 (), 0);
+ }
+ BOOST_UBLAS_INLINE
+ dual_iterator_type end () const {
+ return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
+ }
+ BOOST_UBLAS_INLINE
+ dual_reverse_iterator_type rbegin () const {
+ return dual_reverse_iterator_type (end ());
+ }
+ BOOST_UBLAS_INLINE
+ dual_reverse_iterator_type rend () const {
+ return dual_reverse_iterator_type (begin ());
+ }
+
+ // Assignment
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator1 &operator = (const indexed_const_iterator1 &it) {
+ // FIX: ICC needs full qualification?!
+ // assign (&it ());
+ container_const_reference<C>::assign (&it ());
+ it1_ = it.it1_;
+ it2_ = it.it2_;
+ return *this;
+ }
+
+ // Comparison
+ BOOST_UBLAS_INLINE
+ bool operator == (const indexed_const_iterator1 &it) const {
+ BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
+ BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
+ return it1_ == it.it1_;
+ }
+ BOOST_UBLAS_INLINE
+ bool operator < (const indexed_const_iterator1 &it) const {
+ BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
+ BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
+ return it1_ < it.it1_;
+ }
+
+ private:
+ size_type it1_;
+ size_type it2_;
+
+ friend class indexed_iterator1<container_type, iterator_category>;
+ };
+
+ /** \brief A class implementing an indexed random access iterator
+ * of a matrix.
+ *
+ * \param C the (mutable) container type
+ * \param IC the iterator category
+ *
+ * This class implements a random access iterator. The current
+ * position is stored as two unsigned integers \c it1_ and \c it2_
+ * and the values are accessed via \c operator()(it1_, it2_) of the
+ * container. The iterator changes the second index.
+ *
+ * uBLAS extension: \c index1(), \c index2() and access to the
+ * dual iterator via \c begin(), \c end(), \c rbegin() and \c rend()
+ *
+ * Note: The container has to support the find1(rank, i, j) method
+ */
+ template<class C, class IC>
+ class indexed_iterator2:
+ public container_reference<C>,
+ public random_access_iterator_base<IC,
+ indexed_iterator2<C, IC>,
+ typename C::value_type,
+ typename C::difference_type> {
+ public:
+ typedef C container_type;
+ typedef IC iterator_category;
+ typedef typename container_type::size_type size_type;
+ typedef typename container_type::difference_type difference_type;
+ typedef typename container_type::value_type value_type;
+ typedef typename container_type::reference reference;
+
+ typedef indexed_iterator1<container_type, iterator_category> dual_iterator_type;
+ typedef reverse_iterator_base1<dual_iterator_type> dual_reverse_iterator_type;
+
+ // Construction and destruction
+ BOOST_UBLAS_INLINE
+ indexed_iterator2 ():
+ container_reference<container_type> (), it1_ (), it2_ () {}
+ BOOST_UBLAS_INLINE
+ indexed_iterator2 (container_type &c, size_type it1, size_type it2):
+ container_reference<container_type> (c), it1_ (it1), it2_ (it2) {}
+
+ // Arithmetic
+ BOOST_UBLAS_INLINE
+ indexed_iterator2 &operator ++ () {
+ ++ it2_;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ indexed_iterator2 &operator -- () {
+ -- it2_;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ indexed_iterator2 &operator += (difference_type n) {
+ it2_ += n;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ indexed_iterator2 &operator -= (difference_type n) {
+ it2_ -= n;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ difference_type operator - (const indexed_iterator2 &it) const {
+ BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
+ BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
+ return it2_ - it.it2_;
+ }
+
+ // Dereference
+ BOOST_UBLAS_INLINE
+ reference operator * () const {
+ BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
+ BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
+ return (*this) () (it1_, it2_);
+ }
+ BOOST_UBLAS_INLINE
+ reference operator [] (difference_type n) const {
+ return *((*this) + n);
+ }
+
+ // Index
+ BOOST_UBLAS_INLINE
+ size_type index1 () const {
+ return it1_;
+ }
+ BOOST_UBLAS_INLINE
+ size_type index2 () const {
+ return it2_;
+ }
+
+ BOOST_UBLAS_INLINE
+ dual_iterator_type begin () const {
+ return (*this) ().find1 (1, 0, index2 ());
+ }
+ BOOST_UBLAS_INLINE
+ dual_iterator_type end () const {
+ return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
+ }
+ BOOST_UBLAS_INLINE
+ dual_reverse_iterator_type rbegin () const {
+ return dual_reverse_iterator_type (end ());
+ }
+ BOOST_UBLAS_INLINE
+ dual_reverse_iterator_type rend () const {
+ return dual_reverse_iterator_type (begin ());
+ }
+
+ // Assignment
+ BOOST_UBLAS_INLINE
+ indexed_iterator2 &operator = (const indexed_iterator2 &it) {
+ // FIX: ICC needs full qualification?!
+ // assign (&it ());
+ container_reference<C>::assign (&it ());
+ it1_ = it.it1_;
+ it2_ = it.it2_;
+ return *this;
+ }
+
+ // Comparison
+ BOOST_UBLAS_INLINE
+ bool operator == (const indexed_iterator2 &it) const {
+ BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
+ BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
+ return it2_ == it.it2_;
+ }
+ BOOST_UBLAS_INLINE
+ bool operator < (const indexed_iterator2 &it) const {
+ BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
+ BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
+ return it2_ < it.it2_;
+ }
+
+ private:
+ size_type it1_;
+ size_type it2_;
+ };
+
+ /** \brief A class implementing an indexed random access iterator
+ * of a matrix.
+ *
+ * \param C the (immutable) container type
+ * \param IC the iterator category
+ *
+ * This class implements a random access iterator. The current
+ * position is stored as two unsigned integers \c it1_ and \c it2_
+ * and the values are accessed via \c operator()(it1_, it2_) of the
+ * container. The iterator changes the second index.
+ *
+ * uBLAS extension: \c index1(), \c index2() and access to the
+ * dual iterator via \c begin(), \c end(), \c rbegin() and \c rend()
+ *
+ * Note 1: The container has to support the \c find2(rank, i, j) method
+ *
+ * Note 2: there is an automatic conversion from
+ * \c indexed_iterator2 to \c indexed_const_iterator2
+ */
+
+ template<class C, class IC>
+ class indexed_const_iterator2:
+ public container_const_reference<C>,
+ public random_access_iterator_base<IC,
+ indexed_const_iterator2<C, IC>,
+ typename C::value_type,
+ typename C::difference_type> {
+ public:
+ typedef C container_type;
+ typedef IC iterator_category;
+ typedef typename container_type::size_type size_type;
+ typedef typename container_type::difference_type difference_type;
+ typedef typename container_type::value_type value_type;
+ typedef typename container_type::const_reference reference;
+
+ typedef indexed_iterator2<container_type, iterator_category> iterator_type;
+ typedef indexed_const_iterator1<container_type, iterator_category> dual_iterator_type;
+ typedef reverse_iterator_base1<dual_iterator_type> dual_reverse_iterator_type;
+
+ // Construction and destruction
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator2 ():
+ container_const_reference<container_type> (), it1_ (), it2_ () {}
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator2 (const container_type &c, size_type it1, size_type it2):
+ container_const_reference<container_type> (c), it1_ (it1), it2_ (it2) {}
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator2 (const iterator_type &it):
+ container_const_reference<container_type> (it ()), it1_ (it.index1 ()), it2_ (it.index2 ()) {}
+
+ // Arithmetic
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator2 &operator ++ () {
+ ++ it2_;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator2 &operator -- () {
+ -- it2_;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator2 &operator += (difference_type n) {
+ it2_ += n;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator2 &operator -= (difference_type n) {
+ it2_ -= n;
+ return *this;
+ }
+ BOOST_UBLAS_INLINE
+ difference_type operator - (const indexed_const_iterator2 &it) const {
+ BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
+ BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
+ return it2_ - it.it2_;
+ }
+
+ // Dereference
+ BOOST_UBLAS_INLINE
+ reference operator * () const {
+ BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
+ BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
+ return (*this) () (it1_, it2_);
+ }
+ BOOST_UBLAS_INLINE
+ reference operator [] (difference_type n) const {
+ return *((*this) + n);
+ }
+
+ // Index
+ BOOST_UBLAS_INLINE
+ size_type index1 () const {
+ return it1_;
+ }
+ BOOST_UBLAS_INLINE
+ size_type index2 () const {
+ return it2_;
+ }
+
+ BOOST_UBLAS_INLINE
+ dual_iterator_type begin () const {
+ return (*this) ().find1 (1, 0, index2 ());
+ }
+ BOOST_UBLAS_INLINE
+ dual_iterator_type end () const {
+ return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
+ }
+ BOOST_UBLAS_INLINE
+ dual_reverse_iterator_type rbegin () const {
+ return dual_reverse_iterator_type (end ());
+ }
+ BOOST_UBLAS_INLINE
+ dual_reverse_iterator_type rend () const {
+ return dual_reverse_iterator_type (begin ());
+ }
+
+ // Assignment
+ BOOST_UBLAS_INLINE
+ indexed_const_iterator2 &operator = (const indexed_const_iterator2 &it) {
+ // FIX: ICC needs full qualification?!
+ // assign (&it ());
+ container_const_reference<C>::assign (&it ());
+ it1_ = it.it1_;
+ it2_ = it.it2_;
+ return *this;
+ }
+
+ // Comparison
+ BOOST_UBLAS_INLINE
+ bool operator == (const indexed_const_iterator2 &it) const {
+ BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
+ BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
+ return it2_ == it.it2_;
+ }
+ BOOST_UBLAS_INLINE
+ bool operator < (const indexed_const_iterator2 &it) const {
+ BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
+ BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
+ return it2_ < it.it2_;
+ }
+
+ private:
+ size_type it1_;
+ size_type it2_;
+
+ friend class indexed_iterator2<container_type, iterator_category>;
+ };
+
+}}}
+
+#endif
diff --git a/include/boost/numeric/ublas/detail/matrix_assign.hpp b/include/boost/numeric/ublas/detail/matrix_assign.hpp
new file mode 100644
index 0000000..be172dd
--- /dev/null
+++ b/include/boost/numeric/ublas/detail/matrix_assign.hpp
@@ -0,0 +1,1638 @@
+//
+// Copyright (c) 2000-2002
+// Joerg Walter, Mathias Koch
+//
+// 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)
+//
+// The authors gratefully acknowledge the support of
+// GeNeSys mbH & Co. KG in producing this work.
+//
+
+#ifndef _BOOST_UBLAS_MATRIX_ASSIGN_
+#define _BOOST_UBLAS_MATRIX_ASSIGN_
+
+#include <boost/numeric/ublas/traits.hpp>
+// Required for make_conformant storage
+#include <vector>
+
+// Iterators based on ideas of Jeremy Siek
+
+namespace boost { namespace numeric { namespace ublas {
+namespace detail {
+
+ // Weak equality check - useful to compare equality two arbitary matrix expression results.
+ // Since the actual expressions are unknown, we check for and arbitary error bound
+ // on the relative error.
+ // For a linear expression the infinity norm makes sense as we do not know how the elements will be
+ // combined in the expression. False positive results are inevitable for arbirary expressions!
+ template<class E1, class E2, class S>
+ BOOST_UBLAS_INLINE
+ bool equals (const matrix_expression<E1> &e1, const matrix_expression<E2> &e2, S epsilon, S min_norm) {
+ return norm_inf (e1 - e2) < epsilon *
+ std::max<S> (std::max<S> (norm_inf (e1), norm_inf (e2)), min_norm);
+ }
+
+ template<class E1, class E2>
+ BOOST_UBLAS_INLINE
+ bool expression_type_check (const matrix_expression<E1> &e1, const matrix_expression<E2> &e2) {
+ typedef typename type_traits<typename promote_traits<typename E1::value_type,
+ typename E2::value_type>::promote_type>::real_type real_type;
+ return equals (e1, e2, BOOST_UBLAS_TYPE_CHECK_EPSILON, BOOST_UBLAS_TYPE_CHECK_MIN);
+ }
+
+
+ template<class M, class E, class R>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void make_conformant (M &m, const matrix_expression<E> &e, row_major_tag, R) {
+ BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
+ BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
+ typedef R conformant_restrict_type;
+ typedef typename M::size_type size_type;
+ typedef typename M::difference_type difference_type;
+ typedef typename M::value_type value_type;
+ // FIXME unbounded_array with push_back maybe better
+ std::vector<std::pair<size_type, size_type> > index;
+ typename M::iterator1 it1 (m.begin1 ());
+ typename M::iterator1 it1_end (m.end1 ());
+ typename E::const_iterator1 it1e (e ().begin1 ());
+ typename E::const_iterator1 it1e_end (e ().end1 ());
+ while (it1 != it1_end && it1e != it1e_end) {
+ difference_type compare = it1.index1 () - it1e.index1 ();
+ if (compare == 0) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator2 it2 (it1.begin ());
+ typename M::iterator2 it2_end (it1.end ());
+ typename E::const_iterator2 it2e (it1e.begin ());
+ typename E::const_iterator2 it2e_end (it1e.end ());
+#else
+ typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
+ typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
+ typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
+ typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
+#endif
+ if (it2 != it2_end && it2e != it2e_end) {
+ size_type it2_index = it2.index2 (), it2e_index = it2e.index2 ();
+ while (true) {
+ difference_type compare2 = it2_index - it2e_index;
+ if (compare2 == 0) {
+ ++ it2, ++ it2e;
+ if (it2 != it2_end && it2e != it2e_end) {
+ it2_index = it2.index2 ();
+ it2e_index = it2e.index2 ();
+ } else
+ break;
+ } else if (compare2 < 0) {
+ increment (it2, it2_end, - compare2);
+ if (it2 != it2_end)
+ it2_index = it2.index2 ();
+ else
+ break;
+ } else if (compare2 > 0) {
+ if (conformant_restrict_type::other (it2e.index1 (), it2e.index2 ()))
+ if (static_cast<value_type>(*it2e) != value_type/*zero*/())
+ index.push_back (std::pair<size_type, size_type> (it2e.index1 (), it2e.index2 ()));
+ ++ it2e;
+ if (it2e != it2e_end)
+ it2e_index = it2e.index2 ();
+ else
+ break;
+ }
+ }
+ }
+ while (it2e != it2e_end) {
+ if (conformant_restrict_type::other (it2e.index1 (), it2e.index2 ()))
+ if (static_cast<value_type>(*it2e) != value_type/*zero*/())
+ index.push_back (std::pair<size_type, size_type> (it2e.index1 (), it2e.index2 ()));
+ ++ it2e;
+ }
+ ++ it1, ++ it1e;
+ } else if (compare < 0) {
+ increment (it1, it1_end, - compare);
+ } else if (compare > 0) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename E::const_iterator2 it2e (it1e.begin ());
+ typename E::const_iterator2 it2e_end (it1e.end ());
+#else
+ typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
+ typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
+#endif
+ while (it2e != it2e_end) {
+ if (conformant_restrict_type::other (it2e.index1 (), it2e.index2 ()))
+ if (static_cast<value_type>(*it2e) != value_type/*zero*/())
+ index.push_back (std::pair<size_type, size_type> (it2e.index1 (), it2e.index2 ()));
+ ++ it2e;
+ }
+ ++ it1e;
+ }
+ }
+ while (it1e != it1e_end) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename E::const_iterator2 it2e (it1e.begin ());
+ typename E::const_iterator2 it2e_end (it1e.end ());
+#else
+ typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
+ typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
+#endif
+ while (it2e != it2e_end) {
+ if (conformant_restrict_type::other (it2e.index1 (), it2e.index2 ()))
+ if (static_cast<value_type>(*it2e) != value_type/*zero*/())
+ index.push_back (std::pair<size_type, size_type> (it2e.index1 (), it2e.index2 ()));
+ ++ it2e;
+ }
+ ++ it1e;
+ }
+ // ISSUE proxies require insert_element
+ for (size_type k = 0; k < index.size (); ++ k)
+ m (index [k].first, index [k].second) = value_type/*zero*/();
+ }
+ template<class M, class E, class R>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void make_conformant (M &m, const matrix_expression<E> &e, column_major_tag, R) {
+ BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
+ BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
+ typedef R conformant_restrict_type;
+ typedef typename M::size_type size_type;
+ typedef typename M::difference_type difference_type;
+ typedef typename M::value_type value_type;
+ std::vector<std::pair<size_type, size_type> > index;
+ typename M::iterator2 it2 (m.begin2 ());
+ typename M::iterator2 it2_end (m.end2 ());
+ typename E::const_iterator2 it2e (e ().begin2 ());
+ typename E::const_iterator2 it2e_end (e ().end2 ());
+ while (it2 != it2_end && it2e != it2e_end) {
+ difference_type compare = it2.index2 () - it2e.index2 ();
+ if (compare == 0) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator1 it1 (it2.begin ());
+ typename M::iterator1 it1_end (it2.end ());
+ typename E::const_iterator1 it1e (it2e.begin ());
+ typename E::const_iterator1 it1e_end (it2e.end ());
+#else
+ typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
+ typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
+ typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
+ typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
+#endif
+ if (it1 != it1_end && it1e != it1e_end) {
+ size_type it1_index = it1.index1 (), it1e_index = it1e.index1 ();
+ while (true) {
+ difference_type compare2 = it1_index - it1e_index;
+ if (compare2 == 0) {
+ ++ it1, ++ it1e;
+ if (it1 != it1_end && it1e != it1e_end) {
+ it1_index = it1.index1 ();
+ it1e_index = it1e.index1 ();
+ } else
+ break;
+ } else if (compare2 < 0) {
+ increment (it1, it1_end, - compare2);
+ if (it1 != it1_end)
+ it1_index = it1.index1 ();
+ else
+ break;
+ } else if (compare2 > 0) {
+ if (conformant_restrict_type::other (it1e.index1 (), it1e.index2 ()))
+ if (static_cast<value_type>(*it1e) != value_type/*zero*/())
+ index.push_back (std::pair<size_type, size_type> (it1e.index1 (), it1e.index2 ()));
+ ++ it1e;
+ if (it1e != it1e_end)
+ it1e_index = it1e.index1 ();
+ else
+ break;
+ }
+ }
+ }
+ while (it1e != it1e_end) {
+ if (conformant_restrict_type::other (it1e.index1 (), it1e.index2 ()))
+ if (static_cast<value_type>(*it1e) != value_type/*zero*/())
+ index.push_back (std::pair<size_type, size_type> (it1e.index1 (), it1e.index2 ()));
+ ++ it1e;
+ }
+ ++ it2, ++ it2e;
+ } else if (compare < 0) {
+ increment (it2, it2_end, - compare);
+ } else if (compare > 0) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename E::const_iterator1 it1e (it2e.begin ());
+ typename E::const_iterator1 it1e_end (it2e.end ());
+#else
+ typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
+ typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
+#endif
+ while (it1e != it1e_end) {
+ if (conformant_restrict_type::other (it1e.index1 (), it1e.index2 ()))
+ if (static_cast<value_type>(*it1e) != value_type/*zero*/())
+ index.push_back (std::pair<size_type, size_type> (it1e.index1 (), it1e.index2 ()));
+ ++ it1e;
+ }
+ ++ it2e;
+ }
+ }
+ while (it2e != it2e_end) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename E::const_iterator1 it1e (it2e.begin ());
+ typename E::const_iterator1 it1e_end (it2e.end ());
+#else
+ typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
+ typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
+#endif
+ while (it1e != it1e_end) {
+ if (conformant_restrict_type::other (it1e.index1 (), it1e.index2 ()))
+ if (static_cast<value_type>(*it1e) != value_type/*zero*/())
+ index.push_back (std::pair<size_type, size_type> (it1e.index1 (), it1e.index2 ()));
+ ++ it1e;
+ }
+ ++ it2e;
+ }
+ // ISSUE proxies require insert_element
+ for (size_type k = 0; k < index.size (); ++ k)
+ m (index [k].first, index [k].second) = value_type/*zero*/();
+ }
+
+}//namespace detail
+
+
+ // Explicitly iterating row major
+ template<template <class T1, class T2> class F, class M, class T>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void iterating_matrix_assign_scalar (M &m, const T &t, row_major_tag) {
+ typedef F<typename M::iterator2::reference, T> functor_type;
+ typedef typename M::difference_type difference_type;
+ difference_type size1 (m.size1 ());
+ difference_type size2 (m.size2 ());
+ typename M::iterator1 it1 (m.begin1 ());
+ BOOST_UBLAS_CHECK (size2 == 0 || m.end1 () - it1 == size1, bad_size ());
+ while (-- size1 >= 0) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator2 it2 (it1.begin ());
+#else
+ typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
+#endif
+ BOOST_UBLAS_CHECK (it1.end () - it2 == size2, bad_size ());
+ difference_type temp_size2 (size2);
+#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
+ while (-- temp_size2 >= 0)
+ functor_type::apply (*it2, t), ++ it2;
+#else
+ DD (temp_size2, 4, r, (functor_type::apply (*it2, t), ++ it2));
+#endif
+ ++ it1;
+ }
+ }
+ // Explicitly iterating column major
+ template<template <class T1, class T2> class F, class M, class T>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void iterating_matrix_assign_scalar (M &m, const T &t, column_major_tag) {
+ typedef F<typename M::iterator1::reference, T> functor_type;
+ typedef typename M::difference_type difference_type;
+ difference_type size2 (m.size2 ());
+ difference_type size1 (m.size1 ());
+ typename M::iterator2 it2 (m.begin2 ());
+ BOOST_UBLAS_CHECK (size1 == 0 || m.end2 () - it2 == size2, bad_size ());
+ while (-- size2 >= 0) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator1 it1 (it2.begin ());
+#else
+ typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
+#endif
+ BOOST_UBLAS_CHECK (it2.end () - it1 == size1, bad_size ());
+ difference_type temp_size1 (size1);
+#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
+ while (-- temp_size1 >= 0)
+ functor_type::apply (*it1, t), ++ it1;
+#else
+ DD (temp_size1, 4, r, (functor_type::apply (*it1, t), ++ it1));
+#endif
+ ++ it2;
+ }
+ }
+ // Explicitly indexing row major
+ template<template <class T1, class T2> class F, class M, class T>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void indexing_matrix_assign_scalar (M &m, const T &t, row_major_tag) {
+ typedef F<typename M::reference, T> functor_type;
+ typedef typename M::size_type size_type;
+ size_type size1 (m.size1 ());
+ size_type size2 (m.size2 ());
+ for (size_type i = 0; i < size1; ++ i) {
+#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
+ for (size_type j = 0; j < size2; ++ j)
+ functor_type::apply (m (i, j), t);
+#else
+ size_type j (0);
+ DD (size2, 4, r, (functor_type::apply (m (i, j), t), ++ j));
+#endif
+ }
+ }
+ // Explicitly indexing column major
+ template<template <class T1, class T2> class F, class M, class T>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void indexing_matrix_assign_scalar (M &m, const T &t, column_major_tag) {
+ typedef F<typename M::reference, T> functor_type;
+ typedef typename M::size_type size_type;
+ size_type size2 (m.size2 ());
+ size_type size1 (m.size1 ());
+ for (size_type j = 0; j < size2; ++ j) {
+#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
+ for (size_type i = 0; i < size1; ++ i)
+ functor_type::apply (m (i, j), t);
+#else
+ size_type i (0);
+ DD (size1, 4, r, (functor_type::apply (m (i, j), t), ++ i));
+#endif
+ }
+ }
+
+ // Dense (proxy) case
+ template<template <class T1, class T2> class F, class M, class T, class C>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void matrix_assign_scalar (M &m, const T &t, dense_proxy_tag, C) {
+ typedef C orientation_category;
+#ifdef BOOST_UBLAS_USE_INDEXING
+ indexing_matrix_assign_scalar<F> (m, t, orientation_category ());
+#elif BOOST_UBLAS_USE_ITERATING
+ iterating_matrix_assign_scalar<F> (m, t, orientation_category ());
+#else
+ typedef typename M::size_type size_type;
+ size_type size1 (m.size1 ());
+ size_type size2 (m.size2 ());
+ if (size1 >= BOOST_UBLAS_ITERATOR_THRESHOLD &&
+ size2 >= BOOST_UBLAS_ITERATOR_THRESHOLD)
+ iterating_matrix_assign_scalar<F> (m, t, orientation_category ());
+ else
+ indexing_matrix_assign_scalar<F> (m, t, orientation_category ());
+#endif
+ }
+ // Packed (proxy) row major case
+ template<template <class T1, class T2> class F, class M, class T>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void matrix_assign_scalar (M &m, const T &t, packed_proxy_tag, row_major_tag) {
+ typedef F<typename M::iterator2::reference, T> functor_type;
+ typedef typename M::difference_type difference_type;
+ typename M::iterator1 it1 (m.begin1 ());
+ difference_type size1 (m.end1 () - it1);
+ while (-- size1 >= 0) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator2 it2 (it1.begin ());
+ difference_type size2 (it1.end () - it2);
+#else
+ typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
+ difference_type size2 (end (it1, iterator1_tag ()) - it2);
+#endif
+ while (-- size2 >= 0)
+ functor_type::apply (*it2, t), ++ it2;
+ ++ it1;
+ }
+ }
+ // Packed (proxy) column major case
+ template<template <class T1, class T2> class F, class M, class T>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void matrix_assign_scalar (M &m, const T &t, packed_proxy_tag, column_major_tag) {
+ typedef F<typename M::iterator1::reference, T> functor_type;
+ typedef typename M::difference_type difference_type;
+ typename M::iterator2 it2 (m.begin2 ());
+ difference_type size2 (m.end2 () - it2);
+ while (-- size2 >= 0) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator1 it1 (it2.begin ());
+ difference_type size1 (it2.end () - it1);
+#else
+ typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
+ difference_type size1 (end (it2, iterator2_tag ()) - it1);
+#endif
+ while (-- size1 >= 0)
+ functor_type::apply (*it1, t), ++ it1;
+ ++ it2;
+ }
+ }
+ // Sparse (proxy) row major case
+ template<template <class T1, class T2> class F, class M, class T>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void matrix_assign_scalar (M &m, const T &t, sparse_proxy_tag, row_major_tag) {
+ typedef F<typename M::iterator2::reference, T> functor_type;
+ typename M::iterator1 it1 (m.begin1 ());
+ typename M::iterator1 it1_end (m.end1 ());
+ while (it1 != it1_end) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator2 it2 (it1.begin ());
+ typename M::iterator2 it2_end (it1.end ());
+#else
+ typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
+ typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
+#endif
+ while (it2 != it2_end)
+ functor_type::apply (*it2, t), ++ it2;
+ ++ it1;
+ }
+ }
+ // Sparse (proxy) column major case
+ template<template <class T1, class T2> class F, class M, class T>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void matrix_assign_scalar (M &m, const T &t, sparse_proxy_tag, column_major_tag) {
+ typedef F<typename M::iterator1::reference, T> functor_type;
+ typename M::iterator2 it2 (m.begin2 ());
+ typename M::iterator2 it2_end (m.end2 ());
+ while (it2 != it2_end) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator1 it1 (it2.begin ());
+ typename M::iterator1 it1_end (it2.end ());
+#else
+ typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
+ typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
+#endif
+ while (it1 != it1_end)
+ functor_type::apply (*it1, t), ++ it1;
+ ++ it2;
+ }
+ }
+
+ // Dispatcher
+ template<template <class T1, class T2> class F, class M, class T>
+ BOOST_UBLAS_INLINE
+ void matrix_assign_scalar (M &m, const T &t) {
+ typedef typename M::storage_category storage_category;
+ typedef typename M::orientation_category orientation_category;
+ matrix_assign_scalar<F> (m, t, storage_category (), orientation_category ());
+ }
+
+ template<class SC, bool COMPUTED, class RI1, class RI2>
+ struct matrix_assign_traits {
+ typedef SC storage_category;
+ };
+
+ template<bool COMPUTED>
+ struct matrix_assign_traits<dense_tag, COMPUTED, packed_random_access_iterator_tag, packed_random_access_iterator_tag> {
+ typedef packed_tag storage_category;
+ };
+ template<>
+ struct matrix_assign_traits<dense_tag, false, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
+ typedef sparse_tag storage_category;
+ };
+ template<>
+ struct matrix_assign_traits<dense_tag, true, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
+ typedef sparse_proxy_tag storage_category;
+ };
+
+ template<bool COMPUTED>
+ struct matrix_assign_traits<dense_proxy_tag, COMPUTED, packed_random_access_iterator_tag, packed_random_access_iterator_tag> {
+ typedef packed_proxy_tag storage_category;
+ };
+ template<bool COMPUTED>
+ struct matrix_assign_traits<dense_proxy_tag, COMPUTED, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
+ typedef sparse_proxy_tag storage_category;
+ };
+
+ template<>
+ struct matrix_assign_traits<packed_tag, false, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
+ typedef sparse_tag storage_category;
+ };
+ template<>
+ struct matrix_assign_traits<packed_tag, true, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
+ typedef sparse_proxy_tag storage_category;
+ };
+
+ template<bool COMPUTED>
+ struct matrix_assign_traits<packed_proxy_tag, COMPUTED, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
+ typedef sparse_proxy_tag storage_category;
+ };
+
+ template<>
+ struct matrix_assign_traits<sparse_tag, true, dense_random_access_iterator_tag, dense_random_access_iterator_tag> {
+ typedef sparse_proxy_tag storage_category;
+ };
+ template<>
+ struct matrix_assign_traits<sparse_tag, true, packed_random_access_iterator_tag, packed_random_access_iterator_tag> {
+ typedef sparse_proxy_tag storage_category;
+ };
+ template<>
+ struct matrix_assign_traits<sparse_tag, true, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
+ typedef sparse_proxy_tag storage_category;
+ };
+
+ // Explicitly iterating row major
+ template<template <class T1, class T2> class F, class M, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void iterating_matrix_assign (M &m, const matrix_expression<E> &e, row_major_tag) {
+ typedef F<typename M::iterator2::reference, typename E::value_type> functor_type;
+ typedef typename M::difference_type difference_type;
+ difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), e ().size1 ()));
+ difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), e ().size2 ()));
+ typename M::iterator1 it1 (m.begin1 ());
+ BOOST_UBLAS_CHECK (size2 == 0 || m.end1 () - it1 == size1, bad_size ());
+ typename E::const_iterator1 it1e (e ().begin1 ());
+ BOOST_UBLAS_CHECK (size2 == 0 || e ().end1 () - it1e == size1, bad_size ());
+ while (-- size1 >= 0) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator2 it2 (it1.begin ());
+ typename E::const_iterator2 it2e (it1e.begin ());
+#else
+ typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
+ typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
+#endif
+ BOOST_UBLAS_CHECK (it1.end () - it2 == size2, bad_size ());
+ BOOST_UBLAS_CHECK (it1e.end () - it2e == size2, bad_size ());
+ difference_type temp_size2 (size2);
+#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
+ while (-- temp_size2 >= 0)
+ functor_type::apply (*it2, *it2e), ++ it2, ++ it2e;
+#else
+ DD (temp_size2, 2, r, (functor_type::apply (*it2, *it2e), ++ it2, ++ it2e));
+#endif
+ ++ it1, ++ it1e;
+ }
+ }
+ // Explicitly iterating column major
+ template<template <class T1, class T2> class F, class M, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void iterating_matrix_assign (M &m, const matrix_expression<E> &e, column_major_tag) {
+ typedef F<typename M::iterator1::reference, typename E::value_type> functor_type;
+ typedef typename M::difference_type difference_type;
+ difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), e ().size2 ()));
+ difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), e ().size1 ()));
+ typename M::iterator2 it2 (m.begin2 ());
+ BOOST_UBLAS_CHECK (size1 == 0 || m.end2 () - it2 == size2, bad_size ());
+ typename E::const_iterator2 it2e (e ().begin2 ());
+ BOOST_UBLAS_CHECK (size1 == 0 || e ().end2 () - it2e == size2, bad_size ());
+ while (-- size2 >= 0) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator1 it1 (it2.begin ());
+ typename E::const_iterator1 it1e (it2e.begin ());
+#else
+ typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
+ typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
+#endif
+ BOOST_UBLAS_CHECK (it2.end () - it1 == size1, bad_size ());
+ BOOST_UBLAS_CHECK (it2e.end () - it1e == size1, bad_size ());
+ difference_type temp_size1 (size1);
+#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
+ while (-- temp_size1 >= 0)
+ functor_type::apply (*it1, *it1e), ++ it1, ++ it1e;
+#else
+ DD (temp_size1, 2, r, (functor_type::apply (*it1, *it1e), ++ it1, ++ it1e));
+#endif
+ ++ it2, ++ it2e;
+ }
+ }
+ // Explicitly indexing row major
+ template<template <class T1, class T2> class F, class M, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void indexing_matrix_assign (M &m, const matrix_expression<E> &e, row_major_tag) {
+ typedef F<typename M::reference, typename E::value_type> functor_type;
+ typedef typename M::size_type size_type;
+ size_type size1 (BOOST_UBLAS_SAME (m.size1 (), e ().size1 ()));
+ size_type size2 (BOOST_UBLAS_SAME (m.size2 (), e ().size2 ()));
+ for (size_type i = 0; i < size1; ++ i) {
+#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
+ for (size_type j = 0; j < size2; ++ j)
+ functor_type::apply (m (i, j), e () (i, j));
+#else
+ size_type j (0);
+ DD (size2, 2, r, (functor_type::apply (m (i, j), e () (i, j)), ++ j));
+#endif
+ }
+ }
+ // Explicitly indexing column major
+ template<template <class T1, class T2> class F, class M, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void indexing_matrix_assign (M &m, const matrix_expression<E> &e, column_major_tag) {
+ typedef F<typename M::reference, typename E::value_type> functor_type;
+ typedef typename M::size_type size_type;
+ size_type size2 (BOOST_UBLAS_SAME (m.size2 (), e ().size2 ()));
+ size_type size1 (BOOST_UBLAS_SAME (m.size1 (), e ().size1 ()));
+ for (size_type j = 0; j < size2; ++ j) {
+#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
+ for (size_type i = 0; i < size1; ++ i)
+ functor_type::apply (m (i, j), e () (i, j));
+#else
+ size_type i (0);
+ DD (size1, 2, r, (functor_type::apply (m (i, j), e () (i, j)), ++ i));
+#endif
+ }
+ }
+
+ // Dense (proxy) case
+ template<template <class T1, class T2> class F, class R, class M, class E, class C>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void matrix_assign (M &m, const matrix_expression<E> &e, dense_proxy_tag, C) {
+ // R unnecessary, make_conformant not required
+ typedef C orientation_category;
+#ifdef BOOST_UBLAS_USE_INDEXING
+ indexing_matrix_assign<F> (m, e, orientation_category ());
+#elif BOOST_UBLAS_USE_ITERATING
+ iterating_matrix_assign<F> (m, e, orientation_category ());
+#else
+ typedef typename M::difference_type difference_type;
+ size_type size1 (BOOST_UBLAS_SAME (m.size1 (), e ().size1 ()));
+ size_type size2 (BOOST_UBLAS_SAME (m.size2 (), e ().size2 ()));
+ if (size1 >= BOOST_UBLAS_ITERATOR_THRESHOLD &&
+ size2 >= BOOST_UBLAS_ITERATOR_THRESHOLD)
+ iterating_matrix_assign<F> (m, e, orientation_category ());
+ else
+ indexing_matrix_assign<F> (m, e, orientation_category ());
+#endif
+ }
+ // Packed (proxy) row major case
+ template<template <class T1, class T2> class F, class R, class M, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void matrix_assign (M &m, const matrix_expression<E> &e, packed_proxy_tag, row_major_tag) {
+ typedef typename matrix_traits<E>::value_type expr_value_type;
+ typedef F<typename M::iterator2::reference, expr_value_type> functor_type;
+ // R unnecessary, make_conformant not required
+ typedef typename M::difference_type difference_type;
+
+ BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
+ BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
+
+#if BOOST_UBLAS_TYPE_CHECK
+ typedef typename M::value_type value_type;
+ matrix<value_type, row_major> cm (m.size1 (), m.size2 ());
+ indexing_matrix_assign<scalar_assign> (cm, m, row_major_tag ());
+ indexing_matrix_assign<F> (cm, e, row_major_tag ());
+#endif
+ typename M::iterator1 it1 (m.begin1 ());
+ typename M::iterator1 it1_end (m.end1 ());
+ typename E::const_iterator1 it1e (e ().begin1 ());
+ typename E::const_iterator1 it1e_end (e ().end1 ());
+ difference_type it1_size (it1_end - it1);
+ difference_type it1e_size (it1e_end - it1e);
+ difference_type diff1 (0);
+ if (it1_size > 0 && it1e_size > 0)
+ diff1 = it1.index1 () - it1e.index1 ();
+ if (diff1 != 0) {
+ difference_type size1 = (std::min) (diff1, it1e_size);
+ if (size1 > 0) {
+ it1e += size1;
+ it1e_size -= size1;
+ diff1 -= size1;
+ }
+ size1 = (std::min) (- diff1, it1_size);
+ if (size1 > 0) {
+ it1_size -= size1;
+ if (!functor_type::computed) {
+ while (-- size1 >= 0) { // zeroing
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator2 it2 (it1.begin ());
+ typename M::iterator2 it2_end (it1.end ());
+#else
+ typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
+ typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
+#endif
+ difference_type size2 (it2_end - it2);
+ while (-- size2 >= 0)
+ functor_type::apply (*it2, expr_value_type/*zero*/()), ++ it2;
+ ++ it1;
+ }
+ } else {
+ it1 += size1;
+ }
+ diff1 += size1;
+ }
+ }
+ difference_type size1 ((std::min) (it1_size, it1e_size));
+ it1_size -= size1;
+ it1e_size -= size1;
+ while (-- size1 >= 0) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator2 it2 (it1.begin ());
+ typename M::iterator2 it2_end (it1.end ());
+ typename E::const_iterator2 it2e (it1e.begin ());
+ typename E::const_iterator2 it2e_end (it1e.end ());
+#else
+ typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
+ typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
+ typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
+ typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
+#endif
+ difference_type it2_size (it2_end - it2);
+ difference_type it2e_size (it2e_end - it2e);
+ difference_type diff2 (0);
+ if (it2_size > 0 && it2e_size > 0) {
+ diff2 = it2.index2 () - it2e.index2 ();
+ difference_type size2 = (std::min) (diff2, it2e_size);
+ if (size2 > 0) {
+ it2e += size2;
+ it2e_size -= size2;
+ diff2 -= size2;
+ }
+ size2 = (std::min) (- diff2, it2_size);
+ if (size2 > 0) {
+ it2_size -= size2;
+ if (!functor_type::computed) {
+ while (-- size2 >= 0) // zeroing
+ functor_type::apply (*it2, expr_value_type/*zero*/()), ++ it2;
+ } else {
+ it2 += size2;
+ }
+ diff2 += size2;
+ }
+ }
+ difference_type size2 ((std::min) (it2_size, it2e_size));
+ it2_size -= size2;
+ it2e_size -= size2;
+ while (-- size2 >= 0)
+ functor_type::apply (*it2, *it2e), ++ it2, ++ it2e;
+ size2 = it2_size;
+ if (!functor_type::computed) {
+ while (-- size2 >= 0) // zeroing
+ functor_type::apply (*it2, expr_value_type/*zero*/()), ++ it2;
+ } else {
+ it2 += size2;
+ }
+ ++ it1, ++ it1e;
+ }
+ size1 = it1_size;
+ if (!functor_type::computed) {
+ while (-- size1 >= 0) { // zeroing
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator2 it2 (it1.begin ());
+ typename M::iterator2 it2_end (it1.end ());
+#else
+ typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
+ typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
+#endif
+ difference_type size2 (it2_end - it2);
+ while (-- size2 >= 0)
+ functor_type::apply (*it2, expr_value_type/*zero*/()), ++ it2;
+ ++ it1;
+ }
+ } else {
+ it1 += size1;
+ }
+#if BOOST_UBLAS_TYPE_CHECK
+ if (! disable_type_check<bool>::value)
+ BOOST_UBLAS_CHECK (detail::expression_type_check (m, cm), external_logic ());
+#endif
+ }
+ // Packed (proxy) column major case
+ template<template <class T1, class T2> class F, class R, class M, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void matrix_assign (M &m, const matrix_expression<E> &e, packed_proxy_tag, column_major_tag) {
+ typedef typename matrix_traits<E>::value_type expr_value_type;
+ typedef F<typename M::iterator1::reference, expr_value_type> functor_type;
+ // R unnecessary, make_conformant not required
+ typedef typename M::difference_type difference_type;
+
+ BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
+ BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
+
+#if BOOST_UBLAS_TYPE_CHECK
+ typedef typename M::value_type value_type;
+ matrix<value_type, column_major> cm (m.size1 (), m.size2 ());
+ indexing_matrix_assign<scalar_assign> (cm, m, column_major_tag ());
+ indexing_matrix_assign<F> (cm, e, column_major_tag ());
+#endif
+ typename M::iterator2 it2 (m.begin2 ());
+ typename M::iterator2 it2_end (m.end2 ());
+ typename E::const_iterator2 it2e (e ().begin2 ());
+ typename E::const_iterator2 it2e_end (e ().end2 ());
+ difference_type it2_size (it2_end - it2);
+ difference_type it2e_size (it2e_end - it2e);
+ difference_type diff2 (0);
+ if (it2_size > 0 && it2e_size > 0)
+ diff2 = it2.index2 () - it2e.index2 ();
+ if (diff2 != 0) {
+ difference_type size2 = (std::min) (diff2, it2e_size);
+ if (size2 > 0) {
+ it2e += size2;
+ it2e_size -= size2;
+ diff2 -= size2;
+ }
+ size2 = (std::min) (- diff2, it2_size);
+ if (size2 > 0) {
+ it2_size -= size2;
+ if (!functor_type::computed) {
+ while (-- size2 >= 0) { // zeroing
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator1 it1 (it2.begin ());
+ typename M::iterator1 it1_end (it2.end ());
+#else
+ typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
+ typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
+#endif
+ difference_type size1 (it1_end - it1);
+ while (-- size1 >= 0)
+ functor_type::apply (*it1, expr_value_type/*zero*/()), ++ it1;
+ ++ it2;
+ }
+ } else {
+ it2 += size2;
+ }
+ diff2 += size2;
+ }
+ }
+ difference_type size2 ((std::min) (it2_size, it2e_size));
+ it2_size -= size2;
+ it2e_size -= size2;
+ while (-- size2 >= 0) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator1 it1 (it2.begin ());
+ typename M::iterator1 it1_end (it2.end ());
+ typename E::const_iterator1 it1e (it2e.begin ());
+ typename E::const_iterator1 it1e_end (it2e.end ());
+#else
+ typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
+ typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
+ typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
+ typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
+#endif
+ difference_type it1_size (it1_end - it1);
+ difference_type it1e_size (it1e_end - it1e);
+ difference_type diff1 (0);
+ if (it1_size > 0 && it1e_size > 0) {
+ diff1 = it1.index1 () - it1e.index1 ();
+ difference_type size1 = (std::min) (diff1, it1e_size);
+ if (size1 > 0) {
+ it1e += size1;
+ it1e_size -= size1;
+ diff1 -= size1;
+ }
+ size1 = (std::min) (- diff1, it1_size);
+ if (size1 > 0) {
+ it1_size -= size1;
+ if (!functor_type::computed) {
+ while (-- size1 >= 0) // zeroing
+ functor_type::apply (*it1, expr_value_type/*zero*/()), ++ it1;
+ } else {
+ it1 += size1;
+ }
+ diff1 += size1;
+ }
+ }
+ difference_type size1 ((std::min) (it1_size, it1e_size));
+ it1_size -= size1;
+ it1e_size -= size1;
+ while (-- size1 >= 0)
+ functor_type::apply (*it1, *it1e), ++ it1, ++ it1e;
+ size1 = it1_size;
+ if (!functor_type::computed) {
+ while (-- size1 >= 0) // zeroing
+ functor_type::apply (*it1, expr_value_type/*zero*/()), ++ it1;
+ } else {
+ it1 += size1;
+ }
+ ++ it2, ++ it2e;
+ }
+ size2 = it2_size;
+ if (!functor_type::computed) {
+ while (-- size2 >= 0) { // zeroing
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator1 it1 (it2.begin ());
+ typename M::iterator1 it1_end (it2.end ());
+#else
+ typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
+ typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
+#endif
+ difference_type size1 (it1_end - it1);
+ while (-- size1 >= 0)
+ functor_type::apply (*it1, expr_value_type/*zero*/()), ++ it1;
+ ++ it2;
+ }
+ } else {
+ it2 += size2;
+ }
+#if BOOST_UBLAS_TYPE_CHECK
+ if (! disable_type_check<bool>::value)
+ BOOST_UBLAS_CHECK (detail::expression_type_check (m, cm), external_logic ());
+#endif
+ }
+ // Sparse row major case
+ template<template <class T1, class T2> class F, class R, class M, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void matrix_assign (M &m, const matrix_expression<E> &e, sparse_tag, row_major_tag) {
+ typedef F<typename M::iterator2::reference, typename E::value_type> functor_type;
+ // R unnecessary, make_conformant not required
+ BOOST_STATIC_ASSERT ((!functor_type::computed));
+ BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
+ BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
+ typedef typename M::value_type value_type;
+ // Sparse type has no numeric constraints to check
+
+ m.clear ();
+ typename E::const_iterator1 it1e (e ().begin1 ());
+ typename E::const_iterator1 it1e_end (e ().end1 ());
+ while (it1e != it1e_end) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename E::const_iterator2 it2e (it1e.begin ());
+ typename E::const_iterator2 it2e_end (it1e.end ());
+#else
+ typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
+ typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
+#endif
+ while (it2e != it2e_end) {
+ value_type t (*it2e);
+ if (t != value_type/*zero*/())
+ m.insert_element (it2e.index1 (), it2e.index2 (), t);
+ ++ it2e;
+ }
+ ++ it1e;
+ }
+ }
+ // Sparse column major case
+ template<template <class T1, class T2> class F, class R, class M, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void matrix_assign (M &m, const matrix_expression<E> &e, sparse_tag, column_major_tag) {
+ typedef F<typename M::iterator1::reference, typename E::value_type> functor_type;
+ // R unnecessary, make_conformant not required
+ BOOST_STATIC_ASSERT ((!functor_type::computed));
+ BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
+ BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
+ typedef typename M::value_type value_type;
+ // Sparse type has no numeric constraints to check
+
+ m.clear ();
+ typename E::const_iterator2 it2e (e ().begin2 ());
+ typename E::const_iterator2 it2e_end (e ().end2 ());
+ while (it2e != it2e_end) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename E::const_iterator1 it1e (it2e.begin ());
+ typename E::const_iterator1 it1e_end (it2e.end ());
+#else
+ typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
+ typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
+#endif
+ while (it1e != it1e_end) {
+ value_type t (*it1e);
+ if (t != value_type/*zero*/())
+ m.insert_element (it1e.index1 (), it1e.index2 (), t);
+ ++ it1e;
+ }
+ ++ it2e;
+ }
+ }
+ // Sparse proxy or functional row major case
+ template<template <class T1, class T2> class F, class R, class M, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void matrix_assign (M &m, const matrix_expression<E> &e, sparse_proxy_tag, row_major_tag) {
+ typedef typename matrix_traits<E>::value_type expr_value_type;
+ typedef F<typename M::iterator2::reference, expr_value_type> functor_type;
+ typedef R conformant_restrict_type;
+ typedef typename M::size_type size_type;
+ typedef typename M::difference_type difference_type;
+
+ BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
+ BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
+
+#if BOOST_UBLAS_TYPE_CHECK
+ typedef typename M::value_type value_type;
+ matrix<value_type, row_major> cm (m.size1 (), m.size2 ());
+ indexing_matrix_assign<scalar_assign> (cm, m, row_major_tag ());
+ indexing_matrix_assign<F> (cm, e, row_major_tag ());
+#endif
+ detail::make_conformant (m, e, row_major_tag (), conformant_restrict_type ());
+
+ typename M::iterator1 it1 (m.begin1 ());
+ typename M::iterator1 it1_end (m.end1 ());
+ typename E::const_iterator1 it1e (e ().begin1 ());
+ typename E::const_iterator1 it1e_end (e ().end1 ());
+ while (it1 != it1_end && it1e != it1e_end) {
+ difference_type compare = it1.index1 () - it1e.index1 ();
+ if (compare == 0) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator2 it2 (it1.begin ());
+ typename M::iterator2 it2_end (it1.end ());
+ typename E::const_iterator2 it2e (it1e.begin ());
+ typename E::const_iterator2 it2e_end (it1e.end ());
+#else
+ typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
+ typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
+ typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
+ typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
+#endif
+ if (it2 != it2_end && it2e != it2e_end) {
+ size_type it2_index = it2.index2 (), it2e_index = it2e.index2 ();
+ while (true) {
+ difference_type compare2 = it2_index - it2e_index;
+ if (compare2 == 0) {
+ functor_type::apply (*it2, *it2e);
+ ++ it2, ++ it2e;
+ if (it2 != it2_end && it2e != it2e_end) {
+ it2_index = it2.index2 ();
+ it2e_index = it2e.index2 ();
+ } else
+ break;
+ } else if (compare2 < 0) {
+ if (!functor_type::computed) {
+ functor_type::apply (*it2, expr_value_type/*zero*/());
+ ++ it2;
+ } else
+ increment (it2, it2_end, - compare2);
+ if (it2 != it2_end)
+ it2_index = it2.index2 ();
+ else
+ break;
+ } else if (compare2 > 0) {
+ increment (it2e, it2e_end, compare2);
+ if (it2e != it2e_end)
+ it2e_index = it2e.index2 ();
+ else
+ break;
+ }
+ }
+ }
+ if (!functor_type::computed) {
+ while (it2 != it2_end) { // zeroing
+ functor_type::apply (*it2, expr_value_type/*zero*/());
+ ++ it2;
+ }
+ } else {
+ it2 = it2_end;
+ }
+ ++ it1, ++ it1e;
+ } else if (compare < 0) {
+ if (!functor_type::computed) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator2 it2 (it1.begin ());
+ typename M::iterator2 it2_end (it1.end ());
+#else
+ typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
+ typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
+#endif
+ while (it2 != it2_end) { // zeroing
+ functor_type::apply (*it2, expr_value_type/*zero*/());
+ ++ it2;
+ }
+ ++ it1;
+ } else {
+ increment (it1, it1_end, - compare);
+ }
+ } else if (compare > 0) {
+ increment (it1e, it1e_end, compare);
+ }
+ }
+ if (!functor_type::computed) {
+ while (it1 != it1_end) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator2 it2 (it1.begin ());
+ typename M::iterator2 it2_end (it1.end ());
+#else
+ typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
+ typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
+#endif
+ while (it2 != it2_end) { // zeroing
+ functor_type::apply (*it2, expr_value_type/*zero*/());
+ ++ it2;
+ }
+ ++ it1;
+ }
+ } else {
+ it1 = it1_end;
+ }
+#if BOOST_UBLAS_TYPE_CHECK
+ if (! disable_type_check<bool>::value)
+ BOOST_UBLAS_CHECK (detail::expression_type_check (m, cm), external_logic ());
+#endif
+ }
+ // Sparse proxy or functional column major case
+ template<template <class T1, class T2> class F, class R, class M, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void matrix_assign (M &m, const matrix_expression<E> &e, sparse_proxy_tag, column_major_tag) {
+ typedef typename matrix_traits<E>::value_type expr_value_type;
+ typedef F<typename M::iterator1::reference, expr_value_type> functor_type;
+ typedef R conformant_restrict_type;
+ typedef typename M::size_type size_type;
+ typedef typename M::difference_type difference_type;
+
+ BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
+ BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
+
+#if BOOST_UBLAS_TYPE_CHECK
+ typedef typename M::value_type value_type;
+ matrix<value_type, column_major> cm (m.size1 (), m.size2 ());
+ indexing_matrix_assign<scalar_assign> (cm, m, column_major_tag ());
+ indexing_matrix_assign<F> (cm, e, column_major_tag ());
+#endif
+ detail::make_conformant (m, e, column_major_tag (), conformant_restrict_type ());
+
+ typename M::iterator2 it2 (m.begin2 ());
+ typename M::iterator2 it2_end (m.end2 ());
+ typename E::const_iterator2 it2e (e ().begin2 ());
+ typename E::const_iterator2 it2e_end (e ().end2 ());
+ while (it2 != it2_end && it2e != it2e_end) {
+ difference_type compare = it2.index2 () - it2e.index2 ();
+ if (compare == 0) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator1 it1 (it2.begin ());
+ typename M::iterator1 it1_end (it2.end ());
+ typename E::const_iterator1 it1e (it2e.begin ());
+ typename E::const_iterator1 it1e_end (it2e.end ());
+#else
+ typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
+ typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
+ typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
+ typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
+#endif
+ if (it1 != it1_end && it1e != it1e_end) {
+ size_type it1_index = it1.index1 (), it1e_index = it1e.index1 ();
+ while (true) {
+ difference_type compare2 = it1_index - it1e_index;
+ if (compare2 == 0) {
+ functor_type::apply (*it1, *it1e);
+ ++ it1, ++ it1e;
+ if (it1 != it1_end && it1e != it1e_end) {
+ it1_index = it1.index1 ();
+ it1e_index = it1e.index1 ();
+ } else
+ break;
+ } else if (compare2 < 0) {
+ if (!functor_type::computed) {
+ functor_type::apply (*it1, expr_value_type/*zero*/()); // zeroing
+ ++ it1;
+ } else
+ increment (it1, it1_end, - compare2);
+ if (it1 != it1_end)
+ it1_index = it1.index1 ();
+ else
+ break;
+ } else if (compare2 > 0) {
+ increment (it1e, it1e_end, compare2);
+ if (it1e != it1e_end)
+ it1e_index = it1e.index1 ();
+ else
+ break;
+ }
+ }
+ }
+ if (!functor_type::computed) {
+ while (it1 != it1_end) { // zeroing
+ functor_type::apply (*it1, expr_value_type/*zero*/());
+ ++ it1;
+ }
+ } else {
+ it1 = it1_end;
+ }
+ ++ it2, ++ it2e;
+ } else if (compare < 0) {
+ if (!functor_type::computed) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator1 it1 (it2.begin ());
+ typename M::iterator1 it1_end (it2.end ());
+#else
+ typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
+ typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
+#endif
+ while (it1 != it1_end) { // zeroing
+ functor_type::apply (*it1, expr_value_type/*zero*/());
+ ++ it1;
+ }
+ ++ it2;
+ } else {
+ increment (it2, it2_end, - compare);
+ }
+ } else if (compare > 0) {
+ increment (it2e, it2e_end, compare);
+ }
+ }
+ if (!functor_type::computed) {
+ while (it2 != it2_end) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator1 it1 (it2.begin ());
+ typename M::iterator1 it1_end (it2.end ());
+#else
+ typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
+ typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
+#endif
+ while (it1 != it1_end) { // zeroing
+ functor_type::apply (*it1, expr_value_type/*zero*/());
+ ++ it1;
+ }
+ ++ it2;
+ }
+ } else {
+ it2 = it2_end;
+ }
+#if BOOST_UBLAS_TYPE_CHECK
+ if (! disable_type_check<bool>::value)
+ BOOST_UBLAS_CHECK (detail::expression_type_check (m, cm), external_logic ());
+#endif
+ }
+
+ // Dispatcher
+ template<template <class T1, class T2> class F, class M, class E>
+ BOOST_UBLAS_INLINE
+ void matrix_assign (M &m, const matrix_expression<E> &e) {
+ typedef typename matrix_assign_traits<typename M::storage_category,
+ F<typename M::reference, typename E::value_type>::computed,
+ typename E::const_iterator1::iterator_category,
+ typename E::const_iterator2::iterator_category>::storage_category storage_category;
+ // give preference to matrix M's orientation if known
+ typedef typename boost::mpl::if_<boost::is_same<typename M::orientation_category, unknown_orientation_tag>,
+ typename E::orientation_category ,
+ typename M::orientation_category >::type orientation_category;
+ typedef basic_full<typename M::size_type> unrestricted;
+ matrix_assign<F, unrestricted> (m, e, storage_category (), orientation_category ());
+ }
+ template<template <class T1, class T2> class F, class R, class M, class E>
+ BOOST_UBLAS_INLINE
+ void matrix_assign (M &m, const matrix_expression<E> &e) {
+ typedef R conformant_restrict_type;
+ typedef typename matrix_assign_traits<typename M::storage_category,
+ F<typename M::reference, typename E::value_type>::computed,
+ typename E::const_iterator1::iterator_category,
+ typename E::const_iterator2::iterator_category>::storage_category storage_category;
+ // give preference to matrix M's orientation if known
+ typedef typename boost::mpl::if_<boost::is_same<typename M::orientation_category, unknown_orientation_tag>,
+ typename E::orientation_category ,
+ typename M::orientation_category >::type orientation_category;
+ matrix_assign<F, conformant_restrict_type> (m, e, storage_category (), orientation_category ());
+ }
+
+ template<class SC, class RI1, class RI2>
+ struct matrix_swap_traits {
+ typedef SC storage_category;
+ };
+
+ template<>
+ struct matrix_swap_traits<dense_proxy_tag, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
+ typedef sparse_proxy_tag storage_category;
+ };
+
+ template<>
+ struct matrix_swap_traits<packed_proxy_tag, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
+ typedef sparse_proxy_tag storage_category;
+ };
+
+ // Dense (proxy) row major case
+ template<template <class T1, class T2> class F, class R, class M, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void matrix_swap (M &m, matrix_expression<E> &e, dense_proxy_tag, row_major_tag) {
+ typedef F<typename M::iterator2::reference, typename E::reference> functor_type;
+ // R unnecessary, make_conformant not required
+ //typedef typename M::size_type size_type; // gcc is complaining that this is not used, although this is not right
+ typedef typename M::difference_type difference_type;
+ typename M::iterator1 it1 (m.begin1 ());
+ typename E::iterator1 it1e (e ().begin1 ());
+ difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), typename M::size_type (e ().end1 () - it1e)));
+ while (-- size1 >= 0) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator2 it2 (it1.begin ());
+ typename E::iterator2 it2e (it1e.begin ());
+ difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), typename M::size_type (it1e.end () - it2e)));
+#else
+ typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
+ typename E::iterator2 it2e (begin (it1e, iterator1_tag ()));
+ difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), typename M::size_type (end (it1e, iterator1_tag ()) - it2e)));
+#endif
+ while (-- size2 >= 0)
+ functor_type::apply (*it2, *it2e), ++ it2, ++ it2e;
+ ++ it1, ++ it1e;
+ }
+ }
+ // Dense (proxy) column major case
+ template<template <class T1, class T2> class F, class R, class M, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void matrix_swap (M &m, matrix_expression<E> &e, dense_proxy_tag, column_major_tag) {
+ typedef F<typename M::iterator1::reference, typename E::reference> functor_type;
+ // R unnecessary, make_conformant not required
+ // typedef typename M::size_type size_type; // gcc is complaining that this is not used, although this is not right
+ typedef typename M::difference_type difference_type;
+ typename M::iterator2 it2 (m.begin2 ());
+ typename E::iterator2 it2e (e ().begin2 ());
+ difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), typename M::size_type (e ().end2 () - it2e)));
+ while (-- size2 >= 0) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator1 it1 (it2.begin ());
+ typename E::iterator1 it1e (it2e.begin ());
+ difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), typename M::size_type (it2e.end () - it1e)));
+#else
+ typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
+ typename E::iterator1 it1e (begin (it2e, iterator2_tag ()));
+ difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), typename M::size_type (end (it2e, iterator2_tag ()) - it1e)));
+#endif
+ while (-- size1 >= 0)
+ functor_type::apply (*it1, *it1e), ++ it1, ++ it1e;
+ ++ it2, ++ it2e;
+ }
+ }
+ // Packed (proxy) row major case
+ template<template <class T1, class T2> class F, class R, class M, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void matrix_swap (M &m, matrix_expression<E> &e, packed_proxy_tag, row_major_tag) {
+ typedef F<typename M::iterator2::reference, typename E::reference> functor_type;
+ // R unnecessary, make_conformant not required
+ typedef typename M::difference_type difference_type;
+ typename M::iterator1 it1 (m.begin1 ());
+ typename E::iterator1 it1e (e ().begin1 ());
+ difference_type size1 (BOOST_UBLAS_SAME (m.end1 () - it1, e ().end1 () - it1e));
+ while (-- size1 >= 0) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator2 it2 (it1.begin ());
+ typename E::iterator2 it2e (it1e.begin ());
+ difference_type size2 (BOOST_UBLAS_SAME (it1.end () - it2, it1e.end () - it2e));
+#else
+ typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
+ typename E::iterator2 it2e (begin (it1e, iterator1_tag ()));
+ difference_type size2 (BOOST_UBLAS_SAME (end (it1, iterator1_tag ()) - it2, end (it1e, iterator1_tag ()) - it2e));
+#endif
+ while (-- size2 >= 0)
+ functor_type::apply (*it2, *it2e), ++ it2, ++ it2e;
+ ++ it1, ++ it1e;
+ }
+ }
+ // Packed (proxy) column major case
+ template<template <class T1, class T2> class F, class R, class M, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void matrix_swap (M &m, matrix_expression<E> &e, packed_proxy_tag, column_major_tag) {
+ typedef F<typename M::iterator1::reference, typename E::reference> functor_type;
+ // R unnecessary, make_conformant not required
+ typedef typename M::difference_type difference_type;
+ typename M::iterator2 it2 (m.begin2 ());
+ typename E::iterator2 it2e (e ().begin2 ());
+ difference_type size2 (BOOST_UBLAS_SAME (m.end2 () - it2, e ().end2 () - it2e));
+ while (-- size2 >= 0) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator1 it1 (it2.begin ());
+ typename E::iterator1 it1e (it2e.begin ());
+ difference_type size1 (BOOST_UBLAS_SAME (it2.end () - it1, it2e.end () - it1e));
+#else
+ typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
+ typename E::iterator1 it1e (begin (it2e, iterator2_tag ()));
+ difference_type size1 (BOOST_UBLAS_SAME (end (it2, iterator2_tag ()) - it1, end (it2e, iterator2_tag ()) - it1e));
+#endif
+ while (-- size1 >= 0)
+ functor_type::apply (*it1, *it1e), ++ it1, ++ it1e;
+ ++ it2, ++ it2e;
+ }
+ }
+ // Sparse (proxy) row major case
+ template<template <class T1, class T2> class F, class R, class M, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void matrix_swap (M &m, matrix_expression<E> &e, sparse_proxy_tag, row_major_tag) {
+ typedef F<typename M::iterator2::reference, typename E::reference> functor_type;
+ typedef R conformant_restrict_type;
+ typedef typename M::size_type size_type;
+ typedef typename M::difference_type difference_type;
+ BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
+ BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
+
+ detail::make_conformant (m, e, row_major_tag (), conformant_restrict_type ());
+ // FIXME should be a seperate restriction for E
+ detail::make_conformant (e (), m, row_major_tag (), conformant_restrict_type ());
+
+ typename M::iterator1 it1 (m.begin1 ());
+ typename M::iterator1 it1_end (m.end1 ());
+ typename E::iterator1 it1e (e ().begin1 ());
+ typename E::iterator1 it1e_end (e ().end1 ());
+ while (it1 != it1_end && it1e != it1e_end) {
+ difference_type compare = it1.index1 () - it1e.index1 ();
+ if (compare == 0) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator2 it2 (it1.begin ());
+ typename M::iterator2 it2_end (it1.end ());
+ typename E::iterator2 it2e (it1e.begin ());
+ typename E::iterator2 it2e_end (it1e.end ());
+#else
+ typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
+ typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
+ typename E::iterator2 it2e (begin (it1e, iterator1_tag ()));
+ typename E::iterator2 it2e_end (end (it1e, iterator1_tag ()));
+#endif
+ if (it2 != it2_end && it2e != it2e_end) {
+ size_type it2_index = it2.index2 (), it2e_index = it2e.index2 ();
+ while (true) {
+ difference_type compare2 = it2_index - it2e_index;
+ if (compare2 == 0) {
+ functor_type::apply (*it2, *it2e);
+ ++ it2, ++ it2e;
+ if (it2 != it2_end && it2e != it2e_end) {
+ it2_index = it2.index2 ();
+ it2e_index = it2e.index2 ();
+ } else
+ break;
+ } else if (compare2 < 0) {
+ increment (it2, it2_end, - compare2);
+ if (it2 != it2_end)
+ it2_index = it2.index2 ();
+ else
+ break;
+ } else if (compare2 > 0) {
+ increment (it2e, it2e_end, compare2);
+ if (it2e != it2e_end)
+ it2e_index = it2e.index2 ();
+ else
+ break;
+ }
+ }
+ }
+#if BOOST_UBLAS_TYPE_CHECK
+ increment (it2e, it2e_end);
+ increment (it2, it2_end);
+#endif
+ ++ it1, ++ it1e;
+ } else if (compare < 0) {
+#if BOOST_UBLAS_TYPE_CHECK
+ while (it1.index1 () < it1e.index1 ()) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator2 it2 (it1.begin ());
+ typename M::iterator2 it2_end (it1.end ());
+#else
+ typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
+ typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
+#endif
+ increment (it2, it2_end);
+ ++ it1;
+ }
+#else
+ increment (it1, it1_end, - compare);
+#endif
+ } else if (compare > 0) {
+#if BOOST_UBLAS_TYPE_CHECK
+ while (it1e.index1 () < it1.index1 ()) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename E::iterator2 it2e (it1e.begin ());
+ typename E::iterator2 it2e_end (it1e.end ());
+#else
+ typename E::iterator2 it2e (begin (it1e, iterator1_tag ()));
+ typename E::iterator2 it2e_end (end (it1e, iterator1_tag ()));
+#endif
+ increment (it2e, it2e_end);
+ ++ it1e;
+ }
+#else
+ increment (it1e, it1e_end, compare);
+#endif
+ }
+ }
+#if BOOST_UBLAS_TYPE_CHECK
+ while (it1e != it1e_end) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename E::iterator2 it2e (it1e.begin ());
+ typename E::iterator2 it2e_end (it1e.end ());
+#else
+ typename E::iterator2 it2e (begin (it1e, iterator1_tag ()));
+ typename E::iterator2 it2e_end (end (it1e, iterator1_tag ()));
+#endif
+ increment (it2e, it2e_end);
+ ++ it1e;
+ }
+ while (it1 != it1_end) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator2 it2 (it1.begin ());
+ typename M::iterator2 it2_end (it1.end ());
+#else
+ typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
+ typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
+#endif
+ increment (it2, it2_end);
+ ++ it1;
+ }
+#endif
+ }
+ // Sparse (proxy) column major case
+ template<template <class T1, class T2> class F, class R, class M, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void matrix_swap (M &m, matrix_expression<E> &e, sparse_proxy_tag, column_major_tag) {
+ typedef F<typename M::iterator1::reference, typename E::reference> functor_type;
+ typedef R conformant_restrict_type;
+ typedef typename M::size_type size_type;
+ typedef typename M::difference_type difference_type;
+
+ BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
+ BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
+
+ detail::make_conformant (m, e, column_major_tag (), conformant_restrict_type ());
+ // FIXME should be a seperate restriction for E
+ detail::make_conformant (e (), m, column_major_tag (), conformant_restrict_type ());
+
+ typename M::iterator2 it2 (m.begin2 ());
+ typename M::iterator2 it2_end (m.end2 ());
+ typename E::iterator2 it2e (e ().begin2 ());
+ typename E::iterator2 it2e_end (e ().end2 ());
+ while (it2 != it2_end && it2e != it2e_end) {
+ difference_type compare = it2.index2 () - it2e.index2 ();
+ if (compare == 0) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator1 it1 (it2.begin ());
+ typename M::iterator1 it1_end (it2.end ());
+ typename E::iterator1 it1e (it2e.begin ());
+ typename E::iterator1 it1e_end (it2e.end ());
+#else
+ typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
+ typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
+ typename E::iterator1 it1e (begin (it2e, iterator2_tag ()));
+ typename E::iterator1 it1e_end (end (it2e, iterator2_tag ()));
+#endif
+ if (it1 != it1_end && it1e != it1e_end) {
+ size_type it1_index = it1.index1 (), it1e_index = it1e.index1 ();
+ while (true) {
+ difference_type compare2 = it1_index - it1e_index;
+ if (compare2 == 0) {
+ functor_type::apply (*it1, *it1e);
+ ++ it1, ++ it1e;
+ if (it1 != it1_end && it1e != it1e_end) {
+ it1_index = it1.index1 ();
+ it1e_index = it1e.index1 ();
+ } else
+ break;
+ } else if (compare2 < 0) {
+ increment (it1, it1_end, - compare2);
+ if (it1 != it1_end)
+ it1_index = it1.index1 ();
+ else
+ break;
+ } else if (compare2 > 0) {
+ increment (it1e, it1e_end, compare2);
+ if (it1e != it1e_end)
+ it1e_index = it1e.index1 ();
+ else
+ break;
+ }
+ }
+ }
+#if BOOST_UBLAS_TYPE_CHECK
+ increment (it1e, it1e_end);
+ increment (it1, it1_end);
+#endif
+ ++ it2, ++ it2e;
+ } else if (compare < 0) {
+#if BOOST_UBLAS_TYPE_CHECK
+ while (it2.index2 () < it2e.index2 ()) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator1 it1 (it2.begin ());
+ typename M::iterator1 it1_end (it2.end ());
+#else
+ typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
+ typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
+#endif
+ increment (it1, it1_end);
+ ++ it2;
+ }
+#else
+ increment (it2, it2_end, - compare);
+#endif
+ } else if (compare > 0) {
+#if BOOST_UBLAS_TYPE_CHECK
+ while (it2e.index2 () < it2.index2 ()) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename E::iterator1 it1e (it2e.begin ());
+ typename E::iterator1 it1e_end (it2e.end ());
+#else
+ typename E::iterator1 it1e (begin (it2e, iterator2_tag ()));
+ typename E::iterator1 it1e_end (end (it2e, iterator2_tag ()));
+#endif
+ increment (it1e, it1e_end);
+ ++ it2e;
+ }
+#else
+ increment (it2e, it2e_end, compare);
+#endif
+ }
+ }
+#if BOOST_UBLAS_TYPE_CHECK
+ while (it2e != it2e_end) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename E::iterator1 it1e (it2e.begin ());
+ typename E::iterator1 it1e_end (it2e.end ());
+#else
+ typename E::iterator1 it1e (begin (it2e, iterator2_tag ()));
+ typename E::iterator1 it1e_end (end (it2e, iterator2_tag ()));
+#endif
+ increment (it1e, it1e_end);
+ ++ it2e;
+ }
+ while (it2 != it2_end) {
+#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
+ typename M::iterator1 it1 (it2.begin ());
+ typename M::iterator1 it1_end (it2.end ());
+#else
+ typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
+ typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
+#endif
+ increment (it1, it1_end);
+ ++ it2;
+ }
+#endif
+ }
+
+ // Dispatcher
+ template<template <class T1, class T2> class F, class M, class E>
+ BOOST_UBLAS_INLINE
+ void matrix_swap (M &m, matrix_expression<E> &e) {
+ typedef typename matrix_swap_traits<typename M::storage_category,
+ typename E::const_iterator1::iterator_category,
+ typename E::const_iterator2::iterator_category>::storage_category storage_category;
+ // give preference to matrix M's orientation if known
+ typedef typename boost::mpl::if_<boost::is_same<typename M::orientation_category, unknown_orientation_tag>,
+ typename E::orientation_category ,
+ typename M::orientation_category >::type orientation_category;
+ typedef basic_full<typename M::size_type> unrestricted;
+ matrix_swap<F, unrestricted> (m, e, storage_category (), orientation_category ());
+ }
+ template<template <class T1, class T2> class F, class R, class M, class E>
+ BOOST_UBLAS_INLINE
+ void matrix_swap (M &m, matrix_expression<E> &e) {
+ typedef R conformant_restrict_type;
+ typedef typename matrix_swap_traits<typename M::storage_category,
+ typename E::const_iterator1::iterator_category,
+ typename E::const_iterator2::iterator_category>::storage_category storage_category;
+ // give preference to matrix M's orientation if known
+ typedef typename boost::mpl::if_<boost::is_same<typename M::orientation_category, unknown_orientation_tag>,
+ typename E::orientation_category ,
+ typename M::orientation_category >::type orientation_category;
+ matrix_swap<F, conformant_restrict_type> (m, e, storage_category (), orientation_category ());
+ }
+
+}}}
+
+#endif
diff --git a/include/boost/numeric/ublas/detail/raw.hpp b/include/boost/numeric/ublas/detail/raw.hpp
new file mode 100644
index 0000000..c36c8b0
--- /dev/null
+++ b/include/boost/numeric/ublas/detail/raw.hpp
@@ -0,0 +1,878 @@
+//
+// Copyright (c) 2002-2003
+// Toon Knapen, Kresimir Fresl, Joerg Walter
+//
+// 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_UBLAS_RAW_
+#define _BOOST_UBLAS_RAW_
+
+namespace boost { namespace numeric { namespace ublas { namespace raw {
+
+ // We need data_const() mostly due to MSVC 6.0.
+ // But how shall we write portable code otherwise?
+
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ int size( const V &v ) ;
+
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ int size( const vector_reference<V> &v ) ;
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int size1( const M &m ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int size2( const M &m ) ;
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int size1( const matrix_reference<M> &m ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int size2( const matrix_reference<M> &m ) ;
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int leading_dimension( const M &m, row_major_tag ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int leading_dimension( const M &m, column_major_tag ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int leading_dimension( const M &m ) ;
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int leading_dimension( const matrix_reference<M> &m ) ;
+
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ int stride( const V &v ) ;
+
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ int stride( const vector_range<V> &v ) ;
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ int stride( const vector_slice<V> &v ) ;
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int stride( const matrix_row<M> &v ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int stride( const matrix_column<M> &v ) ;
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int stride1( const M &m ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int stride2( const M &m ) ;
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int stride1( const matrix_reference<M> &m ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int stride2( const matrix_reference<M> &m ) ;
+
+ template < typename T, std::size_t M, std::size_t N >
+ BOOST_UBLAS_INLINE
+ int stride1( const c_matrix<T, M, N> &m ) ;
+ template < typename T, std::size_t M, std::size_t N >
+ BOOST_UBLAS_INLINE
+ int stride2( const c_matrix<T, M, N> &m ) ;
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int stride1( const matrix_range<M> &m ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int stride1( const matrix_slice<M> &m ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int stride2( const matrix_range<M> &m ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int stride2( const matrix_slice<M> &m ) ;
+
+ template < typename MV >
+ BOOST_UBLAS_INLINE
+ typename MV::array_type::array_type::const_pointer data( const MV &mv ) ;
+ template < typename MV >
+ BOOST_UBLAS_INLINE
+ typename MV::array_type::array_type::const_pointer data_const( const MV &mv ) ;
+ template < typename MV >
+ BOOST_UBLAS_INLINE
+ typename MV::array_type::pointer data( MV &mv ) ;
+
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::array_type::const_pointer data( const vector_reference<V> &v ) ;
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::array_type::const_pointer data_const( const vector_reference<V> &v ) ;
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::pointer data( vector_reference<V> &v ) ;
+
+ template < typename T, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_vector<T, N>::array_type::array_type::const_pointer data( const c_vector<T, N> &v ) ;
+ template < typename T, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_vector<T, N>::array_type::array_type::const_pointer data_const( const c_vector<T, N> &v ) ;
+ template < typename T, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_vector<T, N>::pointer data( c_vector<T, N> &v ) ;
+
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::array_type::const_pointer data( const vector_range<V> &v ) ;
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::array_type::const_pointer data( const vector_slice<V> &v ) ;
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::array_type::const_pointer data_const( const vector_range<V> &v ) ;
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::array_type::const_pointer data_const( const vector_slice<V> &v ) ;
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::pointer data( vector_range<V> &v ) ;
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::pointer data( vector_slice<V> &v ) ;
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::array_type::const_pointer data( const matrix_reference<M> &m ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::array_type::const_pointer data_const( const matrix_reference<M> &m ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::pointer data( matrix_reference<M> &m ) ;
+
+ template < typename T, std::size_t M, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_matrix<T, M, N>::array_type::array_type::const_pointer data( const c_matrix<T, M, N> &m ) ;
+ template < typename T, std::size_t M, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_matrix<T, M, N>::array_type::array_type::const_pointer data_const( const c_matrix<T, M, N> &m ) ;
+ template < typename T, std::size_t M, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_matrix<T, M, N>::pointer data( c_matrix<T, M, N> &m ) ;
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::array_type::const_pointer data( const matrix_row<M> &v ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::array_type::const_pointer data( const matrix_column<M> &v ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::array_type::const_pointer data_const( const matrix_row<M> &v ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::array_type::const_pointer data_const( const matrix_column<M> &v ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::pointer data( matrix_row<M> &v ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::pointer data( matrix_column<M> &v ) ;
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::array_type::const_pointer data( const matrix_range<M> &m ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::array_type::const_pointer data( const matrix_slice<M> &m ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::array_type::const_pointer data_const( const matrix_range<M> &m ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::array_type::const_pointer data_const( const matrix_slice<M> &m ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::pointer data( matrix_range<M> &m ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::pointer data( matrix_slice<M> &m ) ;
+
+ template < typename MV >
+ BOOST_UBLAS_INLINE
+ typename MV::array_type::array_type::const_pointer base( const MV &mv ) ;
+
+ template < typename MV >
+ BOOST_UBLAS_INLINE
+ typename MV::array_type::array_type::const_pointer base_const( const MV &mv ) ;
+ template < typename MV >
+ BOOST_UBLAS_INLINE
+ typename MV::array_type::pointer base( MV &mv ) ;
+
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::array_type::const_pointer base( const vector_reference<V> &v ) ;
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::array_type::const_pointer base_const( const vector_reference<V> &v ) ;
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::pointer base( vector_reference<V> &v ) ;
+
+ template < typename T, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_vector<T, N>::array_type::array_type::const_pointer base( const c_vector<T, N> &v ) ;
+ template < typename T, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_vector<T, N>::array_type::array_type::const_pointer base_const( const c_vector<T, N> &v ) ;
+ template < typename T, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_vector<T, N>::pointer base( c_vector<T, N> &v ) ;
+
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::array_type::const_pointer base( const vector_range<V> &v ) ;
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::array_type::const_pointer base( const vector_slice<V> &v ) ;
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::array_type::const_pointer base_const( const vector_range<V> &v ) ;
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::array_type::const_pointer base_const( const vector_slice<V> &v ) ;
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::pointer base( vector_range<V> &v ) ;
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::pointer base( vector_slice<V> &v ) ;
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::array_type::const_pointer base( const matrix_reference<M> &m ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::array_type::const_pointer base_const( const matrix_reference<M> &m ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::pointer base( matrix_reference<M> &m ) ;
+
+ template < typename T, std::size_t M, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_matrix<T, M, N>::array_type::array_type::const_pointer base( const c_matrix<T, M, N> &m ) ;
+ template < typename T, std::size_t M, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_matrix<T, M, N>::array_type::array_type::const_pointer base_const( const c_matrix<T, M, N> &m ) ;
+ template < typename T, std::size_t M, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_matrix<T, M, N>::pointer base( c_matrix<T, M, N> &m ) ;
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::array_type::const_pointer base( const matrix_row<M> &v ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::array_type::const_pointer base( const matrix_column<M> &v ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::array_type::const_pointer base_const( const matrix_row<M> &v ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::array_type::const_pointer base_const( const matrix_column<M> &v ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::pointer base( matrix_row<M> &v ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::pointer base( matrix_column<M> &v ) ;
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::array_type::const_pointer base( const matrix_range<M> &m ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::array_type::const_pointer base( const matrix_slice<M> &m ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::array_type::const_pointer base_const( const matrix_range<M> &m ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::array_type::const_pointer base_const( const matrix_slice<M> &m ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::pointer base( matrix_range<M> &m ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::pointer base( matrix_slice<M> &m ) ;
+
+ template < typename MV >
+ BOOST_UBLAS_INLINE
+ typename MV::size_type start( const MV &mv ) ;
+
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::size_type start( const vector_range<V> &v ) ;
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::size_type start( const vector_slice<V> &v ) ;
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::size_type start( const matrix_row<M> &v ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::size_type start( const matrix_column<M> &v ) ;
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::size_type start( const matrix_range<M> &m ) ;
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::size_type start( const matrix_slice<M> &m ) ;
+
+
+
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ int size( const V &v ) {
+ return v.size() ;
+ }
+
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ int size( const vector_reference<V> &v ) {
+ return size( v ) ;
+ }
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int size1( const M &m ) {
+ return m.size1() ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int size2( const M &m ) {
+ return m.size2() ;
+ }
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int size1( const matrix_reference<M> &m ) {
+ return size1( m.expression() ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int size2( const matrix_reference<M> &m ) {
+ return size2( m.expression() ) ;
+ }
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int leading_dimension( const M &m, row_major_tag ) {
+ return m.size2() ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int leading_dimension( const M &m, column_major_tag ) {
+ return m.size1() ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int leading_dimension( const M &m ) {
+ return leading_dimension( m, typename M::orientation_category() ) ;
+ }
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int leading_dimension( const matrix_reference<M> &m ) {
+ return leading_dimension( m.expression() ) ;
+ }
+
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ int stride( const V &v ) {
+ return 1 ;
+ }
+
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ int stride( const vector_range<V> &v ) {
+ return stride( v.data() ) ;
+ }
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ int stride( const vector_slice<V> &v ) {
+ return v.stride() * stride( v.data() ) ;
+ }
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int stride( const matrix_row<M> &v ) {
+ return stride2( v.data() ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int stride( const matrix_column<M> &v ) {
+ return stride1( v.data() ) ;
+ }
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int stride1( const M &m ) {
+ typedef typename M::functor_type functor_type;
+ return functor_type::one1( m.size1(), m.size2() ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int stride2( const M &m ) {
+ typedef typename M::functor_type functor_type;
+ return functor_type::one2( m.size1(), m.size2() ) ;
+ }
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int stride1( const matrix_reference<M> &m ) {
+ return stride1( m.expression() ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int stride2( const matrix_reference<M> &m ) {
+ return stride2( m.expression() ) ;
+ }
+
+ template < typename T, std::size_t M, std::size_t N >
+ BOOST_UBLAS_INLINE
+ int stride1( const c_matrix<T, M, N> &m ) {
+ return N ;
+ }
+ template < typename T, std::size_t M, std::size_t N >
+ BOOST_UBLAS_INLINE
+ int stride2( const c_matrix<T, M, N> &m ) {
+ return 1 ;
+ }
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int stride1( const matrix_range<M> &m ) {
+ return stride1( m.data() ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int stride1( const matrix_slice<M> &m ) {
+ return m.stride1() * stride1( m.data() ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int stride2( const matrix_range<M> &m ) {
+ return stride2( m.data() ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ int stride2( const matrix_slice<M> &m ) {
+ return m.stride2() * stride2( m.data() ) ;
+ }
+
+ template < typename MV >
+ BOOST_UBLAS_INLINE
+ typename MV::array_type::array_type::array_type::const_pointer data( const MV &mv ) {
+ return &mv.data().begin()[0] ;
+ }
+ template < typename MV >
+ BOOST_UBLAS_INLINE
+ typename MV::array_type::array_type::const_pointer data_const( const MV &mv ) {
+ return &mv.data().begin()[0] ;
+ }
+ template < typename MV >
+ BOOST_UBLAS_INLINE
+ typename MV::array_type::pointer data( MV &mv ) {
+ return &mv.data().begin()[0] ;
+ }
+
+
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::array_type::const_pointer data( const vector_reference<V> &v ) {
+ return data( v.expression () ) ;
+ }
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::array_type::const_pointer data_const( const vector_reference<V> &v ) {
+ return data_const( v.expression () ) ;
+ }
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::pointer data( vector_reference<V> &v ) {
+ return data( v.expression () ) ;
+ }
+
+ template < typename T, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_vector<T, N>::array_type::array_type::const_pointer data( const c_vector<T, N> &v ) {
+ return v.data() ;
+ }
+ template < typename T, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_vector<T, N>::array_type::array_type::const_pointer data_const( const c_vector<T, N> &v ) {
+ return v.data() ;
+ }
+ template < typename T, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_vector<T, N>::pointer data( c_vector<T, N> &v ) {
+ return v.data() ;
+ }
+
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::array_type::const_pointer data( const vector_range<V> &v ) {
+ return data( v.data() ) + v.start() * stride (v.data() ) ;
+ }
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::array_type::const_pointer data( const vector_slice<V> &v ) {
+ return data( v.data() ) + v.start() * stride (v.data() ) ;
+ }
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::array_type::const_pointer data_const( const vector_range<V> &v ) {
+ return data_const( v.data() ) + v.start() * stride (v.data() ) ;
+ }
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::const_pointer data_const( const vector_slice<V> &v ) {
+ return data_const( v.data() ) + v.start() * stride (v.data() ) ;
+ }
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::pointer data( vector_range<V> &v ) {
+ return data( v.data() ) + v.start() * stride (v.data() ) ;
+ }
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::pointer data( vector_slice<V> &v ) {
+ return data( v.data() ) + v.start() * stride (v.data() ) ;
+ }
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::const_pointer data( const matrix_reference<M> &m ) {
+ return data( m.expression () ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::const_pointer data_const( const matrix_reference<M> &m ) {
+ return data_const( m.expression () ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::pointer data( matrix_reference<M> &m ) {
+ return data( m.expression () ) ;
+ }
+
+ template < typename T, std::size_t M, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_matrix<T, M, N>::array_type::const_pointer data( const c_matrix<T, M, N> &m ) {
+ return m.data() ;
+ }
+ template < typename T, std::size_t M, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_matrix<T, M, N>::array_type::const_pointer data_const( const c_matrix<T, M, N> &m ) {
+ return m.data() ;
+ }
+ template < typename T, std::size_t M, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_matrix<T, M, N>::pointer data( c_matrix<T, M, N> &m ) {
+ return m.data() ;
+ }
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::const_pointer data( const matrix_row<M> &v ) {
+ return data( v.data() ) + v.index() * stride1( v.data() ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::const_pointer data( const matrix_column<M> &v ) {
+ return data( v.data() ) + v.index() * stride2( v.data() ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::const_pointer data_const( const matrix_row<M> &v ) {
+ return data_const( v.data() ) + v.index() * stride1( v.data() ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::const_pointer data_const( const matrix_column<M> &v ) {
+ return data_const( v.data() ) + v.index() * stride2( v.data() ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::pointer data( matrix_row<M> &v ) {
+ return data( v.data() ) + v.index() * stride1( v.data() ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::pointer data( matrix_column<M> &v ) {
+ return data( v.data() ) + v.index() * stride2( v.data() ) ;
+ }
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::const_pointer data( const matrix_range<M> &m ) {
+ return data( m.data() ) + m.start1() * stride1( m.data () ) + m.start2() * stride2( m.data () ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::const_pointer data( const matrix_slice<M> &m ) {
+ return data( m.data() ) + m.start1() * stride1( m.data () ) + m.start2() * stride2( m.data () ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::const_pointer data_const( const matrix_range<M> &m ) {
+ return data_const( m.data() ) + m.start1() * stride1( m.data () ) + m.start2() * stride2( m.data () ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::const_pointer data_const( const matrix_slice<M> &m ) {
+ return data_const( m.data() ) + m.start1() * stride1( m.data () ) + m.start2() * stride2( m.data () ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::pointer data( matrix_range<M> &m ) {
+ return data( m.data() ) + m.start1() * stride1( m.data () ) + m.start2() * stride2( m.data () ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::pointer data( matrix_slice<M> &m ) {
+ return data( m.data() ) + m.start1() * stride1( m.data () ) + m.start2() * stride2( m.data () ) ;
+ }
+
+
+ template < typename MV >
+ BOOST_UBLAS_INLINE
+ typename MV::array_type::const_pointer base( const MV &mv ) {
+ return &mv.data().begin()[0] ;
+ }
+ template < typename MV >
+ BOOST_UBLAS_INLINE
+ typename MV::array_type::const_pointer base_const( const MV &mv ) {
+ return &mv.data().begin()[0] ;
+ }
+ template < typename MV >
+ BOOST_UBLAS_INLINE
+ typename MV::array_type::pointer base( MV &mv ) {
+ return &mv.data().begin()[0] ;
+ }
+
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::const_pointer base( const vector_reference<V> &v ) {
+ return base( v.expression () ) ;
+ }
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::const_pointer base_const( const vector_reference<V> &v ) {
+ return base_const( v.expression () ) ;
+ }
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::pointer base( vector_reference<V> &v ) {
+ return base( v.expression () ) ;
+ }
+
+ template < typename T, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_vector<T, N>::array_type::const_pointer base( const c_vector<T, N> &v ) {
+ return v.data() ;
+ }
+ template < typename T, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_vector<T, N>::array_type::const_pointer base_const( const c_vector<T, N> &v ) {
+ return v.data() ;
+ }
+ template < typename T, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_vector<T, N>::pointer base( c_vector<T, N> &v ) {
+ return v.data() ;
+ }
+
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::const_pointer base( const vector_range<V> &v ) {
+ return base( v.data() ) ;
+ }
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::const_pointer base( const vector_slice<V> &v ) {
+ return base( v.data() ) ;
+ }
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::const_pointer base_const( const vector_range<V> &v ) {
+ return base_const( v.data() ) ;
+ }
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::const_pointer base_const( const vector_slice<V> &v ) {
+ return base_const( v.data() ) ;
+ }
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::pointer base( vector_range<V> &v ) {
+ return base( v.data() ) ;
+ }
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::array_type::pointer base( vector_slice<V> &v ) {
+ return base( v.data() ) ;
+ }
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::const_pointer base( const matrix_reference<M> &m ) {
+ return base( m.expression () ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::const_pointer base_const( const matrix_reference<M> &m ) {
+ return base_const( m.expression () ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::pointer base( matrix_reference<M> &m ) {
+ return base( m.expression () ) ;
+ }
+
+ template < typename T, std::size_t M, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_matrix<T, M, N>::array_type::const_pointer base( const c_matrix<T, M, N> &m ) {
+ return m.data() ;
+ }
+ template < typename T, std::size_t M, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_matrix<T, M, N>::array_type::const_pointer base_const( const c_matrix<T, M, N> &m ) {
+ return m.data() ;
+ }
+ template < typename T, std::size_t M, std::size_t N >
+ BOOST_UBLAS_INLINE
+ typename c_matrix<T, M, N>::pointer base( c_matrix<T, M, N> &m ) {
+ return m.data() ;
+ }
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::const_pointer base( const matrix_row<M> &v ) {
+ return base( v.data() ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::const_pointer base( const matrix_column<M> &v ) {
+ return base( v.data() ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::const_pointer base_const( const matrix_row<M> &v ) {
+ return base_const( v.data() ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::const_pointer base_const( const matrix_column<M> &v ) {
+ return base_const( v.data() ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::pointer base( matrix_row<M> &v ) {
+ return base( v.data() ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::pointer base( matrix_column<M> &v ) {
+ return base( v.data() ) ;
+ }
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::const_pointer base( const matrix_range<M> &m ) {
+ return base( m.data() ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::const_pointer base( const matrix_slice<M> &m ) {
+ return base( m.data() ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::const_pointer base_const( const matrix_range<M> &m ) {
+ return base_const( m.data() ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::const_pointer base_const( const matrix_slice<M> &m ) {
+ return base_const( m.data() ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::pointer base( matrix_range<M> &m ) {
+ return base( m.data() ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::array_type::pointer base( matrix_slice<M> &m ) {
+ return base( m.data() ) ;
+ }
+
+ template < typename MV >
+ BOOST_UBLAS_INLINE
+ typename MV::size_type start( const MV &mv ) {
+ return 0 ;
+ }
+
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::size_type start( const vector_range<V> &v ) {
+ return v.start() * stride (v.data() ) ;
+ }
+ template < typename V >
+ BOOST_UBLAS_INLINE
+ typename V::size_type start( const vector_slice<V> &v ) {
+ return v.start() * stride (v.data() ) ;
+ }
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::size_type start( const matrix_row<M> &v ) {
+ return v.index() * stride1( v.data() ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::size_type start( const matrix_column<M> &v ) {
+ return v.index() * stride2( v.data() ) ;
+ }
+
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::size_type start( const matrix_range<M> &m ) {
+ return m.start1() * stride1( m.data () ) + m.start2() * stride2( m.data () ) ;
+ }
+ template < typename M >
+ BOOST_UBLAS_INLINE
+ typename M::size_type start( const matrix_slice<M> &m ) {
+ return m.start1() * stride1( m.data () ) + m.start2() * stride2( m.data () ) ;
+ }
+
+}}}}
+
+#endif
diff --git a/include/boost/numeric/ublas/detail/returntype_deduction.hpp b/include/boost/numeric/ublas/detail/returntype_deduction.hpp
new file mode 100644
index 0000000..030d1f6
--- /dev/null
+++ b/include/boost/numeric/ublas/detail/returntype_deduction.hpp
@@ -0,0 +1,174 @@
+/*
+ * Copyright (c) 2001-2003 Joel de Guzman
+ *
+ * Use, modification and distribution is 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)
+ */
+#ifndef _BOOST_UBLAS_NUMERICTYPE_DEDUCTION_
+#define _BOOST_UBLAS_NUMERICTYPE_DEDUCTION_
+
+// See original in boost-sandbox/boost/utility/type_deduction.hpp for comments
+
+#include <boost/mpl/vector/vector20.hpp>
+#include <boost/mpl/at.hpp>
+#include <boost/mpl/or.hpp>
+#include <boost/mpl/identity.hpp>
+#include <boost/type_traits/remove_cv.hpp>
+#include <boost/type_traits/is_same.hpp>
+#include <boost/utility/enable_if.hpp>
+
+namespace boost { namespace numeric { namespace ublas {
+
+struct error_cant_deduce_type {};
+
+ namespace type_deduction_detail
+ {
+ typedef char(&bool_value_type)[1];
+ typedef char(&float_value_type)[2];
+ typedef char(&double_value_type)[3];
+ typedef char(&long_double_value_type)[4];
+ typedef char(&char_value_type)[5];
+ typedef char(&schar_value_type)[6];
+ typedef char(&uchar_value_type)[7];
+ typedef char(&short_value_type)[8];
+ typedef char(&ushort_value_type)[9];
+ typedef char(&int_value_type)[10];
+ typedef char(&uint_value_type)[11];
+ typedef char(&long_value_type)[12];
+ typedef char(&ulong_value_type)[13];
+
+ typedef char(&x_value_type)[14];
+ typedef char(&y_value_type)[15];
+
+ typedef char(&cant_deduce_type)[16];
+
+ template <typename T, typename PlainT = typename remove_cv<T>::type>
+ struct is_basic
+ : mpl::or_<
+ typename mpl::or_<
+ is_same<PlainT, bool>
+ , is_same<PlainT, float>
+ , is_same<PlainT, double>
+ , is_same<PlainT, long double>
+ > ::type,
+ typename mpl::or_<
+ is_same<PlainT, char>
+ , is_same<PlainT, signed char>
+ , is_same<PlainT, unsigned char>
+ , is_same<PlainT, short>
+ , is_same<PlainT, unsigned short>
+ > ::type,
+ typename mpl::or_<
+ is_same<PlainT, int>
+ , is_same<PlainT, unsigned int>
+ , is_same<PlainT, long>
+ , is_same<PlainT, unsigned long>
+ > ::type
+ > {};
+
+ struct asymmetric;
+
+ template <typename X, typename Y>
+ cant_deduce_type
+ test(...); // The black hole !!!
+
+ template <typename X, typename Y>
+ bool_value_type
+ test(bool const&);
+
+ template <typename X, typename Y>
+ float_value_type
+ test(float const&);
+
+ template <typename X, typename Y>
+ double_value_type
+ test(double const&);
+
+ template <typename X, typename Y>
+ long_double_value_type
+ test(long double const&);
+
+ template <typename X, typename Y>
+ char_value_type
+ test(char const&);
+
+ template <typename X, typename Y>
+ schar_value_type
+ test(signed char const&);
+
+ template <typename X, typename Y>
+ uchar_value_type
+ test(unsigned char const&);
+
+ template <typename X, typename Y>
+ short_value_type
+ test(short const&);
+
+ template <typename X, typename Y>
+ ushort_value_type
+ test(unsigned short const&);
+
+ template <typename X, typename Y>
+ int_value_type
+ test(int const&);
+
+ template <typename X, typename Y>
+ uint_value_type
+ test(unsigned int const&);
+
+ template <typename X, typename Y>
+ long_value_type
+ test(long const&);
+
+ template <typename X, typename Y>
+ ulong_value_type
+ test(unsigned long const&);
+
+ template <typename X, typename Y>
+ typename disable_if<
+ is_basic<X>, x_value_type
+ >::type
+ test(X const&);
+
+ template <typename X, typename Y>
+ typename disable_if<
+ mpl::or_<
+ is_basic<Y>
+ , is_same<Y, asymmetric>
+ , is_same<const X, const Y>
+ >
+ , y_value_type
+ >::type
+ test(Y const&);
+
+ template <typename X, typename Y>
+ struct base_result_of
+ {
+ typedef typename remove_cv<X>::type x_type;
+ typedef typename remove_cv<Y>::type y_type;
+
+ typedef mpl::vector16<
+ mpl::identity<bool>
+ , mpl::identity<float>
+ , mpl::identity<double>
+ , mpl::identity<long double>
+ , mpl::identity<char>
+ , mpl::identity<signed char>
+ , mpl::identity<unsigned char>
+ , mpl::identity<short>
+ , mpl::identity<unsigned short>
+ , mpl::identity<int>
+ , mpl::identity<unsigned int>
+ , mpl::identity<long>
+ , mpl::identity<unsigned long>
+ , mpl::identity<x_type>
+ , mpl::identity<y_type>
+ , mpl::identity<error_cant_deduce_type>
+ >
+ types;
+ };
+
+}}} } // namespace boost::numeric::ublas ::type_deduction_detail
+
+#endif
diff --git a/include/boost/numeric/ublas/detail/temporary.hpp b/include/boost/numeric/ublas/detail/temporary.hpp
new file mode 100644
index 0000000..c2ae468
--- /dev/null
+++ b/include/boost/numeric/ublas/detail/temporary.hpp
@@ -0,0 +1,33 @@
+//
+// Copyright (c) 2000-2002
+// Joerg Walter, Mathias Koch
+//
+// 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)
+//
+// The authors gratefully acknowledge the support of
+// GeNeSys mbH & Co. KG in producing this work.
+//
+
+#ifndef _BOOST_UBLAS_TEMPORARY_
+#define _BOOST_UBLAS_TEMPORARY_
+
+
+namespace boost { namespace numeric { namespace ublas {
+
+/// For the creation of temporary vectors in the assignment of proxies
+template <class M>
+struct vector_temporary_traits {
+ typedef typename M::vector_temporary_type type ;
+};
+
+/// For the creation of temporary vectors in the assignment of proxies
+template <class M>
+struct matrix_temporary_traits {
+ typedef typename M::matrix_temporary_type type ;
+};
+
+} } }
+
+#endif
diff --git a/include/boost/numeric/ublas/detail/vector_assign.hpp b/include/boost/numeric/ublas/detail/vector_assign.hpp
new file mode 100644
index 0000000..e612876
--- /dev/null
+++ b/include/boost/numeric/ublas/detail/vector_assign.hpp
@@ -0,0 +1,570 @@
+//
+// Copyright (c) 2000-2002
+// Joerg Walter, Mathias Koch
+//
+// 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)
+//
+// The authors gratefully acknowledge the support of
+// GeNeSys mbH & Co. KG in producing this work.
+//
+
+#ifndef _BOOST_UBLAS_VECTOR_ASSIGN_
+#define _BOOST_UBLAS_VECTOR_ASSIGN_
+
+#include <boost/numeric/ublas/functional.hpp> // scalar_assign
+// Required for make_conformant storage
+#include <vector>
+
+// Iterators based on ideas of Jeremy Siek
+
+namespace boost { namespace numeric { namespace ublas {
+namespace detail {
+
+ // Weak equality check - useful to compare equality two arbitary vector expression results.
+ // Since the actual expressions are unknown, we check for and arbitary error bound
+ // on the relative error.
+ // For a linear expression the infinity norm makes sense as we do not know how the elements will be
+ // combined in the expression. False positive results are inevitable for arbirary expressions!
+ template<class E1, class E2, class S>
+ BOOST_UBLAS_INLINE
+ bool equals (const vector_expression<E1> &e1, const vector_expression<E2> &e2, S epsilon, S min_norm) {
+ return norm_inf (e1 - e2) <= epsilon *
+ std::max<S> (std::max<S> (norm_inf (e1), norm_inf (e2)), min_norm);
+ }
+
+ template<class E1, class E2>
+ BOOST_UBLAS_INLINE
+ bool expression_type_check (const vector_expression<E1> &e1, const vector_expression<E2> &e2) {
+ typedef typename type_traits<typename promote_traits<typename E1::value_type,
+ typename E2::value_type>::promote_type>::real_type real_type;
+ return equals (e1, e2, BOOST_UBLAS_TYPE_CHECK_EPSILON, BOOST_UBLAS_TYPE_CHECK_MIN);
+ }
+
+
+ // Make sparse proxies conformant
+ template<class V, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void make_conformant (V &v, const vector_expression<E> &e) {
+ BOOST_UBLAS_CHECK (v.size () == e ().size (), bad_size ());
+ typedef typename V::size_type size_type;
+ typedef typename V::difference_type difference_type;
+ typedef typename V::value_type value_type;
+ // FIXME unbounded_array with push_back maybe better
+ std::vector<size_type> index;
+ typename V::iterator it (v.begin ());
+ typename V::iterator it_end (v.end ());
+ typename E::const_iterator ite (e ().begin ());
+ typename E::const_iterator ite_end (e ().end ());
+ if (it != it_end && ite != ite_end) {
+ size_type it_index = it.index (), ite_index = ite.index ();
+ while (true) {
+ difference_type compare = it_index - ite_index;
+ if (compare == 0) {
+ ++ it, ++ ite;
+ if (it != it_end && ite != ite_end) {
+ it_index = it.index ();
+ ite_index = ite.index ();
+ } else
+ break;
+ } else if (compare < 0) {
+ increment (it, it_end, - compare);
+ if (it != it_end)
+ it_index = it.index ();
+ else
+ break;
+ } else if (compare > 0) {
+ if (*ite != value_type/*zero*/())
+ index.push_back (ite.index ());
+ ++ ite;
+ if (ite != ite_end)
+ ite_index = ite.index ();
+ else
+ break;
+ }
+ }
+ }
+
+ while (ite != ite_end) {
+ if (*ite != value_type/*zero*/())
+ index.push_back (ite.index ());
+ ++ ite;
+ }
+ for (size_type k = 0; k < index.size (); ++ k)
+ v (index [k]) = value_type/*zero*/();
+ }
+
+}//namespace detail
+
+
+ // Explicitly iterating
+ template<template <class T1, class T2> class F, class V, class T>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void iterating_vector_assign_scalar (V &v, const T &t) {
+ typedef F<typename V::iterator::reference, T> functor_type;
+ typedef typename V::difference_type difference_type;
+ difference_type size (v.size ());
+ typename V::iterator it (v.begin ());
+ BOOST_UBLAS_CHECK (v.end () - it == size, bad_size ());
+#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
+ while (-- size >= 0)
+ functor_type::apply (*it, t), ++ it;
+#else
+ DD (size, 4, r, (functor_type::apply (*it, t), ++ it));
+#endif
+ }
+ // Explicitly case
+ template<template <class T1, class T2> class F, class V, class T>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void indexing_vector_assign_scalar (V &v, const T &t) {
+ typedef F<typename V::reference, T> functor_type;
+ typedef typename V::size_type size_type;
+ size_type size (v.size ());
+#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
+ for (size_type i = 0; i < size; ++ i)
+ functor_type::apply (v (i), t);
+#else
+ size_type i (0);
+ DD (size, 4, r, (functor_type::apply (v (i), t), ++ i));
+#endif
+ }
+
+ // Dense (proxy) case
+ template<template <class T1, class T2> class F, class V, class T>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void vector_assign_scalar (V &v, const T &t, dense_proxy_tag) {
+#ifdef BOOST_UBLAS_USE_INDEXING
+ indexing_vector_assign_scalar<F> (v, t);
+#elif BOOST_UBLAS_USE_ITERATING
+ iterating_vector_assign_scalar<F> (v, t);
+#else
+ typedef typename V::size_type size_type;
+ size_type size (v.size ());
+ if (size >= BOOST_UBLAS_ITERATOR_THRESHOLD)
+ iterating_vector_assign_scalar<F> (v, t);
+ else
+ indexing_vector_assign_scalar<F> (v, t);
+#endif
+ }
+ // Packed (proxy) case
+ template<template <class T1, class T2> class F, class V, class T>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void vector_assign_scalar (V &v, const T &t, packed_proxy_tag) {
+ typedef F<typename V::iterator::reference, T> functor_type;
+ typedef typename V::difference_type difference_type;
+ typename V::iterator it (v.begin ());
+ difference_type size (v.end () - it);
+ while (-- size >= 0)
+ functor_type::apply (*it, t), ++ it;
+ }
+ // Sparse (proxy) case
+ template<template <class T1, class T2> class F, class V, class T>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void vector_assign_scalar (V &v, const T &t, sparse_proxy_tag) {
+ typedef F<typename V::iterator::reference, T> functor_type;
+ typename V::iterator it (v.begin ());
+ typename V::iterator it_end (v.end ());
+ while (it != it_end)
+ functor_type::apply (*it, t), ++ it;
+ }
+
+ // Dispatcher
+ template<template <class T1, class T2> class F, class V, class T>
+ BOOST_UBLAS_INLINE
+ void vector_assign_scalar (V &v, const T &t) {
+ typedef typename V::storage_category storage_category;
+ vector_assign_scalar<F> (v, t, storage_category ());
+ }
+
+ template<class SC, bool COMPUTED, class RI>
+ struct vector_assign_traits {
+ typedef SC storage_category;
+ };
+
+ template<bool COMPUTED>
+ struct vector_assign_traits<dense_tag, COMPUTED, packed_random_access_iterator_tag> {
+ typedef packed_tag storage_category;
+ };
+ template<>
+ struct vector_assign_traits<dense_tag, false, sparse_bidirectional_iterator_tag> {
+ typedef sparse_tag storage_category;
+ };
+ template<>
+ struct vector_assign_traits<dense_tag, true, sparse_bidirectional_iterator_tag> {
+ typedef sparse_proxy_tag storage_category;
+ };
+
+ template<bool COMPUTED>
+ struct vector_assign_traits<dense_proxy_tag, COMPUTED, packed_random_access_iterator_tag> {
+ typedef packed_proxy_tag storage_category;
+ };
+ template<>
+ struct vector_assign_traits<dense_proxy_tag, false, sparse_bidirectional_iterator_tag> {
+ typedef sparse_proxy_tag storage_category;
+ };
+ template<>
+ struct vector_assign_traits<dense_proxy_tag, true, sparse_bidirectional_iterator_tag> {
+ typedef sparse_proxy_tag storage_category;
+ };
+
+ template<>
+ struct vector_assign_traits<packed_tag, false, sparse_bidirectional_iterator_tag> {
+ typedef sparse_tag storage_category;
+ };
+ template<>
+ struct vector_assign_traits<packed_tag, true, sparse_bidirectional_iterator_tag> {
+ typedef sparse_proxy_tag storage_category;
+ };
+
+ template<bool COMPUTED>
+ struct vector_assign_traits<packed_proxy_tag, COMPUTED, sparse_bidirectional_iterator_tag> {
+ typedef sparse_proxy_tag storage_category;
+ };
+
+ template<>
+ struct vector_assign_traits<sparse_tag, true, dense_random_access_iterator_tag> {
+ typedef sparse_proxy_tag storage_category;
+ };
+ template<>
+ struct vector_assign_traits<sparse_tag, true, packed_random_access_iterator_tag> {
+ typedef sparse_proxy_tag storage_category;
+ };
+ template<>
+ struct vector_assign_traits<sparse_tag, true, sparse_bidirectional_iterator_tag> {
+ typedef sparse_proxy_tag storage_category;
+ };
+
+ // Explicitly iterating
+ template<template <class T1, class T2> class F, class V, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void iterating_vector_assign (V &v, const vector_expression<E> &e) {
+ typedef F<typename V::iterator::reference, typename E::value_type> functor_type;
+ typedef typename V::difference_type difference_type;
+ difference_type size (BOOST_UBLAS_SAME (v.size (), e ().size ()));
+ typename V::iterator it (v.begin ());
+ BOOST_UBLAS_CHECK (v.end () - it == size, bad_size ());
+ typename E::const_iterator ite (e ().begin ());
+ BOOST_UBLAS_CHECK (e ().end () - ite == size, bad_size ());
+#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
+ while (-- size >= 0)
+ functor_type::apply (*it, *ite), ++ it, ++ ite;
+#else
+ DD (size, 2, r, (functor_type::apply (*it, *ite), ++ it, ++ ite));
+#endif
+ }
+ // Explicitly indexing
+ template<template <class T1, class T2> class F, class V, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void indexing_vector_assign (V &v, const vector_expression<E> &e) {
+ typedef F<typename V::reference, typename E::value_type> functor_type;
+ typedef typename V::size_type size_type;
+ size_type size (BOOST_UBLAS_SAME (v.size (), e ().size ()));
+#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
+ for (size_type i = 0; i < size; ++ i)
+ functor_type::apply (v (i), e () (i));
+#else
+ size_type i (0);
+ DD (size, 2, r, (functor_type::apply (v (i), e () (i)), ++ i));
+#endif
+ }
+
+ // Dense (proxy) case
+ template<template <class T1, class T2> class F, class V, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void vector_assign (V &v, const vector_expression<E> &e, dense_proxy_tag) {
+#ifdef BOOST_UBLAS_USE_INDEXING
+ indexing_vector_assign<F> (v, e);
+#elif BOOST_UBLAS_USE_ITERATING
+ iterating_vector_assign<F> (v, e);
+#else
+ typedef typename V::size_type size_type;
+ size_type size (BOOST_UBLAS_SAME (v.size (), e ().size ()));
+ if (size >= BOOST_UBLAS_ITERATOR_THRESHOLD)
+ iterating_vector_assign<F> (v, e);
+ else
+ indexing_vector_assign<F> (v, e);
+#endif
+ }
+ // Packed (proxy) case
+ template<template <class T1, class T2> class F, class V, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void vector_assign (V &v, const vector_expression<E> &e, packed_proxy_tag) {
+ BOOST_UBLAS_CHECK (v.size () == e ().size (), bad_size ());
+ typedef F<typename V::iterator::reference, typename E::value_type> functor_type;
+ typedef typename V::difference_type difference_type;
+ typedef typename V::value_type value_type;
+#if BOOST_UBLAS_TYPE_CHECK
+ vector<value_type> cv (v.size ());
+ indexing_vector_assign<scalar_assign> (cv, v);
+ indexing_vector_assign<F> (cv, e);
+#endif
+ typename V::iterator it (v.begin ());
+ typename V::iterator it_end (v.end ());
+ typename E::const_iterator ite (e ().begin ());
+ typename E::const_iterator ite_end (e ().end ());
+ difference_type it_size (it_end - it);
+ difference_type ite_size (ite_end - ite);
+ if (it_size > 0 && ite_size > 0) {
+ difference_type size ((std::min) (difference_type (it.index () - ite.index ()), ite_size));
+ if (size > 0) {
+ ite += size;
+ ite_size -= size;
+ }
+ }
+ if (it_size > 0 && ite_size > 0) {
+ difference_type size ((std::min) (difference_type (ite.index () - it.index ()), it_size));
+ if (size > 0) {
+ it_size -= size;
+ if (!functor_type::computed) {
+ while (-- size >= 0) // zeroing
+ functor_type::apply (*it, value_type/*zero*/()), ++ it;
+ } else {
+ it += size;
+ }
+ }
+ }
+ difference_type size ((std::min) (it_size, ite_size));
+ it_size -= size;
+ ite_size -= size;
+ while (-- size >= 0)
+ functor_type::apply (*it, *ite), ++ it, ++ ite;
+ size = it_size;
+ if (!functor_type::computed) {
+ while (-- size >= 0) // zeroing
+ functor_type::apply (*it, value_type/*zero*/()), ++ it;
+ } else {
+ it += size;
+ }
+#if BOOST_UBLAS_TYPE_CHECK
+ if (! disable_type_check<bool>::value)
+ BOOST_UBLAS_CHECK (detail::expression_type_check (v, cv),
+ external_logic ("external logic or bad condition of inputs"));
+#endif
+ }
+ // Sparse case
+ template<template <class T1, class T2> class F, class V, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void vector_assign (V &v, const vector_expression<E> &e, sparse_tag) {
+ BOOST_UBLAS_CHECK (v.size () == e ().size (), bad_size ());
+ typedef F<typename V::iterator::reference, typename E::value_type> functor_type;
+ BOOST_STATIC_ASSERT ((!functor_type::computed));
+ typedef typename V::value_type value_type;
+#if BOOST_UBLAS_TYPE_CHECK
+ vector<value_type> cv (v.size ());
+ indexing_vector_assign<scalar_assign> (cv, v);
+ indexing_vector_assign<F> (cv, e);
+#endif
+ v.clear ();
+ typename E::const_iterator ite (e ().begin ());
+ typename E::const_iterator ite_end (e ().end ());
+ while (ite != ite_end) {
+ value_type t (*ite);
+ if (t != value_type/*zero*/())
+ v.insert_element (ite.index (), t);
+ ++ ite;
+ }
+#if BOOST_UBLAS_TYPE_CHECK
+ if (! disable_type_check<bool>::value)
+ BOOST_UBLAS_CHECK (detail::expression_type_check (v, cv),
+ external_logic ("external logic or bad condition of inputs"));
+#endif
+ }
+ // Sparse proxy or functional case
+ template<template <class T1, class T2> class F, class V, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void vector_assign (V &v, const vector_expression<E> &e, sparse_proxy_tag) {
+ BOOST_UBLAS_CHECK (v.size () == e ().size (), bad_size ());
+ typedef F<typename V::iterator::reference, typename E::value_type> functor_type;
+ typedef typename V::size_type size_type;
+ typedef typename V::difference_type difference_type;
+ typedef typename V::value_type value_type;
+
+#if BOOST_UBLAS_TYPE_CHECK
+ vector<value_type> cv (v.size ());
+ indexing_vector_assign<scalar_assign> (cv, v);
+ indexing_vector_assign<F> (cv, e);
+#endif
+ detail::make_conformant (v, e);
+
+ typename V::iterator it (v.begin ());
+ typename V::iterator it_end (v.end ());
+ typename E::const_iterator ite (e ().begin ());
+ typename E::const_iterator ite_end (e ().end ());
+ if (it != it_end && ite != ite_end) {
+ size_type it_index = it.index (), ite_index = ite.index ();
+ while (true) {
+ difference_type compare = it_index - ite_index;
+ if (compare == 0) {
+ functor_type::apply (*it, *ite);
+ ++ it, ++ ite;
+ if (it != it_end && ite != ite_end) {
+ it_index = it.index ();
+ ite_index = ite.index ();
+ } else
+ break;
+ } else if (compare < 0) {
+ if (!functor_type::computed) {
+ functor_type::apply (*it, value_type/*zero*/());
+ ++ it;
+ } else
+ increment (it, it_end, - compare);
+ if (it != it_end)
+ it_index = it.index ();
+ else
+ break;
+ } else if (compare > 0) {
+ increment (ite, ite_end, compare);
+ if (ite != ite_end)
+ ite_index = ite.index ();
+ else
+ break;
+ }
+ }
+ }
+
+ if (!functor_type::computed) {
+ while (it != it_end) { // zeroing
+ functor_type::apply (*it, value_type/*zero*/());
+ ++ it;
+ }
+ } else {
+ it = it_end;
+ }
+#if BOOST_UBLAS_TYPE_CHECK
+ if (! disable_type_check<bool>::value)
+ BOOST_UBLAS_CHECK (detail::expression_type_check (v, cv),
+ external_logic ("external logic or bad condition of inputs"));
+#endif
+ }
+
+ // Dispatcher
+ template<template <class T1, class T2> class F, class V, class E>
+ BOOST_UBLAS_INLINE
+ void vector_assign (V &v, const vector_expression<E> &e) {
+ typedef typename vector_assign_traits<typename V::storage_category,
+ F<typename V::reference, typename E::value_type>::computed,
+ typename E::const_iterator::iterator_category>::storage_category storage_category;
+ vector_assign<F> (v, e, storage_category ());
+ }
+
+ template<class SC, class RI>
+ struct vector_swap_traits {
+ typedef SC storage_category;
+ };
+
+ template<>
+ struct vector_swap_traits<dense_proxy_tag, sparse_bidirectional_iterator_tag> {
+ typedef sparse_proxy_tag storage_category;
+ };
+
+ template<>
+ struct vector_swap_traits<packed_proxy_tag, sparse_bidirectional_iterator_tag> {
+ typedef sparse_proxy_tag storage_category;
+ };
+
+ // Dense (proxy) case
+ template<template <class T1, class T2> class F, class V, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void vector_swap (V &v, vector_expression<E> &e, dense_proxy_tag) {
+ typedef F<typename V::iterator::reference, typename E::iterator::reference> functor_type;
+ typedef typename V::difference_type difference_type;
+ difference_type size (BOOST_UBLAS_SAME (v.size (), e ().size ()));
+ typename V::iterator it (v.begin ());
+ typename E::iterator ite (e ().begin ());
+ while (-- size >= 0)
+ functor_type::apply (*it, *ite), ++ it, ++ ite;
+ }
+ // Packed (proxy) case
+ template<template <class T1, class T2> class F, class V, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void vector_swap (V &v, vector_expression<E> &e, packed_proxy_tag) {
+ typedef F<typename V::iterator::reference, typename E::iterator::reference> functor_type;
+ typedef typename V::difference_type difference_type;
+ typename V::iterator it (v.begin ());
+ typename V::iterator it_end (v.end ());
+ typename E::iterator ite (e ().begin ());
+ typename E::iterator ite_end (e ().end ());
+ difference_type it_size (it_end - it);
+ difference_type ite_size (ite_end - ite);
+ if (it_size > 0 && ite_size > 0) {
+ difference_type size ((std::min) (difference_type (it.index () - ite.index ()), ite_size));
+ if (size > 0) {
+ ite += size;
+ ite_size -= size;
+ }
+ }
+ if (it_size > 0 && ite_size > 0) {
+ difference_type size ((std::min) (difference_type (ite.index () - it.index ()), it_size));
+ if (size > 0)
+ it_size -= size;
+ }
+ difference_type size ((std::min) (it_size, ite_size));
+ it_size -= size;
+ ite_size -= size;
+ while (-- size >= 0)
+ functor_type::apply (*it, *ite), ++ it, ++ ite;
+ }
+ // Sparse proxy case
+ template<template <class T1, class T2> class F, class V, class E>
+ // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
+ void vector_swap (V &v, vector_expression<E> &e, sparse_proxy_tag) {
+ BOOST_UBLAS_CHECK (v.size () == e ().size (), bad_size ());
+ typedef F<typename V::iterator::reference, typename E::iterator::reference> functor_type;
+ typedef typename V::size_type size_type;
+ typedef typename V::difference_type difference_type;
+
+ detail::make_conformant (v, e);
+ // FIXME should be a seperate restriction for E
+ detail::make_conformant (e (), v);
+
+ typename V::iterator it (v.begin ());
+ typename V::iterator it_end (v.end ());
+ typename E::iterator ite (e ().begin ());
+ typename E::iterator ite_end (e ().end ());
+ if (it != it_end && ite != ite_end) {
+ size_type it_index = it.index (), ite_index = ite.index ();
+ while (true) {
+ difference_type compare = it_index - ite_index;
+ if (compare == 0) {
+ functor_type::apply (*it, *ite);
+ ++ it, ++ ite;
+ if (it != it_end && ite != ite_end) {
+ it_index = it.index ();
+ ite_index = ite.index ();
+ } else
+ break;
+ } else if (compare < 0) {
+ increment (it, it_end, - compare);
+ if (it != it_end)
+ it_index = it.index ();
+ else
+ break;
+ } else if (compare > 0) {
+ increment (ite, ite_end, compare);
+ if (ite != ite_end)
+ ite_index = ite.index ();
+ else
+ break;
+ }
+ }
+ }
+
+#if BOOST_UBLAS_TYPE_CHECK
+ increment (ite, ite_end);
+ increment (it, it_end);
+#endif
+ }
+
+ // Dispatcher
+ template<template <class T1, class T2> class F, class V, class E>
+ BOOST_UBLAS_INLINE
+ void vector_swap (V &v, vector_expression<E> &e) {
+ typedef typename vector_swap_traits<typename V::storage_category,
+ typename E::const_iterator::iterator_category>::storage_category storage_category;
+ vector_swap<F> (v, e, storage_category ());
+ }
+
+}}}
+
+#endif