diff --git a/include/boost/numeric/ublas/detail/concepts.hpp b/include/boost/numeric/ublas/detail/concepts.hpp index b8b0ef1e..df552b3d 100644 --- a/include/boost/numeric/ublas/detail/concepts.hpp +++ b/include/boost/numeric/ublas/detail/concepts.hpp @@ -898,6 +898,9 @@ namespace boost { namespace numeric { namespace ublas { #define INTERNAL_EXPRESSION #endif + // TODO enable this for development + // #define VIEW_CONCEPTS + // Element value type for tests typedef float T; @@ -969,6 +972,17 @@ namespace boost { namespace numeric { namespace ublas { } #endif +#ifdef VIEW_CONCEPTS + // read only vectors + { + typedef vector_view container_model; + function_requires< RandomAccessContainerConcept >(); + function_requires< VectorConcept >(); + function_requires< IndexedRandomAccess1DIteratorConcept >(); + function_requires< IndexedRandomAccess1DIteratorConcept >(); + } +#endif + // Vector #if defined (INTERNAL_VECTOR) || defined (INTERNAL_VECTOR_DENSE) { diff --git a/include/boost/numeric/ublas/experimental/sparse_view.hpp b/include/boost/numeric/ublas/experimental/sparse_view.hpp new file mode 100644 index 00000000..688d5145 --- /dev/null +++ b/include/boost/numeric/ublas/experimental/sparse_view.hpp @@ -0,0 +1,316 @@ +// +// Copyright (c) 2009 +// Gunter Winkler +// +// 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_SPARSE_VIEW_ +#define _BOOST_UBLAS_SPARSE_VIEW_ + +#include +#include +#if BOOST_UBLAS_TYPE_CHECK +#include +#endif + +#include +#include + +namespace boost { namespace numeric { namespace ublas { + + // view a chunk of memory as ublas array + + template < class T > + class c_array_view + : public storage_array< c_array_view > { + private: + typedef c_array_view self_type; + typedef T * pointer; + + public: + // TODO: think about a const pointer + typedef const pointer array_type; + + typedef std::size_t size_type; + typedef std::ptrdiff_t difference_type; + + typedef T value_type; + typedef const T &const_reference; + typedef const T *const_pointer; + + typedef const_pointer const_iterator; + typedef std::reverse_iterator const_reverse_iterator; + + // + // typedefs required by vector concept + // + + typedef dense_tag storage_category; + typedef const vector_reference const_closure_type; + + c_array_view(size_type size, array_type data) : + size_(size), data_(data) + {} + + ~c_array_view() + {} + + // + // immutable methods of container concept + // + + BOOST_UBLAS_INLINE + size_type size () const { + return size_; + } + + BOOST_UBLAS_INLINE + const_reference operator [] (size_type i) const { + BOOST_UBLAS_CHECK (i < size_, bad_index ()); + return data_ [i]; + } + + BOOST_UBLAS_INLINE + const_iterator begin () const { + return data_; + } + BOOST_UBLAS_INLINE + const_iterator end () const { + return data_ + size_; + } + + BOOST_UBLAS_INLINE + const_reverse_iterator rbegin () const { + return const_reverse_iterator (end ()); + } + BOOST_UBLAS_INLINE + const_reverse_iterator rend () const { + return const_reverse_iterator (begin ()); + } + + private: + size_type size_; + array_type data_; + }; + + + /** \brief Present existing arrays as compressed array based + * sparse matrix. + * This class provides CRS / CCS storage layout. + * + * see also http://www.netlib.org/utk/papers/templates/node90.html + * + * \param L layout type, either row_major or column_major + * \param IB index base, use 0 for C indexing and 1 for + * FORTRAN indexing of the internal index arrays. This + * does not affect the operator()(int,int) where the first + * row/column has always index 0. + * \param IA index array type, e.g., int[] + * \param TA value array type, e.g., double[] + */ + template + class compressed_matrix_view: + public matrix_expression > { + + public: + typedef typename vector_view_traits::value_type value_type; + + private: + typedef value_type &true_reference; + typedef value_type *pointer; + typedef const value_type *const_pointer; + typedef L layout_type; + typedef compressed_matrix_view self_type; + + public: +#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS + using matrix_expression::operator (); +#endif + // ISSUE require type consistency check + // is_convertable (IA::size_type, TA::size_type) + typedef typename boost::remove_cv::value_type>::type index_type; + // for compatibility, should be removed some day ... + typedef index_type size_type; + // size_type for the data arrays. + typedef typename vector_view_traits::size_type array_size_type; + typedef typename vector_view_traits::difference_type difference_type; + typedef const value_type & const_reference; + + // do NOT define reference type, because class is read only + // typedef value_type & reference; + + typedef IA rowptr_array_type; + typedef JA index_array_type; + typedef TA value_array_type; + typedef const matrix_reference const_closure_type; + typedef matrix_reference closure_type; + + // FIXME: define a corresponding temporary type + // typedef compressed_vector vector_temporary_type; + + // FIXME: define a corresponding temporary type + // typedef self_type matrix_temporary_type; + + typedef sparse_tag storage_category; + typedef typename L::orientation_category orientation_category; + + // + // private types for internal use + // + + private: + typedef typename vector_view_traits::const_iterator const_subiterator_type; + + // + // Construction and destruction + // + private: + /// private default constructor because data must be filled by caller + BOOST_UBLAS_INLINE + compressed_matrix_view () { } + + public: + BOOST_UBLAS_INLINE + compressed_matrix_view (index_type n_rows, index_type n_cols, array_size_type nnz + , const rowptr_array_type & iptr + , const index_array_type & jptr + , const value_array_type & values): + matrix_expression (), + size1_ (n_rows), size2_ (n_cols), + nnz_ (nnz), + index1_data_ (iptr), + index2_data_ (jptr), + value_data_ (values) { + storage_invariants (); + } + + BOOST_UBLAS_INLINE + compressed_matrix_view(const compressed_matrix_view& o) : + size1_(size1_), size2_(size2_), + nnz_(nnz_), + index1_data_(index1_data_), + index2_data_(index2_data_), + value_data_(value_data_) + {} + + // + // implement immutable iterator types + // + + class const_iterator1 {}; + class const_iterator2 {}; + + typedef reverse_iterator_base1 const_reverse_iterator1; + typedef reverse_iterator_base2 const_reverse_iterator2; + + // + // implement all read only methods for the matrix expression concept + // + + //! return the number of rows + index_type size1() const { + return size1_; + } + + //! return the number of columns + index_type size2() const { + return size2_; + } + + //! return value at position (i,j) + value_type operator()(index_type i, index_type j) const { + const_pointer p = find_element(i,j); + if (!p) { + return zero_; + } else { + return *p; + } + } + + + private: + // + // private helper functions + // + + const_pointer find_element (index_type i, index_type j) const { + index_type element1 (layout_type::index_M (i, j)); + index_type element2 (layout_type::index_m (i, j)); + + const array_size_type itv = zero_based( index1_data_[element1] ); + const array_size_type itv_next = zero_based( index1_data_[element1+1] ); + + const_subiterator_type it_start = boost::next(vector_view_traits::begin(index2_data_),itv); + const_subiterator_type it_end = boost::next(vector_view_traits::begin(index2_data_),itv_next); + const_subiterator_type it = find_index_in_row(it_start, it_end, element2) ; + + if (it == it_end || *it != k_based (element2)) + return 0; + return &value_data_ [it - vector_view_traits::begin(index2_data_)]; + } + + const_subiterator_type find_index_in_row(const_subiterator_type it_start + , const_subiterator_type it_end + , index_type index) const { + return std::lower_bound( it_start + , it_end + , k_based (index) ); + } + + + private: + void storage_invariants () const { + BOOST_UBLAS_CHECK (index1_data_ [layout_type::size_M (size1_, size2_)] == k_based (nnz_), external_logic ()); + } + + index_type size1_; + index_type size2_; + + array_size_type nnz_; + + const rowptr_array_type & index1_data_; + const index_array_type & index2_data_; + const value_array_type & value_data_; + + static const value_type zero_; + + BOOST_UBLAS_INLINE + static index_type zero_based (index_type k_based_index) { + return k_based_index - IB; + } + BOOST_UBLAS_INLINE + static index_type k_based (index_type zero_based_index) { + return zero_based_index + IB; + } + + friend class iterator1; + friend class iterator2; + friend class const_iterator1; + friend class const_iterator2; + }; + + template + const typename compressed_matrix_view::value_type + compressed_matrix_view::zero_ = value_type/*zero*/(); + + + template + compressed_matrix_view + make_compressed_matrix_view(typename vector_view_traits::value_type n_rows + , typename vector_view_traits::value_type n_cols + , typename vector_view_traits::size_type nnz + , const IA & ia + , const JA & ja + , const TA & ta) { + + return compressed_matrix_view(n_rows, n_cols, nnz, ia, ja, ta); + + } + +}}} + +#endif diff --git a/include/boost/numeric/ublas/traits.hpp b/include/boost/numeric/ublas/traits.hpp index 75f6925b..e08c517c 100644 --- a/include/boost/numeric/ublas/traits.hpp +++ b/include/boost/numeric/ublas/traits.hpp @@ -512,11 +512,11 @@ namespace boost { namespace numeric { namespace ublas { } - /** \brief Traits class to extract type information from a matrix or vector CONTAINER. + /** \brief Traits class to extract type information from a constant matrix or vector CONTAINER. * */ template < class E > - struct container_traits { + struct container_view_traits { /// type of indices typedef typename E::size_type size_type; /// type of differences of indices @@ -527,36 +527,127 @@ namespace boost { namespace numeric { namespace ublas { /// type of elements typedef typename E::value_type value_type; - /// reference to an element - typedef typename E::reference reference; /// const reference to an element typedef typename E::const_reference const_reference; - /// type used in expressions to mark a reference to this class (usually a container_reference or the class itself) - typedef typename E::closure_type closure_type; /// type used in expressions to mark a reference to this class (usually a const container_reference or the class itself) typedef typename E::const_closure_type const_closure_type; }; + /** \brief Traits class to extract additional type information from a mutable matrix or vector CONTAINER. + * + */ + template < class E > + struct mutable_container_traits { + /// reference to an element + typedef typename E::reference reference; + + /// type used in expressions to mark a reference to this class (usually a container_reference or the class itself) + typedef typename E::closure_type closure_type; + }; + + /** \brief Traits class to extract type information from a matrix or vector CONTAINER. + * + */ + template < class E > + struct container_traits + : container_view_traits, mutable_container_traits { + + }; + + + /** \brief Traits class to extract type information from a constant MATRIX. + * + */ + template < class MATRIX > + struct matrix_view_traits : container_view_traits { + + /// orientation of the matrix, either \c row_major_tag, \c column_major_tag or \c unknown_orientation_tag + typedef typename MATRIX::orientation_category orientation_category; + + /// row iterator for the matrix + typedef typename MATRIX::const_iterator1 const_iterator1; + + /// column iterator for the matrix + typedef typename MATRIX::const_iterator2 const_iterator2; + }; + + /** \brief Traits class to extract additional type information from a mutable MATRIX. + * + */ + template < class MATRIX > + struct mutable_matrix_traits + : mutable_container_traits { + + /// row iterator for the matrix + typedef typename MATRIX::iterator1 iterator1; + + /// column iterator for the matrix + typedef typename MATRIX::iterator2 iterator2; + }; + + /** \brief Traits class to extract type information from a MATRIX. * */ template < class MATRIX > - struct matrix_traits : container_traits { - - /// orientation of the matrix, either \c row_major_tag, \c column_major_tag or \c unknown_orientation_tag - typedef typename MATRIX::orientation_category orientation_category; - + struct matrix_traits + : matrix_view_traits , mutable_matrix_traits { }; /** \brief Traits class to extract type information from a VECTOR. * */ template < class VECTOR > - struct vector_traits : container_traits { + struct vector_view_traits : container_view_traits { + + /// iterator for the VECTOR + typedef typename VECTOR::const_iterator const_iterator; + + /// iterator pointing to the first element + static + const_iterator begin(const VECTOR & v) { + return v.begin(); + } + /// iterator pointing behind the last element + static + const_iterator end(const VECTOR & v) { + return v.end(); + } }; + /** \brief Traits class to extract type information from a VECTOR. + * + */ + template < class VECTOR > + struct mutable_vector_traits : mutable_container_traits { + /// iterator for the VECTOR + typedef typename VECTOR::iterator iterator; + + /// iterator pointing to the first element + static + iterator begin(VECTOR & v) { + return v.begin(); + } + + /// iterator pointing behind the last element + static + iterator end(VECTOR & v) { + return v.end(); + } + }; + + /** \brief Traits class to extract type information from a VECTOR. + * + */ + template < class VECTOR > + struct vector_traits + : vector_view_traits , mutable_vector_traits { + }; + + + template < class T, int M, int N > struct matrix_traits < T[M][N] > { typedef T matrix_type[M][N]; @@ -574,10 +665,12 @@ namespace boost { namespace numeric { namespace ublas { // \todo { define correct wrapper } typedef matrix_reference closure_type; typedef const matrix_reference const_closure_type; + + // \todo { define appropriate iterators } }; template < class T, int N > - struct vector_traits < T[N] > { + struct vector_view_traits < T[N] > { typedef T vector_type[N]; typedef std::size_t size_type; @@ -586,12 +679,31 @@ namespace boost { namespace numeric { namespace ublas { typedef dense_tag storage_category; typedef T value_type; - typedef T *reference; typedef const T *const_reference; // \todo { define correct wrapper } - typedef vector_reference closure_type; typedef const vector_reference const_closure_type; + + typedef const_reference const_iterator; + + /// iterator pointing to the first element + static + const_iterator begin(const vector_type & v) { + return & (v[0]); + } + /// iterator pointing behind the last element + static + const_iterator end(const vector_type & v) { + return & (v[N]); + } + }; + + template < class T, int N > + struct mutable_vector_traits < T[N] > { + + typedef T *reference; + typedef vector_reference< T[N] > closure_type; + }; }}}