text: //
  //  Copyright (c) 2000-2007
  //  Joerg Walter, Mathias Koch, 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)
  //
  //  The authors gratefully acknowledge the support of
  //  GeNeSys mbH & Co. KG in producing this work.
  //
  
  #ifndef _BOOST_UBLAS_MATRIX_SPARSE_
  #define _BOOST_UBLAS_MATRIX_SPARSE_
  
  #include <boost/numeric/ublas/vector_sparse.hpp>
  #include <boost/numeric/ublas/matrix_expression.hpp>
  #include <boost/numeric/ublas/detail/matrix_assign.hpp>
  #if BOOST_UBLAS_TYPE_CHECK
  #include <boost/numeric/ublas/matrix.hpp>
  #endif
  
  // Iterators based on ideas of Jeremy Siek
  
  namespace boost { namespace numeric { namespace ublas {
  
  #ifdef BOOST_UBLAS_STRICT_MATRIX_SPARSE
  
      template<class M>
      class sparse_matrix_element:
         public container_reference<M> {
      public:
          typedef M matrix_type;
          typedef typename M::size_type size_type;
          typedef typename M::value_type value_type;
          typedef const value_type &const_reference;
          typedef value_type *pointer;
          typedef const value_type *const_pointer;
  
      private:
          // Proxied element operations
          void get_d () const {
              const_pointer p = (*this) ().find_element (i_, j_);
              if (p)
                  d_ = *p;
              else
                  d_ = value_type/*zero*/();
          }
  
          void set (const value_type &s) const {
              pointer p = (*this) ().find_element (i_, j_);
              if (!p)
                  (*this) ().insert_element (i_, j_, s);
              else
                  *p = s;
          }
          
      public:
          // Construction and destruction
          BOOST_UBLAS_INLINE
          sparse_matrix_element (matrix_type &m, size_type i, size_type j):
              container_reference<matrix_type> (m), i_ (i), j_ (j) {
          }
          BOOST_UBLAS_INLINE
          sparse_matrix_element (const sparse_matrix_element &p):
              container_reference<matrix_type> (p), i_ (p.i_), j_ (p.j_) {}
          BOOST_UBLAS_INLINE
          ~sparse_matrix_element () {
          }
  
          // Assignment
          BOOST_UBLAS_INLINE
          sparse_matrix_element &operator = (const sparse_matrix_element &p) {
              // Overide the implict copy assignment
              p.get_d ();
              set (p.d_);
              return *this;
          }
          template<class D>
          BOOST_UBLAS_INLINE
          sparse_matrix_element &operator = (const D &d) {
              set (d);
              return *this;
          }
          template<class D>
          BOOST_UBLAS_INLINE
          sparse_matrix_element &operator += (const D &d) {
              get_d ();
              d_ += d;
              set (d_);
              return *this;
          }
          template<class D>
          BOOST_UBLAS_INLINE
          sparse_matrix_element &operator -= (const D &d) {
              get_d ();
              d_ -= d;
              set (d_);
              return *this;
          }
          template<class D>
          BOOST_UBLAS_INLINE
          sparse_matrix_element &operator *= (const D &d) {
              get_d ();
              d_ *= d;
              set (d_);
              return *this;
          }
          template<class D>
          BOOST_UBLAS_INLINE
          sparse_matrix_element &operator /= (const D &d) {
              get_d ();
              d_ /= d;
              set (d_);
              return *this;
          }
  
          // Comparison
          template<class D>
          BOOST_UBLAS_INLINE
          bool operator == (const D &d) const {
              get_d ();
              return d_ == d;
          }
          template<class D>
          BOOST_UBLAS_INLINE
          bool operator != (const D &d) const {
              get_d ();
              return d_ != d;
          }
  
          // Conversion - weak link in proxy as d_ is not a perfect alias for the element
          BOOST_UBLAS_INLINE
          operator const_reference () const {
              get_d ();
              return d_;
          }
  
          // Conversion to reference - may be invalidated
          BOOST_UBLAS_INLINE
          value_type& ref () const {
              const pointer p = (*this) ().find_element (i_, j_);
              if (!p)
                  return (*this) ().insert_element (i_, j_, value_type/*zero*/());
              else
                  return *p;
          }
  
      private:
          size_type i_;
          size_type j_;
          mutable value_type d_;
      };
  
      /*
       * Generalise explicit reference access
       */
      namespace detail {
          template <class V>
          struct element_reference<sparse_matrix_element<V> > {
              typedef typename V::value_type& reference;
              static reference get_reference (const sparse_matrix_element<V>& sve)
              {
                  return sve.ref ();
              }
          };
      }
  
  
      template<class M>
      struct type_traits<sparse_matrix_element<M> > {
          typedef typename M::value_type element_type;
          typedef type_traits<sparse_matrix_element<M> > self_type;
          typedef typename type_traits<element_type>::value_type value_type;
          typedef typename type_traits<element_type>::const_reference const_reference;
          typedef sparse_matrix_element<M> reference;
          typedef typename type_traits<element_type>::real_type real_type;
          typedef typename type_traits<element_type>::precision_type precision_type;
  
          static const unsigned plus_complexity = type_traits<element_type>::plus_complexity;
          static const unsigned multiplies_complexity = type_traits<element_type>::multiplies_complexity;
  
          static
          BOOST_UBLAS_INLINE
          real_type real (const_reference t) {
              return type_traits<element_type>::real (t);
          }
          static
          BOOST_UBLAS_INLINE
          real_type imag (const_reference t) {
              return type_traits<element_type>::imag (t);
          }
          static
          BOOST_UBLAS_INLINE
          value_type conj (const_reference t) {
              return type_traits<element_type>::conj (t);
          }
  
          static
          BOOST_UBLAS_INLINE
          real_type type_abs (const_reference t) {
              return type_traits<element_type>::type_abs (t);
          }
          static
          BOOST_UBLAS_INLINE
          value_type type_sqrt (const_reference t) {
              return type_traits<element_type>::type_sqrt (t);
          }
  
          static
          BOOST_UBLAS_INLINE
          real_type norm_1 (const_reference t) {
              return type_traits<element_type>::norm_1 (t);
          }
          static
          BOOST_UBLAS_INLINE
          real_type norm_2 (const_reference t) {
              return type_traits<element_type>::norm_2 (t);
          }
          static
          BOOST_UBLAS_INLINE
          real_type norm_inf (const_reference t) {
              return type_traits<element_type>::norm_inf (t);
          }
  
          static
          BOOST_UBLAS_INLINE
          bool equals (const_reference t1, const_reference t2) {
              return type_traits<element_type>::equals (t1, t2);
          }
      };
  
      template<class M1, class T2>
      struct promote_traits<sparse_matrix_element<M1>, T2> {
          typedef typename promote_traits<typename sparse_matrix_element<M1>::value_type, T2>::promote_type promote_type;
      };
      template<class T1, class M2>
      struct promote_traits<T1, sparse_matrix_element<M2> > {
          typedef typename promote_traits<T1, typename sparse_matrix_element<M2>::value_type>::promote_type promote_type;
      };
      template<class M1, class M2>
      struct promote_traits<sparse_matrix_element<M1>, sparse_matrix_element<M2> > {
          typedef typename promote_traits<typename sparse_matrix_element<M1>::value_type,
                                          typename sparse_matrix_element<M2>::value_type>::promote_type promote_type;
      };
  
  #endif
  
     /** \brief Index map based sparse matrix of values of type \c T
      *
      * This class represents a matrix by using a \c key to value mapping. The default type is
      * \code template<class T, class L = row_major, class A =  map_std<std::size_t, T> > class mapped_matrix; \endcode
      * So, by default a STL map container is used to associate keys and values. The key is computed depending on 
      * the layout type \c L as \code key = layout_type::element(i, size1_, j, size2_); \endcode
      * which means \code key = (i*size2+j) \endcode for a row major matrix.
      * Limitations: The matrix size must not exceed \f$(size1*size2) < \f$ \code std::limits<std::size_t> \endcode. 
      * The \ref find1() and \ref find2() operations have a complexity of at least \f$\mathcal{O}(log(nnz))\f$, depending
      * on the efficiency of \c std::lower_bound on the key set of the map.
      * Orientation and storage can also be specified, otherwise a row major orientation is used. 
      * It is \b not required by the storage to initialize elements of the matrix. By default, the orientation is \c row_major. 
      *
      * \sa fwd.hpp, storage_sparse.hpp
      *
      * \tparam T the type of object stored in the matrix (like double, float, complex, etc...)
      * \tparam L the storage organization. It can be either \c row_major or \c column_major. By default it is \c row_major
      */
      template<class T, class L, class A>
      class mapped_matrix:
          public matrix_container<mapped_matrix<T, L, A> > {
  
          typedef T &true_reference;
          typedef T *pointer;
          typedef const T * const_pointer;
          typedef L layout_type;
          typedef mapped_matrix<T, L, A> self_type;
      public:
  #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
          using matrix_container<self_type>::operator ();
  #endif
          typedef typename A::size_type size_type;
          typedef typename A::difference_type difference_type;
          typedef T value_type;
          typedef A array_type;
          typedef const T &const_reference;
  #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
          typedef typename detail::map_traits<A, T>::reference reference;
  #else
          typedef sparse_matrix_element<self_type> reference;
  #endif
          typedef const matrix_reference<const self_type> const_closure_type;
          typedef matrix_reference<self_type> closure_type;
          typedef mapped_vector<T, A> vector_temporary_type;
          typedef self_type matrix_temporary_type;
          typedef sparse_tag storage_category;
          typedef typename L::orientation_category orientation_category;
  
          // Construction and destruction
          BOOST_UBLAS_INLINE
          mapped_matrix ():
              matrix_container<self_type> (),
              size1_ (0), size2_ (0), data_ () {}
          BOOST_UBLAS_INLINE
          mapped_matrix (size_type size1, size_type size2, size_type non_zeros = 0):
              matrix_container<self_type> (),
              size1_ (size1), size2_ (size2), data_ () {
              detail::map_reserve (data (), restrict_capacity (non_zeros));
          }
          BOOST_UBLAS_INLINE
          mapped_matrix (const mapped_matrix &m):
              matrix_container<self_type> (),
              size1_ (m.size1_), size2_ (m.size2_), data_ (m.data_) {}
          template<class AE>
          BOOST_UBLAS_INLINE
          mapped_matrix (const matrix_expression<AE> &ae, size_type non_zeros = 0):
              matrix_container<self_type> (),
              size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), data_ () {
              detail::map_reserve (data (), restrict_capacity (non_zeros));
              matrix_assign<scalar_assign> (*this, ae);
          }
  
          // Accessors
          BOOST_UBLAS_INLINE
          size_type size1 () const {
              return size1_;
          }
          BOOST_UBLAS_INLINE
          size_type size2 () const {
              return size2_;
          }
          BOOST_UBLAS_INLINE
          size_type nnz_capacity () const {
              return detail::map_capacity (data ());
          }
          BOOST_UBLAS_INLINE
          size_type nnz () const {
              return data (). size ();
          }
  
          // Storage accessors
          BOOST_UBLAS_INLINE
          const array_type &data () const {
              return data_;
          }
          BOOST_UBLAS_INLINE
          array_type &data () {
              return data_;
          }
  
          // Resizing
      private:
          BOOST_UBLAS_INLINE
          size_type restrict_capacity (size_type non_zeros) const {
              // Guarding against overflow - thanks to Alexei Novakov for the hint.
              // non_zeros = (std::min) (non_zeros, size1_ * size2_);
              if (size1_ > 0 && non_zeros / size1_ >= size2_)
                  non_zeros = size1_ * size2_;
              return non_zeros;
          }
      public:
          BOOST_UBLAS_INLINE
          void resize (size_type size1, size_type size2, bool preserve = true) {
              // FIXME preserve unimplemented
              BOOST_UBLAS_CHECK (!preserve, internal_logic ());
              size1_ = size1;
              size2_ = size2;
              data ().clear ();
          }
  
          // Reserving
          BOOST_UBLAS_INLINE
          void reserve (size_type non_zeros, bool preserve = true) {
              detail::map_reserve (data (), restrict_capacity (non_zeros));
          }
  
          // Element support
          BOOST_UBLAS_INLINE
          pointer find_element (size_type i, size_type j) {
              return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i, j));
          }
          BOOST_UBLAS_INLINE
          const_pointer find_element (size_type i, size_type j) const {
              const size_type element = layout_type::element (i, size1_, j, size2_);
              const_subiterator_type it (data ().find (element));
              if (it == data ().end ())
                  return 0;
              BOOST_UBLAS_CHECK ((*it).first == element, internal_logic ());   // broken map
              return &(*it).second;
          }
  
          // Element access
          BOOST_UBLAS_INLINE
          const_reference operator () (size_type i, size_type j) const {
              const size_type element = layout_type::element (i, size1_, j, size2_);
              const_subiterator_type it (data ().find (element));
              if (it == data ().end ())
                  return zero_;
              BOOST_UBLAS_CHECK ((*it).first == element, internal_logic ());   // broken map
              return (*it).second;
          }
          BOOST_UBLAS_INLINE
          reference operator () (size_type i, size_type j) {
  #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
              const size_type element = layout_type::element (i, size1_, j, size2_);
              std::pair<subiterator_type, bool> ii (data ().insert (typename array_type::value_type (element, value_type/*zero*/())));
              BOOST_UBLAS_CHECK ((ii.first)->first == element, internal_logic ());   // broken map
              return (ii.first)->second;
  #else
              return reference (*this, i, j);
  #endif
          }
  
          // Element assingment
          BOOST_UBLAS_INLINE
          true_reference insert_element (size_type i, size_type j, const_reference t) {
              BOOST_UBLAS_CHECK (!find_element (i, j), bad_index ());        // duplicate element
              const size_type element = layout_type::element (i, size1_, j, size2_);
              std::pair<subiterator_type, bool> ii (data ().insert (typename array_type::value_type (element, t)));
              BOOST_UBLAS_CHECK ((ii.first)->first == element, internal_logic ());   // broken map
              if (!ii.second)     // existing element
                  (ii.first)->second = t;
              return (ii.first)->second;
          }
          BOOST_UBLAS_INLINE
          void erase_element (size_type i, size_type j) {
              subiterator_type it = data ().find (layout_type::element (i, size1_, j, size2_));
              if (it == data ().end ())
                  return;
              data ().erase (it);
          }
          
          // Zeroing
          BOOST_UBLAS_INLINE
          void clear () {
              data ().clear ();
          }
  
          // Assignment
          BOOST_UBLAS_INLINE
          mapped_matrix &operator = (const mapped_matrix &m) {
              if (this != &m) {
                  size1_ = m.size1_;
                  size2_ = m.size2_;
                  data () = m.data ();
              }
              return *this;
          }
          template<class C>          // Container assignment without temporary
          BOOST_UBLAS_INLINE
          mapped_matrix &operator = (const matrix_container<C> &m) {
              resize (m ().size1 (), m ().size2 (), false);
              assign (m);
              return *this;
          }
          BOOST_UBLAS_INLINE
          mapped_matrix &assign_temporary (mapped_matrix &m) {
              swap (m);
              return *this;
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          mapped_matrix &operator = (const matrix_expression<AE> &ae) {
              self_type temporary (ae, detail::map_capacity (data ()));
              return assign_temporary (temporary);
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          mapped_matrix &assign (const matrix_expression<AE> &ae) {
              matrix_assign<scalar_assign> (*this, ae);
              return *this;
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          mapped_matrix& operator += (const matrix_expression<AE> &ae) {
              self_type temporary (*this + ae, detail::map_capacity (data ()));
              return assign_temporary (temporary);
          }
          template<class C>          // Container assignment without temporary
          BOOST_UBLAS_INLINE
          mapped_matrix &operator += (const matrix_container<C> &m) {
              plus_assign (m);
              return *this;
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          mapped_matrix &plus_assign (const matrix_expression<AE> &ae) {
              matrix_assign<scalar_plus_assign> (*this, ae);
              return *this;
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          mapped_matrix& operator -= (const matrix_expression<AE> &ae) {
              self_type temporary (*this - ae, detail::map_capacity (data ()));
              return assign_temporary (temporary);
          }
          template<class C>          // Container assignment without temporary
          BOOST_UBLAS_INLINE
          mapped_matrix &operator -= (const matrix_container<C> &m) {
              minus_assign (m);
              return *this;
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          mapped_matrix &minus_assign (const matrix_expression<AE> &ae) {
              matrix_assign<scalar_minus_assign> (*this, ae);
              return *this;
          }
          template<class AT>
          BOOST_UBLAS_INLINE
          mapped_matrix& operator *= (const AT &at) {
              matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
              return *this;
          }
          template<class AT>
          BOOST_UBLAS_INLINE
          mapped_matrix& operator /= (const AT &at) {
              matrix_assign_scalar<scalar_divides_assign> (*this, at);
              return *this;
          }
  
          // Swapping
          BOOST_UBLAS_INLINE
          void swap (mapped_matrix &m) {
              if (this != &m) {
                  std::swap (size1_, m.size1_);
                  std::swap (size2_, m.size2_);
                  data ().swap (m.data ());
              }
          }
          BOOST_UBLAS_INLINE
          friend void swap (mapped_matrix &m1, mapped_matrix &m2) {
              m1.swap (m2);
          }
  
          // Iterator types
      private:
          // Use storage iterator
          typedef typename A::const_iterator const_subiterator_type;
          typedef typename A::iterator subiterator_type;
  
          BOOST_UBLAS_INLINE
          true_reference at_element (size_type i, size_type j) {
              const size_type element = layout_type::element (i, size1_, j, size2_);
              subiterator_type it (data ().find (element));
              BOOST_UBLAS_CHECK (it != data ().end(), bad_index ());
              BOOST_UBLAS_CHECK ((*it).first == element, internal_logic ());   // broken map
              return it->second;
          }
  
      public:
          class const_iterator1;
          class iterator1;
          class const_iterator2;
          class iterator2;
          typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
          typedef reverse_iterator_base1<iterator1> reverse_iterator1;
          typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
          typedef reverse_iterator_base2<iterator2> reverse_iterator2;
  
          // Element lookup
          // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
          const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const {
              const_subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_)));
              const_subiterator_type it_end (data ().end ());
              size_type index1 = size_type (-1);
              size_type index2 = size_type (-1);
              while (rank == 1 && it != it_end) {
                  index1 = layout_type::index_i ((*it).first, size1_, size2_);
                  index2 = layout_type::index_j ((*it).first, size1_, size2_);
                  if (direction > 0) {
                      if ((index1 >= i && index2 == j) || (i >= size1_))
                          break;
                      ++ i;
                  } else /* if (direction < 0) */ {
                      if ((index1 <= i && index2 == j) || (i == 0))
                          break;
                      -- i;
                  }
                  it = data ().lower_bound (layout_type::address (i, size1_, j, size2_));
              }
              if (rank == 1 && index2 != j) {
                  if (direction > 0)
                      i = size1_;
                  else /* if (direction < 0) */
                      i = 0;
                  rank = 0;
              }
              return const_iterator1 (*this, rank, i, j, it);
          }
          // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
          iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) {
              subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_)));
              subiterator_type it_end (data ().end ());
              size_type index1 = size_type (-1);
              size_type index2 = size_type (-1);
              while (rank == 1 && it != it_end) {
                  index1 = layout_type::index_i ((*it).first, size1_, size2_);
                  index2 = layout_type::index_j ((*it).first, size1_, size2_);
                  if (direction > 0) {
                      if ((index1 >= i && index2 == j) || (i >= size1_))
                          break;
                      ++ i;
                  } else /* if (direction < 0) */ {
                      if ((index1 <= i && index2 == j) || (i == 0))
                          break;
                      -- i;
                  }
                  it = data ().lower_bound (layout_type::address (i, size1_, j, size2_));
              }
              if (rank == 1 && index2 != j) {
                  if (direction > 0)
                      i = size1_;
                  else /* if (direction < 0) */
                      i = 0;
                  rank = 0;
              }
              return iterator1 (*this, rank, i, j, it);
          }
          // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
          const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const {
              const_subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_)));
              const_subiterator_type it_end (data ().end ());
              size_type index1 = size_type (-1);
              size_type index2 = size_type (-1);
              while (rank == 1 && it != it_end) {
                  index1 = layout_type::index_i ((*it).first, size1_, size2_);
                  index2 = layout_type::index_j ((*it).first, size1_, size2_);
                  if (direction > 0) {
                      if ((index2 >= j && index1 == i) || (j >= size2_))
                          break;
                      ++ j;
                  } else /* if (direction < 0) */ {
                      if ((index2 <= j && index1 == i) || (j == 0))
                          break;
                      -- j;
                  }
                  it = data ().lower_bound (layout_type::address (i, size1_, j, size2_));
              }
              if (rank == 1 && index1 != i) {
                  if (direction > 0)
                      j = size2_;
                  else /* if (direction < 0) */
                      j = 0;
                  rank = 0;
              }
              return const_iterator2 (*this, rank, i, j, it);
          }
          // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
          iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) {
              subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_)));
              subiterator_type it_end (data ().end ());
              size_type index1 = size_type (-1);
              size_type index2 = size_type (-1);
              while (rank == 1 && it != it_end) {
                  index1 = layout_type::index_i ((*it).first, size1_, size2_);
                  index2 = layout_type::index_j ((*it).first, size1_, size2_);
                  if (direction > 0) {
                      if ((index2 >= j && index1 == i) || (j >= size2_))
                          break;
                      ++ j;
                  } else /* if (direction < 0) */ {
                      if ((index2 <= j && index1 == i) || (j == 0))
                          break;
                      -- j;
                  }
                  it = data ().lower_bound (layout_type::address (i, size1_, j, size2_));
              }
              if (rank == 1 && index1 != i) {
                  if (direction > 0)
                      j = size2_;
                  else /* if (direction < 0) */
                      j = 0;
                  rank = 0;
              }
              return iterator2 (*this, rank, i, j, it);
          }
  
  
          class const_iterator1:
              public container_const_reference<mapped_matrix>,
              public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
                                                 const_iterator1, value_type> {
          public:
              typedef typename mapped_matrix::value_type value_type;
              typedef typename mapped_matrix::difference_type difference_type;
              typedef typename mapped_matrix::const_reference reference;
              typedef const typename mapped_matrix::pointer pointer;
  
              typedef const_iterator2 dual_iterator_type;
              typedef const_reverse_iterator2 dual_reverse_iterator_type;
  
              // Construction and destruction
              BOOST_UBLAS_INLINE
              const_iterator1 ():
                  container_const_reference<self_type> (), rank_ (), i_ (), j_ (), it_ () {}
              BOOST_UBLAS_INLINE
              const_iterator1 (const self_type &m, int rank, size_type i, size_type j, const const_subiterator_type &it):
                  container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), it_ (it) {}
              BOOST_UBLAS_INLINE
              const_iterator1 (const iterator1 &it):
                  container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), it_ (it.it_) {}
  
              // Arithmetic
              BOOST_UBLAS_INLINE
              const_iterator1 &operator ++ () {
                  if (rank_ == 1 && layout_type::fast_i ())
                      ++ it_;
                  else
                      *this = (*this) ().find1 (rank_, index1 () + 1, j_, 1);
                  return *this;
              }
              BOOST_UBLAS_INLINE
              const_iterator1 &operator -- () {
                  if (rank_ == 1 && layout_type::fast_i ())
                      -- it_;
                  else
                      *this = (*this) ().find1 (rank_, index1 () - 1, j_, -1);
                  return *this;
              }
  
              // Dereference
              BOOST_UBLAS_INLINE
              const_reference operator * () const {
                  BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                  BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
                  if (rank_ == 1) {
                      return (*it_).second;
                  } else {
                      return (*this) () (i_, j_);
                  }
              }
  
  #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_iterator2 begin () const {
                  const self_type &m = (*this) ();
                  return m.find2 (1, index1 (), 0);
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_iterator2 end () const {
                  const self_type &m = (*this) ();
                  return m.find2 (1, index1 (), m.size2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_reverse_iterator2 rbegin () const {
                  return const_reverse_iterator2 (end ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_reverse_iterator2 rend () const {
                  return const_reverse_iterator2 (begin ());
              }
  #endif
  
              // Indices
              BOOST_UBLAS_INLINE
              size_type index1 () const {
                  BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
                  if (rank_ == 1) {
                      const self_type &m = (*this) ();
                      BOOST_UBLAS_CHECK (layout_type::index_i ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size1 (), bad_index ());
                      return layout_type::index_i ((*it_).first, m.size1 (), m.size2 ());
                  } else {
                      return i_;
                  }
              }
              BOOST_UBLAS_INLINE
              size_type index2 () const {
                  if (rank_ == 1) {
                      const self_type &m = (*this) ();
                      BOOST_UBLAS_CHECK (layout_type::index_j ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size2 (), bad_index ());
                      return layout_type::index_j ((*it_).first, m.size1 (), m.size2 ());
                  } else {
                      return j_;
                  }
              }
  
              // Assignment
              BOOST_UBLAS_INLINE
              const_iterator1 &operator = (const const_iterator1 &it) {
                  container_const_reference<self_type>::assign (&it ());
                  rank_ = it.rank_;
                  i_ = it.i_;
                  j_ = it.j_;
                  it_ = it.it_;
                  return *this;
              }
  
              // Comparison
              BOOST_UBLAS_INLINE
              bool operator == (const const_iterator1 &it) const {
                  BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
                  // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
                  if (rank_ == 1 || it.rank_ == 1) {
                      return it_ == it.it_;
                  } else {
                      return i_ == it.i_ && j_ == it.j_;
                  }
              }
  
          private:
              int rank_;
              size_type i_;
              size_type j_;
              const_subiterator_type it_;
          };
  
          BOOST_UBLAS_INLINE
          const_iterator1 begin1 () const {
              return find1 (0, 0, 0);
          }
          BOOST_UBLAS_INLINE
          const_iterator1 end1 () const {
              return find1 (0, size1_, 0);
          }
  
          class iterator1:
              public container_reference<mapped_matrix>,
              public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
                                                 iterator1, value_type> {
          public:
              typedef typename mapped_matrix::value_type value_type;
              typedef typename mapped_matrix::difference_type difference_type;
              typedef typename mapped_matrix::true_reference reference;
              typedef typename mapped_matrix::pointer pointer;
  
              typedef iterator2 dual_iterator_type;
              typedef reverse_iterator2 dual_reverse_iterator_type;
  
              // Construction and destruction
              BOOST_UBLAS_INLINE
              iterator1 ():
                  container_reference<self_type> (), rank_ (), i_ (), j_ (), it_ () {}
              BOOST_UBLAS_INLINE
              iterator1 (self_type &m, int rank, size_type i, size_type j, const subiterator_type &it):
                  container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), it_ (it) {}
  
              // Arithmetic
              BOOST_UBLAS_INLINE
              iterator1 &operator ++ () {
                  if (rank_ == 1 && layout_type::fast_i ())
                      ++ it_;
                  else
                      *this = (*this) ().find1 (rank_, index1 () + 1, j_, 1);
                  return *this;
              }
              BOOST_UBLAS_INLINE
              iterator1 &operator -- () {
                  if (rank_ == 1 && layout_type::fast_i ())
                      -- it_;
                  else
                      *this = (*this) ().find1 (rank_, index1 () - 1, j_, -1);
                  return *this;
              }
  
              // Dereference
              BOOST_UBLAS_INLINE
              reference operator * () const {
                  BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                  BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
                  if (rank_ == 1) {
                      return (*it_).second;
                  } else {
                      return (*this) ().at_element (i_, j_);
                  }
              }
  
  #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              iterator2 begin () const {
                  self_type &m = (*this) ();
                  return m.find2 (1, index1 (), 0);
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              iterator2 end () const {
                  self_type &m = (*this) ();
                  return m.find2 (1, index1 (), m.size2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              reverse_iterator2 rbegin () const {
                  return reverse_iterator2 (end ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              reverse_iterator2 rend () const {
                  return reverse_iterator2 (begin ());
              }
  #endif
  
              // Indices
              BOOST_UBLAS_INLINE
              size_type index1 () const {
                  BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
                  if (rank_ == 1) {
                      const self_type &m = (*this) ();
                      BOOST_UBLAS_CHECK (layout_type::index_i ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size1 (), bad_index ());
                      return layout_type::index_i ((*it_).first, m.size1 (), m.size2 ());
                  } else {
                      return i_;
                  }
              }
              BOOST_UBLAS_INLINE
              size_type index2 () const {
                  if (rank_ == 1) {
                      const self_type &m = (*this) ();
                      BOOST_UBLAS_CHECK (layout_type::index_j ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size2 (), bad_index ());
                      return layout_type::index_j ((*it_).first, m.size1 (), m.size2 ());
                  } else {
                      return j_;
                  }
              }
  
              // Assignment
              BOOST_UBLAS_INLINE
              iterator1 &operator = (const iterator1 &it) {
                  container_reference<self_type>::assign (&it ());
                  rank_ = it.rank_;
                  i_ = it.i_;
                  j_ = it.j_;
                  it_ = it.it_;
                  return *this;
              }
  
              // Comparison
              BOOST_UBLAS_INLINE
              bool operator == (const iterator1 &it) const {
                  BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
                  // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
                  if (rank_ == 1 || it.rank_ == 1) {
                      return it_ == it.it_;
                  } else {
                      return i_ == it.i_ && j_ == it.j_;
                  }
              }
  
          private:
              int rank_;
              size_type i_;
              size_type j_;
              subiterator_type it_;
  
              friend class const_iterator1;
          };
  
          BOOST_UBLAS_INLINE
          iterator1 begin1 () {
              return find1 (0, 0, 0);
          }
          BOOST_UBLAS_INLINE
          iterator1 end1 () {
              return find1 (0, size1_, 0);
          }
  
          class const_iterator2:
              public container_const_reference<mapped_matrix>,
              public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
                                                 const_iterator2, value_type> {
          public:
              typedef typename mapped_matrix::value_type value_type;
              typedef typename mapped_matrix::difference_type difference_type;
              typedef typename mapped_matrix::const_reference reference;
              typedef const typename mapped_matrix::pointer pointer;
  
              typedef const_iterator1 dual_iterator_type;
              typedef const_reverse_iterator1 dual_reverse_iterator_type;
  
              // Construction and destruction
              BOOST_UBLAS_INLINE
              const_iterator2 ():
                  container_const_reference<self_type> (), rank_ (), i_ (), j_ (), it_ () {}
              BOOST_UBLAS_INLINE
              const_iterator2 (const self_type &m, int rank, size_type i, size_type j, const const_subiterator_type &it):
                  container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), it_ (it) {}
              BOOST_UBLAS_INLINE
              const_iterator2 (const iterator2 &it):
                  container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), it_ (it.it_) {}
  
              // Arithmetic
              BOOST_UBLAS_INLINE
              const_iterator2 &operator ++ () {
                  if (rank_ == 1 && layout_type::fast_j ())
                      ++ it_;
                  else
                      *this = (*this) ().find2 (rank_, i_, index2 () + 1, 1);
                  return *this;
              }
              BOOST_UBLAS_INLINE
              const_iterator2 &operator -- () {
                  if (rank_ == 1 && layout_type::fast_j ())
                      -- it_;
                  else
                      *this = (*this) ().find2 (rank_, i_, index2 () - 1, -1);
                  return *this;
              }
  
              // Dereference
              BOOST_UBLAS_INLINE
              const_reference operator * () const {
                  BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                  BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
                  if (rank_ == 1) {
                      return (*it_).second;
                  } else {
                      return (*this) () (i_, j_);
                  }
              }
  
  #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_iterator1 begin () const {
                  const self_type &m = (*this) ();
                  return m.find1 (1, 0, index2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_iterator1 end () const {
                  const self_type &m = (*this) ();
                  return m.find1 (1, m.size1 (), index2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_reverse_iterator1 rbegin () const {
                  return const_reverse_iterator1 (end ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_reverse_iterator1 rend () const {
                  return const_reverse_iterator1 (begin ());
              }
  #endif
  
              // Indices
              BOOST_UBLAS_INLINE
              size_type index1 () const {
                  if (rank_ == 1) {
                      const self_type &m = (*this) ();
                      BOOST_UBLAS_CHECK (layout_type::index_i ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size1 (), bad_index ());
                      return layout_type::index_i ((*it_).first, m.size1 (), m.size2 ());
                  } else {
                      return i_;
                  }
              }
              BOOST_UBLAS_INLINE
              size_type index2 () const {
                  BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
                  if (rank_ == 1) {
                      const self_type &m = (*this) ();
                      BOOST_UBLAS_CHECK (layout_type::index_j ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size2 (), bad_index ());
                      return layout_type::index_j ((*it_).first, m.size1 (), m.size2 ());
                  } else {
                      return j_;
                  }
              }
  
              // Assignment
              BOOST_UBLAS_INLINE
              const_iterator2 &operator = (const const_iterator2 &it) {
                  container_const_reference<self_type>::assign (&it ());
                  rank_ = it.rank_;
                  i_ = it.i_;
                  j_ = it.j_;
                  it_ = it.it_;
                  return *this;
              }
  
              // Comparison
              BOOST_UBLAS_INLINE
              bool operator == (const const_iterator2 &it) const {
                  BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
                  // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
                  if (rank_ == 1 || it.rank_ == 1) {
                      return it_ == it.it_;
                  } else {
                      return i_ == it.i_ && j_ == it.j_;
                  }
              }
  
          private:
              int rank_;
              size_type i_;
              size_type j_;
              const_subiterator_type it_;
          };
  
          BOOST_UBLAS_INLINE
          const_iterator2 begin2 () const {
              return find2 (0, 0, 0);
          }
          BOOST_UBLAS_INLINE
          const_iterator2 end2 () const {
              return find2 (0, 0, size2_);
          }
  
          class iterator2:
              public container_reference<mapped_matrix>,
              public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
                                                 iterator2, value_type> {
          public:
              typedef typename mapped_matrix::value_type value_type;
              typedef typename mapped_matrix::difference_type difference_type;
              typedef typename mapped_matrix::true_reference reference;
              typedef typename mapped_matrix::pointer pointer;
  
              typedef iterator1 dual_iterator_type;
              typedef reverse_iterator1 dual_reverse_iterator_type;
  
              // Construction and destruction
              BOOST_UBLAS_INLINE
              iterator2 ():
                  container_reference<self_type> (), rank_ (), i_ (), j_ (), it_ () {}
              BOOST_UBLAS_INLINE
              iterator2 (self_type &m, int rank, size_type i, size_type j, const subiterator_type &it):
                  container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), it_ (it) {}
  
              // Arithmetic
              BOOST_UBLAS_INLINE
              iterator2 &operator ++ () {
                  if (rank_ == 1 && layout_type::fast_j ())
                      ++ it_;
                  else
                      *this = (*this) ().find2 (rank_, i_, index2 () + 1, 1);
                  return *this;
              }
              BOOST_UBLAS_INLINE
              iterator2 &operator -- () {
                  if (rank_ == 1 && layout_type::fast_j ())
                      -- it_;
                  else
                      *this = (*this) ().find2 (rank_, i_, index2 () - 1, -1);
                  return *this;
              }
  
              // Dereference
              BOOST_UBLAS_INLINE
              reference operator * () const {
                  BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                  BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
                  if (rank_ == 1) {
                      return (*it_).second;
                  } else {
                      return (*this) ().at_element (i_, j_);
                  }
              }
  
  #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              iterator1 begin () const {
                  self_type &m = (*this) ();
                  return m.find1 (1, 0, index2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              iterator1 end () const {
                  self_type &m = (*this) ();
                  return m.find1 (1, m.size1 (), index2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              reverse_iterator1 rbegin () const {
                  return reverse_iterator1 (end ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              reverse_iterator1 rend () const {
                  return reverse_iterator1 (begin ());
              }
  #endif
  
              // Indices
              BOOST_UBLAS_INLINE
              size_type index1 () const {
                  if (rank_ == 1) {
                      const self_type &m = (*this) ();
                      BOOST_UBLAS_CHECK (layout_type::index_i ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size1 (), bad_index ());
                      return layout_type::index_i ((*it_).first, m.size1 (), m.size2 ());
                  } else {
                      return i_;
                  }
              }
              BOOST_UBLAS_INLINE
              size_type index2 () const {
                  BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
                  if (rank_ == 1) {
                      const self_type &m = (*this) ();
                      BOOST_UBLAS_CHECK (layout_type::index_j ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size2 (), bad_index ());
                      return layout_type::index_j ((*it_).first, m.size1 (), m.size2 ());
                  } else {
                      return j_;
                  }
              }
  
              // Assignment
              BOOST_UBLAS_INLINE
              iterator2 &operator = (const iterator2 &it) {
                  container_reference<self_type>::assign (&it ());
                  rank_ = it.rank_;
                  i_ = it.i_;
                  j_ = it.j_;
                  it_ = it.it_;
                  return *this;
              }
  
              // Comparison
              BOOST_UBLAS_INLINE
              bool operator == (const iterator2 &it) const {
                  BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
                  // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
                  if (rank_ == 1 || it.rank_ == 1) {
                      return it_ == it.it_;
                  } else {
                      return i_ == it.i_ && j_ == it.j_;
                  }
              }
  
          private:
              int rank_;
              size_type i_;
              size_type j_;
              subiterator_type it_;
  
              friend class const_iterator2;
          };
  
          BOOST_UBLAS_INLINE
          iterator2 begin2 () {
              return find2 (0, 0, 0);
          }
          BOOST_UBLAS_INLINE
          iterator2 end2 () {
              return find2 (0, 0, size2_);
          }
  
          // Reverse iterators
  
          BOOST_UBLAS_INLINE
          const_reverse_iterator1 rbegin1 () const {
              return const_reverse_iterator1 (end1 ());
          }
          BOOST_UBLAS_INLINE
          const_reverse_iterator1 rend1 () const {
              return const_reverse_iterator1 (begin1 ());
          }
  
          BOOST_UBLAS_INLINE
          reverse_iterator1 rbegin1 () {
              return reverse_iterator1 (end1 ());
          }
          BOOST_UBLAS_INLINE
          reverse_iterator1 rend1 () {
              return reverse_iterator1 (begin1 ());
          }
  
          BOOST_UBLAS_INLINE
          const_reverse_iterator2 rbegin2 () const {
              return const_reverse_iterator2 (end2 ());
          }
          BOOST_UBLAS_INLINE
          const_reverse_iterator2 rend2 () const {
              return const_reverse_iterator2 (begin2 ());
          }
  
          BOOST_UBLAS_INLINE
          reverse_iterator2 rbegin2 () {
              return reverse_iterator2 (end2 ());
          }
          BOOST_UBLAS_INLINE
          reverse_iterator2 rend2 () {
              return reverse_iterator2 (begin2 ());
          }
  
           // Serialization
          template<class Archive>
          void serialize(Archive & ar, const unsigned int /* file_version */){
              serialization::collection_size_type s1 (size1_);
              serialization::collection_size_type s2 (size2_);
              ar & serialization::make_nvp("size1",s1);
              ar & serialization::make_nvp("size2",s2);
              if (Archive::is_loading::value) {
                  size1_ = s1;
                  size2_ = s2;
              }
              ar & serialization::make_nvp("data", data_);
          }
  
      private:
          size_type size1_;
          size_type size2_;
          array_type data_;
          static const value_type zero_;
      };
  
      template<class T, class L, class A>
      const typename mapped_matrix<T, L, A>::value_type mapped_matrix<T, L, A>::zero_ = value_type/*zero*/();
  
  
      // Vector index map based sparse matrix class
      template<class T, class L, class A>
      class mapped_vector_of_mapped_vector:
          public matrix_container<mapped_vector_of_mapped_vector<T, L, A> > {
  
          typedef T &true_reference;
          typedef T *pointer;
          typedef const T *const_pointer;
          typedef A array_type;
          typedef const A const_array_type;
          typedef L layout_type;
          typedef mapped_vector_of_mapped_vector<T, L, A> self_type;
      public:
  #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
          using matrix_container<self_type>::operator ();
  #endif
          typedef typename A::size_type size_type;
          typedef typename A::difference_type difference_type;
          typedef T value_type;
          typedef const T &const_reference;
  #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
          typedef typename detail::map_traits<typename A::data_value_type, T>::reference reference;
  #else
          typedef sparse_matrix_element<self_type> reference;
  #endif
          typedef const matrix_reference<const self_type> const_closure_type;
          typedef matrix_reference<self_type> closure_type;
          typedef mapped_vector<T, typename A::value_type> vector_temporary_type;
          typedef self_type matrix_temporary_type;
          typedef typename A::value_type::second_type vector_data_value_type;
          typedef sparse_tag storage_category;
          typedef typename L::orientation_category orientation_category;
  
          // Construction and destruction
          BOOST_UBLAS_INLINE
          mapped_vector_of_mapped_vector ():
              matrix_container<self_type> (),
              size1_ (0), size2_ (0), data_ () {
              data_ [layout_type::size_M (size1_, size2_)] = vector_data_value_type ();
          }
          BOOST_UBLAS_INLINE
          mapped_vector_of_mapped_vector (size_type size1, size_type size2, size_type non_zeros = 0):
              matrix_container<self_type> (),
              size1_ (size1), size2_ (size2), data_ () {
              data_ [layout_type::size_M (size1_, size2_)] = vector_data_value_type ();
          }
          BOOST_UBLAS_INLINE
          mapped_vector_of_mapped_vector (const mapped_vector_of_mapped_vector &m):
              matrix_container<self_type> (),
              size1_ (m.size1_), size2_ (m.size2_), data_ (m.data_) {}
          template<class AE>
          BOOST_UBLAS_INLINE
          mapped_vector_of_mapped_vector (const matrix_expression<AE> &ae, size_type non_zeros = 0):
              matrix_container<self_type> (),
              size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), data_ () {
              data_ [layout_type::size_M (size1_, size2_)] = vector_data_value_type ();
              matrix_assign<scalar_assign> (*this, ae);
          }
  
          // Accessors
          BOOST_UBLAS_INLINE
          size_type size1 () const {
              return size1_;
          }
          BOOST_UBLAS_INLINE
          size_type size2 () const {
              return size2_;
          }
          BOOST_UBLAS_INLINE
          size_type nnz_capacity () const {
              size_type non_zeros = 0;
              for (vector_const_subiterator_type itv = data_ ().begin (); itv != data_ ().end (); ++ itv)
                  non_zeros += detail::map_capacity (*itv);
              return non_zeros;
          }
          BOOST_UBLAS_INLINE
          size_type nnz () const {
              size_type filled = 0;
              for (vector_const_subiterator_type itv = data_ ().begin (); itv != data_ ().end (); ++ itv)
                  filled += (*itv).size ();
              return filled;
          }
  
          // Storage accessors
          BOOST_UBLAS_INLINE
          const_array_type &data () const {
              return data_;
          }
          BOOST_UBLAS_INLINE
          array_type &data () {
              return data_;
          }
  
          // Resizing
          BOOST_UBLAS_INLINE
          void resize (size_type size1, size_type size2, bool preserve = true) {
              // FIXME preserve unimplemented
              BOOST_UBLAS_CHECK (!preserve, internal_logic ());
              size1_ = size1;
              size2_ = size2;
              data ().clear ();
              data () [layout_type::size_M (size1_, size2_)] = vector_data_value_type ();
          }
  
          // Element support
          BOOST_UBLAS_INLINE
          pointer find_element (size_type i, size_type j) {
              return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i, j));
          }
          BOOST_UBLAS_INLINE
          const_pointer find_element (size_type i, size_type j) const {
              const size_type element1 = layout_type::index_M (i, j);
              const size_type element2 = layout_type::index_m (i, j);
              vector_const_subiterator_type itv (data ().find (element1));
              if (itv == data ().end ())
                  return 0;
              BOOST_UBLAS_CHECK ((*itv).first == element1, internal_logic ());   // broken map
              const_subiterator_type it ((*itv).second.find (element2));
              if (it == (*itv).second.end ())
                  return 0;
              BOOST_UBLAS_CHECK ((*it).first == element2, internal_logic ());   // broken map
              return &(*it).second;
          }
  
          // Element access
          BOOST_UBLAS_INLINE
          const_reference operator () (size_type i, size_type j) const {
              const size_type element1 = layout_type::index_M (i, j);
              const size_type element2 = layout_type::index_m (i, j);
              vector_const_subiterator_type itv (data ().find (element1));
              if (itv == data ().end ())
                  return zero_;
              BOOST_UBLAS_CHECK ((*itv).first == element1, internal_logic ());   // broken map
              const_subiterator_type it ((*itv).second.find (element2));
              if (it == (*itv).second.end ())
                  return zero_;
              BOOST_UBLAS_CHECK ((*itv).first == element1, internal_logic ());   // broken map
              return (*it).second;
          }
          BOOST_UBLAS_INLINE
          reference operator () (size_type i, size_type j) {
  #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
              const size_type element1 = layout_type::index_M (i, j);
              const size_type element2 = layout_type::index_m (i, j);
              vector_data_value_type& vd (data () [element1]);
              std::pair<subiterator_type, bool> ii (vd.insert (typename array_type::value_type::second_type::value_type (element2, value_type/*zero*/())));
              BOOST_UBLAS_CHECK ((ii.first)->first == element2, internal_logic ());   // broken map
              return (ii.first)->second;
  #else
              return reference (*this, i, j);
  #endif
          }
  
          // Element assignment
          BOOST_UBLAS_INLINE
          true_reference insert_element (size_type i, size_type j, const_reference t) {
              BOOST_UBLAS_CHECK (!find_element (i, j), bad_index ());        // duplicate element
              const size_type element1 = layout_type::index_M (i, j);
              const size_type element2 = layout_type::index_m (i, j);
  
              vector_data_value_type& vd (data () [element1]);
              std::pair<subiterator_type, bool> ii (vd.insert (typename vector_data_value_type::value_type (element2, t)));
              BOOST_UBLAS_CHECK ((ii.first)->first == element2, internal_logic ());   // broken map
              if (!ii.second)     // existing element
                  (ii.first)->second = t;
              return (ii.first)->second;
          }
          BOOST_UBLAS_INLINE
          void erase_element (size_type i, size_type j) {
              vector_subiterator_type itv (data ().find (layout_type::index_M (i, j)));
              if (itv == data ().end ())
                  return;
              subiterator_type it ((*itv).second.find (layout_type::index_m (i, j)));
              if (it == (*itv).second.end ())
                  return;
              (*itv).second.erase (it);
          }
          
          // Zeroing
          BOOST_UBLAS_INLINE
          void clear () {
              data ().clear ();
              data_ [layout_type::size_M (size1_, size2_)] = vector_data_value_type ();
          }
  
          // Assignment
          BOOST_UBLAS_INLINE
          mapped_vector_of_mapped_vector &operator = (const mapped_vector_of_mapped_vector &m) {
              if (this != &m) {
                  size1_ = m.size1_;
                  size2_ = m.size2_;
                  data () = m.data ();
              }
              return *this;
          }
          template<class C>          // Container assignment without temporary
          BOOST_UBLAS_INLINE
          mapped_vector_of_mapped_vector &operator = (const matrix_container<C> &m) {
              resize (m ().size1 (), m ().size2 ());
              assign (m);
              return *this;
          }
          BOOST_UBLAS_INLINE
          mapped_vector_of_mapped_vector &assign_temporary (mapped_vector_of_mapped_vector &m) {
              swap (m);
              return *this;
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          mapped_vector_of_mapped_vector &operator = (const matrix_expression<AE> &ae) {
              self_type temporary (ae);
              return assign_temporary (temporary);
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          mapped_vector_of_mapped_vector &assign (const matrix_expression<AE> &ae) {
              matrix_assign<scalar_assign> (*this, ae);
              return *this;
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          mapped_vector_of_mapped_vector& operator += (const matrix_expression<AE> &ae) {
              self_type temporary (*this + ae);
              return assign_temporary (temporary);
          }
          template<class C>          // Container assignment without temporary
          BOOST_UBLAS_INLINE
          mapped_vector_of_mapped_vector &operator += (const matrix_container<C> &m) {
              plus_assign (m);
              return *this;
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          mapped_vector_of_mapped_vector &plus_assign (const matrix_expression<AE> &ae) {
              matrix_assign<scalar_plus_assign> (*this, ae);
              return *this;
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          mapped_vector_of_mapped_vector& operator -= (const matrix_expression<AE> &ae) {
              self_type temporary (*this - ae);
              return assign_temporary (temporary);
          }
          template<class C>          // Container assignment without temporary
          BOOST_UBLAS_INLINE
          mapped_vector_of_mapped_vector &operator -= (const matrix_container<C> &m) {
              minus_assign (m);
              return *this;
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          mapped_vector_of_mapped_vector &minus_assign (const matrix_expression<AE> &ae) {
              matrix_assign<scalar_minus_assign> (*this, ae);
              return *this;
          }
          template<class AT>
          BOOST_UBLAS_INLINE
          mapped_vector_of_mapped_vector& operator *= (const AT &at) {
              matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
              return *this;
          }
          template<class AT>
          BOOST_UBLAS_INLINE
          mapped_vector_of_mapped_vector& operator /= (const AT &at) {
              matrix_assign_scalar<scalar_divides_assign> (*this, at);
              return *this;
          }
  
          // Swapping
          BOOST_UBLAS_INLINE
          void swap (mapped_vector_of_mapped_vector &m) {
              if (this != &m) {
                  std::swap (size1_, m.size1_);
                  std::swap (size2_, m.size2_);
                  data ().swap (m.data ());
              }
          }
          BOOST_UBLAS_INLINE
          friend void swap (mapped_vector_of_mapped_vector &m1, mapped_vector_of_mapped_vector &m2) {
              m1.swap (m2);
          }
  
          // Iterator types
      private:
          // Use storage iterators
          typedef typename A::const_iterator vector_const_subiterator_type;
          typedef typename A::iterator vector_subiterator_type;
          typedef typename A::value_type::second_type::const_iterator const_subiterator_type;
          typedef typename A::value_type::second_type::iterator subiterator_type;
  
          BOOST_UBLAS_INLINE
          true_reference at_element (size_type i, size_type j) {
              const size_type element1 = layout_type::index_M (i, j);
              const size_type element2 = layout_type::index_m (i, j);
              vector_subiterator_type itv (data ().find (element1));
              BOOST_UBLAS_CHECK (itv != data ().end(), bad_index ());
              BOOST_UBLAS_CHECK ((*itv).first == element1, internal_logic ());   // broken map
              subiterator_type it ((*itv).second.find (element2));
              BOOST_UBLAS_CHECK (it != (*itv).second.end (), bad_index ());
              BOOST_UBLAS_CHECK ((*it).first == element2, internal_logic ());    // broken map
              
              return it->second;
          }
  
      public:
          class const_iterator1;
          class iterator1;
          class const_iterator2;
          class iterator2;
          typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
          typedef reverse_iterator_base1<iterator1> reverse_iterator1;
          typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
          typedef reverse_iterator_base2<iterator2> reverse_iterator2;
  
          // Element lookup
          // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
          const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const {
              BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ());
              for (;;) {
                  vector_const_subiterator_type itv (data ().lower_bound (layout_type::index_M (i, j)));
                  vector_const_subiterator_type itv_end (data ().end ());
                  if (itv == itv_end)
                      return const_iterator1 (*this, rank, i, j, itv_end, (*(-- itv)).second.end ());
  
                  const_subiterator_type it ((*itv).second.lower_bound (layout_type::index_m (i, j)));
                  const_subiterator_type it_end ((*itv).second.end ());
                  if (rank == 0) {
                      // advance to the first available major index
                      size_type M = itv->first;
                      size_type m;
                      if (it != it_end) { 
                          m = it->first; 
                      } else {
                          m = layout_type::size_m(size1_, size2_);
                      }
                      size_type first_i = layout_type::index_M(M,m);
                      return const_iterator1 (*this, rank, first_i, j, itv, it);
                  }
                  if (it != it_end && (*it).first == layout_type::index_m (i, j))
                      return const_iterator1 (*this, rank, i, j, itv, it);
                  if (direction > 0) {
                      if (layout_type::fast_i ()) {
                          if (it == it_end)
                              return const_iterator1 (*this, rank, i, j, itv, it);
                          i = (*it).first;
                      } else {
                          if (i >= size1_)
                              return const_iterator1 (*this, rank, i, j, itv, it);
                          ++ i;
                      }
                  } else /* if (direction < 0)  */ {
                      if (layout_type::fast_i ()) {
                          if (it == (*itv).second.begin ())
                              return const_iterator1 (*this, rank, i, j, itv, it);
                          -- it;
                          i = (*it).first;
                      } else {
                          if (i == 0)
                              return const_iterator1 (*this, rank, i, j, itv, it);
                          -- i;
                      }
                  }
              }
          }
          // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
          iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) {
              BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ());
              for (;;) {
                  vector_subiterator_type itv (data ().lower_bound (layout_type::index_M (i, j)));
                  vector_subiterator_type itv_end (data ().end ());
                  if (itv == itv_end)
                      return iterator1 (*this, rank, i, j, itv_end, (*(-- itv)).second.end ());
  
                  subiterator_type it ((*itv).second.lower_bound (layout_type::index_m (i, j)));
                  subiterator_type it_end ((*itv).second.end ());
                  if (rank == 0) {
                      // advance to the first available major index
                      size_type M = itv->first;
                      size_type m;
                      if (it != it_end) { 
                          m = it->first; 
                      } else {
                          m = layout_type::size_m(size1_, size2_);
                      }
                      size_type first_i = layout_type::index_M(M,m);
                      return iterator1 (*this, rank, first_i, j, itv, it);
                  }
                  if (it != it_end && (*it).first == layout_type::index_m (i, j))
                      return iterator1 (*this, rank, i, j, itv, it);
                  if (direction > 0) {
                      if (layout_type::fast_i ()) {
                          if (it == it_end)
                              return iterator1 (*this, rank, i, j, itv, it);
                          i = (*it).first;
                      } else {
                          if (i >= size1_)
                              return iterator1 (*this, rank, i, j, itv, it);
                          ++ i;
                      }
                  } else /* if (direction < 0)  */ {
                      if (layout_type::fast_i ()) {
                          if (it == (*itv).second.begin ())
                              return iterator1 (*this, rank, i, j, itv, it);
                          -- it;
                          i = (*it).first;
                      } else {
                          if (i == 0)
                              return iterator1 (*this, rank, i, j, itv, it);
                          -- i;
                      }
                  }
              }
          }
          // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
          const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const {
              BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ());
              for (;;) {
                  vector_const_subiterator_type itv (data ().lower_bound (layout_type::index_M (i, j)));
                  vector_const_subiterator_type itv_end (data ().end ());
                  if (itv == itv_end)
                      return const_iterator2 (*this, rank, i, j, itv_end, (*(-- itv)).second.end ());
  
                  const_subiterator_type it ((*itv).second.lower_bound (layout_type::index_m (i, j)));
                  const_subiterator_type it_end ((*itv).second.end ());
                  if (rank == 0) {
                      // advance to the first available major index
                      size_type M = itv->first;
                      size_type m;
                      if (it != it_end) { 
                          m = it->first; 
                      } else {
                          m = layout_type::size_m(size1_, size2_);
                      }
                      size_type first_j = layout_type::index_m(M,m);
                      return const_iterator2 (*this, rank, i, first_j, itv, it);
                  }
                  if (it != it_end && (*it).first == layout_type::index_m (i, j))
                      return const_iterator2 (*this, rank, i, j, itv, it);
                  if (direction > 0) {
                      if (layout_type::fast_j ()) {
                          if (it == it_end)
                              return const_iterator2 (*this, rank, i, j, itv, it);
                          j = (*it).first;
                      } else {
                          if (j >= size2_)
                              return const_iterator2 (*this, rank, i, j, itv, it);
                          ++ j;
                      }
                  } else /* if (direction < 0)  */ {
                      if (layout_type::fast_j ()) {
                          if (it == (*itv).second.begin ())
                              return const_iterator2 (*this, rank, i, j, itv, it);
                          -- it;
                          j = (*it).first;
                      } else {
                          if (j == 0)
                              return const_iterator2 (*this, rank, i, j, itv, it);
                          -- j;
                      }
                  }
              }
          }
          // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
          iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) {
              BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ());
              for (;;) {
                  vector_subiterator_type itv (data ().lower_bound (layout_type::index_M (i, j)));
                  vector_subiterator_type itv_end (data ().end ());
                  if (itv == itv_end)
                      return iterator2 (*this, rank, i, j, itv_end, (*(-- itv)).second.end ());
  
                  subiterator_type it ((*itv).second.lower_bound (layout_type::index_m (i, j)));
                  subiterator_type it_end ((*itv).second.end ());
                  if (rank == 0) {
                      // advance to the first available major index
                      size_type M = itv->first;
                      size_type m;
                      if (it != it_end) { 
                          m = it->first; 
                      } else {
                          m = layout_type::size_m(size1_, size2_);
                      }
                      size_type first_j = layout_type::index_m(M,m);
                      return iterator2 (*this, rank, i, first_j, itv, it);
                  }
                  if (it != it_end && (*it).first == layout_type::index_m (i, j))
                      return iterator2 (*this, rank, i, j, itv, it);
                  if (direction > 0) {
                      if (layout_type::fast_j ()) {
                          if (it == it_end)
                              return iterator2 (*this, rank, i, j, itv, it);
                          j = (*it).first;
                      } else {
                          if (j >= size2_)
                              return iterator2 (*this, rank, i, j, itv, it);
                          ++ j;
                      }
                  } else /* if (direction < 0)  */ {
                      if (layout_type::fast_j ()) {
                          if (it == (*itv).second.begin ())
                              return iterator2 (*this, rank, i, j, itv, it);
                          -- it;
                          j = (*it).first;
                      } else {
                          if (j == 0)
                              return iterator2 (*this, rank, i, j, itv, it);
                          -- j;
                      }
                  }
              }
          }
  
          class const_iterator1:
              public container_const_reference<mapped_vector_of_mapped_vector>,
              public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
                                                 const_iterator1, value_type> {
          public:
              typedef typename mapped_vector_of_mapped_vector::value_type value_type;
              typedef typename mapped_vector_of_mapped_vector::difference_type difference_type;
              typedef typename mapped_vector_of_mapped_vector::const_reference reference;
              typedef const typename mapped_vector_of_mapped_vector::pointer pointer;
  
              typedef const_iterator2 dual_iterator_type;
              typedef const_reverse_iterator2 dual_reverse_iterator_type;
  
              // Construction and destruction
              BOOST_UBLAS_INLINE
              const_iterator1 ():
                  container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
              BOOST_UBLAS_INLINE
              const_iterator1 (const self_type &m, int rank, size_type i, size_type j, const vector_const_subiterator_type &itv, const const_subiterator_type &it):
                  container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
              BOOST_UBLAS_INLINE
              const_iterator1 (const iterator1 &it):
                  container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {}
  
              // Arithmetic
              BOOST_UBLAS_INLINE
              const_iterator1 &operator ++ () {
                  if (rank_ == 1 && layout_type::fast_i ())
                      ++ it_;
                  else {
                      const self_type &m = (*this) ();
                      if (rank_ == 0) {
                          ++ itv_;
                          i_ = itv_->first;
                      } else {
                          i_ = index1 () + 1;
                      }
                      if (rank_ == 1 && ++ itv_ == m.end1 ().itv_)
                          *this = m.find1 (rank_, i_, j_, 1);
                      else if (rank_ == 1) {
                          it_ = (*itv_).second.begin ();
                          if (it_ == (*itv_).second.end () || index2 () != j_)
                              *this = m.find1 (rank_, i_, j_, 1);
                      }
                  }
                  return *this;
              }
              BOOST_UBLAS_INLINE
              const_iterator1 &operator -- () {
                  if (rank_ == 1 && layout_type::fast_i ())
                      -- it_;
                  else {
                      const self_type &m = (*this) ();
                      if (rank_ == 0) {
                          -- itv_;
                          i_ = itv_->first;
                      } else {
                          i_ = index1 () - 1;
                      }
                      // FIXME: this expression should never become true!
                      if (rank_ == 1 && -- itv_ == m.end1 ().itv_)
                          *this = m.find1 (rank_, i_, j_, -1);
                      else if (rank_ == 1) {
                          it_ = (*itv_).second.begin ();
                          if (it_ == (*itv_).second.end () || index2 () != j_)
                              *this = m.find1 (rank_, i_, j_, -1);
                      }
                  }
                  return *this;
              }
  
              // Dereference
              BOOST_UBLAS_INLINE
              const_reference operator * () const {
                  BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                  BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
                  if (rank_ == 1) {
                      return (*it_).second;
                  } else {
                      return (*this) () (i_, j_);
                  }
              }
  
  #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_iterator2 begin () const {
                  const self_type &m = (*this) ();
                  return m.find2 (1, index1 (), 0);
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_iterator2 end () const {
                  const self_type &m = (*this) ();
                  return m.find2 (1, index1 (), m.size2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_reverse_iterator2 rbegin () const {
                  return const_reverse_iterator2 (end ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_reverse_iterator2 rend () const {
                  return const_reverse_iterator2 (begin ());
              }
  #endif
  
              // Indices
              BOOST_UBLAS_INLINE
              size_type index1 () const {
                  BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_M ((*itv_).first, (*it_).first) < (*this) ().size1 (), bad_index ());
                      return layout_type::index_M ((*itv_).first, (*it_).first);
                  } else {
                      return i_;
                  }
              }
              BOOST_UBLAS_INLINE
              size_type index2 () const {
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_m ((*itv_).first, (*it_).first) < (*this) ().size2 (), bad_index ());
                      return layout_type::index_m ((*itv_).first, (*it_).first);
                  } else {
                      return j_;
                  }
              }
  
              // Assignment
              BOOST_UBLAS_INLINE
              const_iterator1 &operator = (const const_iterator1 &it) {
                  container_const_reference<self_type>::assign (&it ());
                  rank_ = it.rank_;
                  i_ = it.i_;
                  j_ = it.j_;
                  itv_ = it.itv_;
                  it_ = it.it_;
                  return *this;
              }
  
              // Comparison
              BOOST_UBLAS_INLINE
              bool operator == (const const_iterator1 &it) const {
                  BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
                  // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
                  if (rank_ == 1 || it.rank_ == 1) {
                      return it_ == it.it_;
                  } else {
                      return i_ == it.i_ && j_ == it.j_;
                  }
              }
  
          private:
              int rank_;
              size_type i_;
              size_type j_;
              vector_const_subiterator_type itv_;
              const_subiterator_type it_;
          };
  
          BOOST_UBLAS_INLINE
          const_iterator1 begin1 () const {
              return find1 (0, 0, 0);
          }
          BOOST_UBLAS_INLINE
          const_iterator1 end1 () const {
              return find1 (0, size1_, 0);
          }
  
          class iterator1:
              public container_reference<mapped_vector_of_mapped_vector>,
              public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
                                                 iterator1, value_type> {
          public:
              typedef typename mapped_vector_of_mapped_vector::value_type value_type;
              typedef typename mapped_vector_of_mapped_vector::difference_type difference_type;
              typedef typename mapped_vector_of_mapped_vector::true_reference reference;
              typedef typename mapped_vector_of_mapped_vector::pointer pointer;
  
              typedef iterator2 dual_iterator_type;
              typedef reverse_iterator2 dual_reverse_iterator_type;
  
              // Construction and destruction
              BOOST_UBLAS_INLINE
              iterator1 ():
                  container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
              BOOST_UBLAS_INLINE
              iterator1 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it):
                  container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
  
              // Arithmetic
              BOOST_UBLAS_INLINE
              iterator1 &operator ++ () {
                  if (rank_ == 1 && layout_type::fast_i ())
                      ++ it_;
                  else {
                      self_type &m = (*this) ();
                      if (rank_ == 0) {
                          ++ itv_;
                          i_ = itv_->first;
                      } else {
                          i_ = index1 () + 1;
                      }
                      if (rank_ == 1 && ++ itv_ == m.end1 ().itv_)
                          *this = m.find1 (rank_, i_, j_, 1);
                      else if (rank_ == 1) {
                          it_ = (*itv_).second.begin ();
                          if (it_ == (*itv_).second.end () || index2 () != j_)
                              *this = m.find1 (rank_, i_, j_, 1);
                      }
                  }
                  return *this;
              }
              BOOST_UBLAS_INLINE
              iterator1 &operator -- () {
                  if (rank_ == 1 && layout_type::fast_i ())
                      -- it_;
                  else {
                      self_type &m = (*this) ();
                      if (rank_ == 0) {
                          -- itv_;
                          i_ = itv_->first;
                      } else {
                          i_ = index1 () - 1;
                      }
                      // FIXME: this expression should never become true!
                      if (rank_ == 1 && -- itv_ == m.end1 ().itv_)
                          *this = m.find1 (rank_, i_, j_, -1);
                      else if (rank_ == 1) {
                          it_ = (*itv_).second.begin ();
                          if (it_ == (*itv_).second.end () || index2 () != j_)
                              *this = m.find1 (rank_, i_, j_, -1);
                      }
                  }
                  return *this;
              }
  
              // Dereference
              BOOST_UBLAS_INLINE
              reference operator * () const {
                  BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                  BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
                  if (rank_ == 1) {
                      return (*it_).second;
                  } else {
                      return (*this) ().at_element (i_, j_);
                  }
              }
  
  #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              iterator2 begin () const {
                  self_type &m = (*this) ();
                  return m.find2 (1, index1 (), 0);
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              iterator2 end () const {
                  self_type &m = (*this) ();
                  return m.find2 (1, index1 (), m.size2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              reverse_iterator2 rbegin () const {
                  return reverse_iterator2 (end ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              reverse_iterator2 rend () const {
                  return reverse_iterator2 (begin ());
              }
  #endif
  
              // Indices
              BOOST_UBLAS_INLINE
              size_type index1 () const {
                  BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_M ((*itv_).first, (*it_).first) < (*this) ().size1 (), bad_index ());
                      return layout_type::index_M ((*itv_).first, (*it_).first);
                  } else {
                      return i_;
                  }
              }
              BOOST_UBLAS_INLINE
              size_type index2 () const {
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_m ((*itv_).first, (*it_).first) < (*this) ().size2 (), bad_index ());
                      return layout_type::index_m ((*itv_).first, (*it_).first);
                  } else {
                      return j_;
                  }
              }
  
              // Assignment
              BOOST_UBLAS_INLINE
              iterator1 &operator = (const iterator1 &it) {
                  container_reference<self_type>::assign (&it ());
                  rank_ = it.rank_;
                  i_ = it.i_;
                  j_ = it.j_;
                  itv_ = it.itv_;
                  it_ = it.it_;
                  return *this;
              }
  
              // Comparison
              BOOST_UBLAS_INLINE
              bool operator == (const iterator1 &it) const {
                  BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
                  // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
                  if (rank_ == 1 || it.rank_ == 1) {
                      return it_ == it.it_;
                  } else {
                      return i_ == it.i_ && j_ == it.j_;
                  }
              }
  
          private:
              int rank_;
              size_type i_;
              size_type j_;
              vector_subiterator_type itv_;
              subiterator_type it_;
  
              friend class const_iterator1;
          };
  
          BOOST_UBLAS_INLINE
          iterator1 begin1 () {
              return find1 (0, 0, 0);
          }
          BOOST_UBLAS_INLINE
          iterator1 end1 () {
              return find1 (0, size1_, 0);
          }
  
          class const_iterator2:
              public container_const_reference<mapped_vector_of_mapped_vector>,
              public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
                                                 const_iterator2, value_type> {
          public:
              typedef typename mapped_vector_of_mapped_vector::value_type value_type;
              typedef typename mapped_vector_of_mapped_vector::difference_type difference_type;
              typedef typename mapped_vector_of_mapped_vector::const_reference reference;
              typedef const typename mapped_vector_of_mapped_vector::pointer pointer;
  
              typedef const_iterator1 dual_iterator_type;
              typedef const_reverse_iterator1 dual_reverse_iterator_type;
  
              // Construction and destruction
              BOOST_UBLAS_INLINE
              const_iterator2 ():
                  container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
              BOOST_UBLAS_INLINE
              const_iterator2 (const self_type &m, int rank, size_type i, size_type j, const vector_const_subiterator_type &itv, const const_subiterator_type &it):
                  container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
              BOOST_UBLAS_INLINE
              const_iterator2 (const iterator2 &it):
                  container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {}
  
              // Arithmetic
              BOOST_UBLAS_INLINE
              const_iterator2 &operator ++ () {
                  if (rank_ == 1 && layout_type::fast_j ())
                      ++ it_;
                  else {
                      const self_type &m = (*this) ();
                      if (rank_ == 0) {
                          ++ itv_;
                          j_ = itv_->first;
                      } else {
                          j_ = index2 () + 1;
                      }
                      if (rank_ == 1 && ++ itv_ == m.end2 ().itv_)
                          *this = m.find2 (rank_, i_, j_, 1);
                      else if (rank_ == 1) {
                          it_ = (*itv_).second.begin ();
                          if (it_ == (*itv_).second.end () || index1 () != i_)
                              *this = m.find2 (rank_, i_, j_, 1);
                      }
                  }
                  return *this;
              }
              BOOST_UBLAS_INLINE
              const_iterator2 &operator -- () {
                  if (rank_ == 1 && layout_type::fast_j ())
                      -- it_;
                  else {
                      const self_type &m = (*this) ();
                      if (rank_ == 0) {
                          -- itv_;
                          j_ = itv_->first;
                      } else {
                          j_ = index2 () - 1;
                      }
                      // FIXME: this expression should never become true!
                      if (rank_ == 1 && -- itv_ == m.end2 ().itv_)
                          *this = m.find2 (rank_, i_, j_, -1);
                      else if (rank_ == 1) {
                          it_ = (*itv_).second.begin ();
                          if (it_ == (*itv_).second.end () || index1 () != i_)
                              *this = m.find2 (rank_, i_, j_, -1);
                      }
                  }
                  return *this;
              }
  
              // Dereference
              BOOST_UBLAS_INLINE
              const_reference operator * () const {
                  BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                  BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
                  if (rank_ == 1) {
                      return (*it_).second;
                  } else {
                      return (*this) () (i_, j_);
                  }
              }
  
  #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_iterator1 begin () const {
                  const self_type &m = (*this) ();
                  return m.find1 (1, 0, index2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_iterator1 end () const {
                  const self_type &m = (*this) ();
                  return m.find1 (1, m.size1 (), index2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_reverse_iterator1 rbegin () const {
                  return const_reverse_iterator1 (end ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_reverse_iterator1 rend () const {
                  return const_reverse_iterator1 (begin ());
              }
  #endif
  
              // Indices
              BOOST_UBLAS_INLINE
              size_type index1 () const {
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_M ((*itv_).first, (*it_).first) < (*this) ().size1 (), bad_index ());
                      return layout_type::index_M ((*itv_).first, (*it_).first);
                  } else {
                      return i_;
                  }
              }
              BOOST_UBLAS_INLINE
              size_type index2 () const {
                  BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_m ((*itv_).first, (*it_).first) < (*this) ().size2 (), bad_index ());
                      return layout_type::index_m ((*itv_).first, (*it_).first);
                  } else {
                      return j_;
                  }
              }
  
              // Assignment
              BOOST_UBLAS_INLINE
              const_iterator2 &operator = (const const_iterator2 &it) {
                  container_const_reference<self_type>::assign (&it ());
                  rank_ = it.rank_;
                  i_ = it.i_;
                  j_ = it.j_;
                  itv_ = it.itv_;
                  it_ = it.it_;
                  return *this;
              }
  
              // Comparison
              BOOST_UBLAS_INLINE
              bool operator == (const const_iterator2 &it) const {
                  BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
                  // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
                  if (rank_ == 1 || it.rank_ == 1) {
                      return it_ == it.it_;
                  } else {
                      return i_ == it.i_ && j_ == it.j_;
                  }
              }
  
          private:
              int rank_;
              size_type i_;
              size_type j_;
              vector_const_subiterator_type itv_;
              const_subiterator_type it_;
          };
  
          BOOST_UBLAS_INLINE
          const_iterator2 begin2 () const {
              return find2 (0, 0, 0);
          }
          BOOST_UBLAS_INLINE
          const_iterator2 end2 () const {
              return find2 (0, 0, size2_);
          }
  
          class iterator2:
              public container_reference<mapped_vector_of_mapped_vector>,
              public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
                                                 iterator2, value_type> {
          public:
              typedef typename mapped_vector_of_mapped_vector::value_type value_type;
              typedef typename mapped_vector_of_mapped_vector::difference_type difference_type;
              typedef typename mapped_vector_of_mapped_vector::true_reference reference;
              typedef typename mapped_vector_of_mapped_vector::pointer pointer;
  
              typedef iterator1 dual_iterator_type;
              typedef reverse_iterator1 dual_reverse_iterator_type;
  
              // Construction and destruction
              BOOST_UBLAS_INLINE
              iterator2 ():
                  container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
              BOOST_UBLAS_INLINE
              iterator2 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it):
                  container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
  
              // Arithmetic
              BOOST_UBLAS_INLINE
              iterator2 &operator ++ () {
                  if (rank_ == 1 && layout_type::fast_j ())
                      ++ it_;
                  else {
                      self_type &m = (*this) ();
                      if (rank_ == 0) {
                          ++ itv_;
                          j_ = itv_->first;
                      } else {
                          j_ = index2 () + 1;
                      }
                      if (rank_ == 1 && ++ itv_ == m.end2 ().itv_)
                          *this = m.find2 (rank_, i_, j_, 1);
                      else if (rank_ == 1) {
                          it_ = (*itv_).second.begin ();
                          if (it_ == (*itv_).second.end () || index1 () != i_)
                              *this = m.find2 (rank_, i_, j_, 1);
                      }
                  }
                  return *this;
              }
              BOOST_UBLAS_INLINE
              iterator2 &operator -- () {
                  if (rank_ == 1 && layout_type::fast_j ())
                      -- it_;
                  else {
                      self_type &m = (*this) ();
                      if (rank_ == 0) {
                          -- itv_;
                          j_ = itv_->first;
                      } else {
                          j_ = index2 () - 1;
                      }
                      // FIXME: this expression should never become true!
                      if (rank_ == 1 && -- itv_ == m.end2 ().itv_)
                          *this = m.find2 (rank_, i_, j_, -1);
                      else if (rank_ == 1) {
                          it_ = (*itv_).second.begin ();
                          if (it_ == (*itv_).second.end () || index1 () != i_)
                              *this = m.find2 (rank_, i_, j_, -1);
                      }
                  }
                  return *this;
              }
  
              // Dereference
              BOOST_UBLAS_INLINE
              reference operator * () const {
                  BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                  BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
                  if (rank_ == 1) {
                      return (*it_).second;
                  } else {
                      return (*this) ().at_element (i_, j_);
                  }
              }
  
  #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              iterator1 begin () const {
                  self_type &m = (*this) ();
                  return m.find1 (1, 0, index2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              iterator1 end () const {
                  self_type &m = (*this) ();
                  return m.find1 (1, m.size1 (), index2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              reverse_iterator1 rbegin () const {
                  return reverse_iterator1 (end ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              reverse_iterator1 rend () const {
                  return reverse_iterator1 (begin ());
              }
  #endif
  
              // Indices
              BOOST_UBLAS_INLINE
              size_type index1 () const {
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_M ((*itv_).first, (*it_).first) < (*this) ().size1 (), bad_index ());
                      return layout_type::index_M ((*itv_).first, (*it_).first);
                  } else {
                      return i_;
                  }
              }
              BOOST_UBLAS_INLINE
              size_type index2 () const {
                  BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_m ((*itv_).first, (*it_).first) < (*this) ().size2 (), bad_index ());
                      return layout_type::index_m ((*itv_).first, (*it_).first);
                  } else {
                      return j_;
                  }
              }
  
              // Assignment
              BOOST_UBLAS_INLINE
              iterator2 &operator = (const iterator2 &it) {
                  container_reference<self_type>::assign (&it ());
                  rank_ = it.rank_;
                  i_ = it.i_;
                  j_ = it.j_;
                  itv_ = it.itv_;
                  it_ = it.it_;
                  return *this;
              }
  
              // Comparison
              BOOST_UBLAS_INLINE
              bool operator == (const iterator2 &it) const {
                  BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
                  // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
                  if (rank_ == 1 || it.rank_ == 1) {
                      return it_ == it.it_;
                  } else {
                      return i_ == it.i_ && j_ == it.j_;
                  }
              }
  
          private:
              int rank_;
              size_type i_;
              size_type j_;
              vector_subiterator_type itv_;
              subiterator_type it_;
  
              friend class const_iterator2;
          };
  
          BOOST_UBLAS_INLINE
          iterator2 begin2 () {
              return find2 (0, 0, 0);
          }
          BOOST_UBLAS_INLINE
          iterator2 end2 () {
              return find2 (0, 0, size2_);
          }
  
          // Reverse iterators
  
          BOOST_UBLAS_INLINE
          const_reverse_iterator1 rbegin1 () const {
              return const_reverse_iterator1 (end1 ());
          }
          BOOST_UBLAS_INLINE
          const_reverse_iterator1 rend1 () const {
              return const_reverse_iterator1 (begin1 ());
          }
  
          BOOST_UBLAS_INLINE
          reverse_iterator1 rbegin1 () {
              return reverse_iterator1 (end1 ());
          }
          BOOST_UBLAS_INLINE
          reverse_iterator1 rend1 () {
              return reverse_iterator1 (begin1 ());
          }
  
          BOOST_UBLAS_INLINE
          const_reverse_iterator2 rbegin2 () const {
              return const_reverse_iterator2 (end2 ());
          }
          BOOST_UBLAS_INLINE
          const_reverse_iterator2 rend2 () const {
              return const_reverse_iterator2 (begin2 ());
          }
  
          BOOST_UBLAS_INLINE
          reverse_iterator2 rbegin2 () {
              return reverse_iterator2 (end2 ());
          }
          BOOST_UBLAS_INLINE
          reverse_iterator2 rend2 () {
              return reverse_iterator2 (begin2 ());
          }
  
           // Serialization
          template<class Archive>
          void serialize(Archive & ar, const unsigned int /* file_version */){
              serialization::collection_size_type s1 (size1_);
              serialization::collection_size_type s2 (size2_);
              ar & serialization::make_nvp("size1",s1);
              ar & serialization::make_nvp("size2",s2);
              if (Archive::is_loading::value) {
                  size1_ = s1;
                  size2_ = s2;
              }
              ar & serialization::make_nvp("data", data_);
          }
  
      private:
          size_type size1_;
          size_type size2_;
          array_type data_;
          static const value_type zero_;
      };
  
      template<class T, class L, class A>
      const typename mapped_vector_of_mapped_vector<T, L, A>::value_type mapped_vector_of_mapped_vector<T, L, A>::zero_ = value_type/*zero*/();
  
  
      // Comperssed array based sparse matrix class
      // Thanks to Kresimir Fresl for extending this to cover different index bases.
      template<class T, class L, std::size_t IB, class IA, class TA>
      class compressed_matrix:
          public matrix_container<compressed_matrix<T, L, IB, IA, TA> > {
  
          typedef T &true_reference;
          typedef T *pointer;
          typedef const T *const_pointer;
          typedef L layout_type;
          typedef compressed_matrix<T, L, IB, IA, TA> self_type;
      public:
  #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
          using matrix_container<self_type>::operator ();
  #endif
          // ISSUE require type consistency check
          // is_convertable (IA::size_type, TA::size_type)
          typedef typename IA::value_type size_type;
          // size_type for the data arrays.
          typedef typename IA::size_type array_size_type;
          // FIXME difference type for sparse storage iterators should it be in the container?
          typedef typename IA::difference_type difference_type;
          typedef T value_type;
          typedef const T &const_reference;
  #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
          typedef T &reference;
  #else
          typedef sparse_matrix_element<self_type> reference;
  #endif
          typedef IA index_array_type;
          typedef TA value_array_type;
          typedef const matrix_reference<const self_type> const_closure_type;
          typedef matrix_reference<self_type> closure_type;
          typedef compressed_vector<T, IB, IA, TA> vector_temporary_type;
          typedef self_type matrix_temporary_type;
          typedef sparse_tag storage_category;
          typedef typename L::orientation_category orientation_category;
  
          // Construction and destruction
          BOOST_UBLAS_INLINE
          compressed_matrix ():
              matrix_container<self_type> (),
              size1_ (0), size2_ (0), capacity_ (restrict_capacity (0)),
              filled1_ (1), filled2_ (0),
              index1_data_ (layout_type::size_M (size1_, size2_) + 1), index2_data_ (capacity_), value_data_ (capacity_) {
              index1_data_ [filled1_ - 1] = k_based (filled2_);
              storage_invariants ();
          }
          BOOST_UBLAS_INLINE
          compressed_matrix (size_type size1, size_type size2, size_type non_zeros = 0):
              matrix_container<self_type> (),
              size1_ (size1), size2_ (size2), capacity_ (restrict_capacity (non_zeros)),
              filled1_ (1), filled2_ (0),
              index1_data_ (layout_type::size_M (size1_, size2_) + 1), index2_data_ (capacity_), value_data_ (capacity_) {
              index1_data_ [filled1_ - 1] = k_based (filled2_);
              storage_invariants ();
          }
          BOOST_UBLAS_INLINE
          compressed_matrix (const compressed_matrix &m):
              matrix_container<self_type> (),
              size1_ (m.size1_), size2_ (m.size2_), capacity_ (m.capacity_),
              filled1_ (m.filled1_), filled2_ (m.filled2_),
              index1_data_ (m.index1_data_), index2_data_ (m.index2_data_), value_data_ (m.value_data_) {
              storage_invariants ();
          }
           
          BOOST_UBLAS_INLINE
          compressed_matrix (const coordinate_matrix<T, L, IB, IA, TA> &m):
              matrix_container<self_type> (),
              size1_ (m.size1()), size2_ (m.size2()),
              index1_data_ (layout_type::size_M (size1_, size2_) + 1)
          {
              m.sort();
              reserve(m.nnz(), false);
              filled2_ = m.nnz();
              const_subiterator_type  i_start = m.index1_data().begin();
              const_subiterator_type  i_end   = (i_start + filled2_);
              const_subiterator_type  i = i_start;
              size_type r = 1;
              for (; (r < layout_type::size_M (size1_, size2_)) && (i != i_end); ++r) {
                  i = std::lower_bound(i, i_end, r);
                  index1_data_[r] = k_based( i - i_start );
              }
              filled1_ = r + 1;
              std::copy( m.index2_data().begin(), m.index2_data().begin() + filled2_, index2_data_.begin());
              std::copy( m.value_data().begin(), m.value_data().begin() + filled2_, value_data_.begin());
              index1_data_ [filled1_ - 1] = k_based(filled2_);
              storage_invariants ();
          }
  
         template<class AE>
         BOOST_UBLAS_INLINE
         compressed_matrix (const matrix_expression<AE> &ae, size_type non_zeros = 0):
              matrix_container<self_type> (),
              size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), capacity_ (restrict_capacity (non_zeros)),
              filled1_ (1), filled2_ (0),
              index1_data_ (layout_type::size_M (ae ().size1 (), ae ().size2 ()) + 1),
              index2_data_ (capacity_), value_data_ (capacity_) {
              index1_data_ [filled1_ - 1] = k_based (filled2_);
              storage_invariants ();
              matrix_assign<scalar_assign> (*this, ae);
          }
  
          // Accessors
          BOOST_UBLAS_INLINE
          size_type size1 () const {
              return size1_;
          }
          BOOST_UBLAS_INLINE
          size_type size2 () const {
              return size2_;
          }
          BOOST_UBLAS_INLINE
          size_type nnz_capacity () const {
              return capacity_;
          }
          BOOST_UBLAS_INLINE
          size_type nnz () const {
              return filled2_;
          }
          
          // Storage accessors
          BOOST_UBLAS_INLINE
          static size_type index_base () {
              return IB;
          }
          BOOST_UBLAS_INLINE
          array_size_type filled1 () const {
              return filled1_;
          }
          BOOST_UBLAS_INLINE
          array_size_type filled2 () const {
              return filled2_;
          }
          BOOST_UBLAS_INLINE
          const index_array_type &index1_data () const {
              return index1_data_;
          }
          BOOST_UBLAS_INLINE
          const index_array_type &index2_data () const {
              return index2_data_;
          }
          BOOST_UBLAS_INLINE
          const value_array_type &value_data () const {
              return value_data_;
          }
          BOOST_UBLAS_INLINE
          void set_filled (const array_size_type& filled1, const array_size_type& filled2) {
              filled1_ = filled1;
              filled2_ = filled2;
              storage_invariants ();
          }
          BOOST_UBLAS_INLINE
          index_array_type &index1_data () {
              return index1_data_;
          }
          BOOST_UBLAS_INLINE
          index_array_type &index2_data () {
              return index2_data_;
          }
          BOOST_UBLAS_INLINE
          value_array_type &value_data () {
              return value_data_;
          }
          BOOST_UBLAS_INLINE
          void complete_index1_data () {
              while (filled1_ <= layout_type::size_M (size1_, size2_)) {
                  this->index1_data_ [filled1_] = k_based (filled2_);
                  ++ this->filled1_;
              }
          }
  
          // Resizing
      private:
          BOOST_UBLAS_INLINE
          size_type restrict_capacity (size_type non_zeros) const {
              non_zeros = (std::max) (non_zeros, (std::min) (size1_, size2_));
              // Guarding against overflow - Thanks to Alexei Novakov for the hint.
              // non_zeros = (std::min) (non_zeros, size1_ * size2_);
              if (size1_ > 0 && non_zeros / size1_ >= size2_)
                  non_zeros = size1_ * size2_;
              return non_zeros;
          }
      public:
          BOOST_UBLAS_INLINE
          void resize (size_type size1, size_type size2, bool preserve = true) {
              // FIXME preserve unimplemented
              BOOST_UBLAS_CHECK (!preserve, internal_logic ());
              size1_ = size1;
              size2_ = size2;
              capacity_ = restrict_capacity (capacity_);
              filled1_ = 1;
              filled2_ = 0;
              index1_data_.resize (layout_type::size_M (size1_, size2_) + 1);
              index2_data_.resize (capacity_);
              value_data_.resize (capacity_);
              index1_data_ [filled1_ - 1] = k_based (filled2_);
              storage_invariants ();
          }
  
          // Reserving
          BOOST_UBLAS_INLINE
          void reserve (size_type non_zeros, bool preserve = true) {
              capacity_ = restrict_capacity (non_zeros);
              if (preserve) {
                  index2_data_.resize (capacity_, size_type ());
                  value_data_.resize (capacity_, value_type ());
                  filled2_ = (std::min) (capacity_, filled2_);
              }
              else {
                  index2_data_.resize (capacity_);
                  value_data_.resize (capacity_);
                  filled1_ = 1;
                  filled2_ = 0;
                  index1_data_ [filled1_ - 1] = k_based (filled2_);
              }
              storage_invariants ();
         }
  
          // Element support
          BOOST_UBLAS_INLINE
          pointer find_element (size_type i, size_type j) {
              return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i, j));
          }
          BOOST_UBLAS_INLINE
          const_pointer find_element (size_type i, size_type j) const {
              size_type element1 (layout_type::index_M (i, j));
              size_type element2 (layout_type::index_m (i, j));
              if (filled1_ <= element1 + 1)
                  return 0;
              vector_const_subiterator_type itv (index1_data_.begin () + element1);
              const_subiterator_type it_begin (index2_data_.begin () + zero_based (*itv));
              const_subiterator_type it_end (index2_data_.begin () + zero_based (*(itv + 1)));
              const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (element2), std::less<size_type> ()));
              if (it == it_end || *it != k_based (element2))
                  return 0;
              return &value_data_ [it - index2_data_.begin ()];
          }
  
          // Element access
          BOOST_UBLAS_INLINE
          const_reference operator () (size_type i, size_type j) const {
              const_pointer p = find_element (i, j);
              if (p)
                  return *p;
              else
                  return zero_;
          }
          BOOST_UBLAS_INLINE
          reference operator () (size_type i, size_type j) {
  #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
              size_type element1 (layout_type::index_M (i, j));
              size_type element2 (layout_type::index_m (i, j));
              if (filled1_ <= element1 + 1)
                  return insert_element (i, j, value_type/*zero*/());
              pointer p = find_element (i, j);
              if (p)
                  return *p;
              else
                  return insert_element (i, j, value_type/*zero*/());
  #else
              return reference (*this, i, j);
  #endif
          }
  
          // Element assignment
          BOOST_UBLAS_INLINE
          true_reference insert_element (size_type i, size_type j, const_reference t) {
              BOOST_UBLAS_CHECK (!find_element (i, j), bad_index ());        // duplicate element
              if (filled2_ >= capacity_)
                  reserve (2 * filled2_, true);
              BOOST_UBLAS_CHECK (filled2_ < capacity_, internal_logic ());
              size_type element1 = layout_type::index_M (i, j);
              size_type element2 = layout_type::index_m (i, j);
              while (filled1_ <= element1 + 1) {
                  index1_data_ [filled1_] = k_based (filled2_);
                  ++ filled1_;
              }
              vector_subiterator_type itv (index1_data_.begin () + element1);
              subiterator_type it_begin (index2_data_.begin () + zero_based (*itv));
              subiterator_type it_end (index2_data_.begin () + zero_based (*(itv + 1)));
              subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (element2), std::less<size_type> ()));
              typename std::iterator_traits<subiterator_type>::difference_type n = it - index2_data_.begin ();
              BOOST_UBLAS_CHECK (it == it_end || *it != k_based (element2), internal_logic ());   // duplicate bound by lower_bound
              ++ filled2_;
              it = index2_data_.begin () + n;
              std::copy_backward (it, index2_data_.begin () + filled2_ - 1, index2_data_.begin () + filled2_);
              *it = k_based (element2);
              typename value_array_type::iterator itt (value_data_.begin () + n);
              std::copy_backward (itt, value_data_.begin () + filled2_ - 1, value_data_.begin () + filled2_);
              *itt = t;
              while (element1 + 1 < filled1_) {
                  ++ index1_data_ [element1 + 1];
                  ++ element1;
              }
              storage_invariants ();
              return *itt;
          }
          BOOST_UBLAS_INLINE
          void erase_element (size_type i, size_type j) {
              size_type element1 = layout_type::index_M (i, j);
              size_type element2 = layout_type::index_m (i, j);
              if (element1 + 1 >= filled1_)
                  return;
              vector_subiterator_type itv (index1_data_.begin () + element1);
              subiterator_type it_begin (index2_data_.begin () + zero_based (*itv));
              subiterator_type it_end (index2_data_.begin () + zero_based (*(itv + 1)));
              subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (element2), std::less<size_type> ()));
              if (it != it_end && *it == k_based (element2)) {
                  typename std::iterator_traits<subiterator_type>::difference_type n = it - index2_data_.begin ();
                  std::copy (it + 1, index2_data_.begin () + filled2_, it);
                  typename value_array_type::iterator itt (value_data_.begin () + n);
                  std::copy (itt + 1, value_data_.begin () + filled2_, itt);
                  -- filled2_;
                  while (index1_data_ [filled1_ - 2] > k_based (filled2_)) {
                      index1_data_ [filled1_ - 1] = 0;
                      -- filled1_;
                  }
                  while (element1 + 1 < filled1_) {
                      -- index1_data_ [element1 + 1];
                      ++ element1;
                  }
              }
              storage_invariants ();
          }
          
          // Zeroing
          BOOST_UBLAS_INLINE
          void clear () {
              filled1_ = 1;
              filled2_ = 0;
              index1_data_ [filled1_ - 1] = k_based (filled2_);
              storage_invariants ();
          }
  
          // Assignment
          BOOST_UBLAS_INLINE
          compressed_matrix &operator = (const compressed_matrix &m) {
              if (this != &m) {
                  size1_ = m.size1_;
                  size2_ = m.size2_;
                  capacity_ = m.capacity_;
                  filled1_ = m.filled1_;
                  filled2_ = m.filled2_;
                  index1_data_ = m.index1_data_;
                  index2_data_ = m.index2_data_;
                  value_data_ = m.value_data_;
              }
              storage_invariants ();
              return *this;
          }
          template<class C>          // Container assignment without temporary
          BOOST_UBLAS_INLINE
          compressed_matrix &operator = (const matrix_container<C> &m) {
              resize (m ().size1 (), m ().size2 (), false);
              assign (m);
              return *this;
          }
          BOOST_UBLAS_INLINE
          compressed_matrix &assign_temporary (compressed_matrix &m) {
              swap (m);
              return *this;
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          compressed_matrix &operator = (const matrix_expression<AE> &ae) {
              self_type temporary (ae, capacity_);
              return assign_temporary (temporary);
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          compressed_matrix &assign (const matrix_expression<AE> &ae) {
              matrix_assign<scalar_assign> (*this, ae);
              return *this;
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          compressed_matrix& operator += (const matrix_expression<AE> &ae) {
              self_type temporary (*this + ae, capacity_);
              return assign_temporary (temporary);
          }
          template<class C>          // Container assignment without temporary
          BOOST_UBLAS_INLINE
          compressed_matrix &operator += (const matrix_container<C> &m) {
              plus_assign (m);
              return *this;
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          compressed_matrix &plus_assign (const matrix_expression<AE> &ae) {
              matrix_assign<scalar_plus_assign> (*this, ae);
              return *this;
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          compressed_matrix& operator -= (const matrix_expression<AE> &ae) {
              self_type temporary (*this - ae, capacity_);
              return assign_temporary (temporary);
          }
          template<class C>          // Container assignment without temporary
          BOOST_UBLAS_INLINE
          compressed_matrix &operator -= (const matrix_container<C> &m) {
              minus_assign (m);
              return *this;
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          compressed_matrix &minus_assign (const matrix_expression<AE> &ae) {
              matrix_assign<scalar_minus_assign> (*this, ae);
              return *this;
          }
          template<class AT>
          BOOST_UBLAS_INLINE
          compressed_matrix& operator *= (const AT &at) {
              matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
              return *this;
          }
          template<class AT>
          BOOST_UBLAS_INLINE
          compressed_matrix& operator /= (const AT &at) {
              matrix_assign_scalar<scalar_divides_assign> (*this, at);
              return *this;
          }
  
          // Swapping
          BOOST_UBLAS_INLINE
          void swap (compressed_matrix &m) {
              if (this != &m) {
                  std::swap (size1_, m.size1_);
                  std::swap (size2_, m.size2_);
                  std::swap (capacity_, m.capacity_);
                  std::swap (filled1_, m.filled1_);
                  std::swap (filled2_, m.filled2_);
                  index1_data_.swap (m.index1_data_);
                  index2_data_.swap (m.index2_data_);
                  value_data_.swap (m.value_data_);
              }
              storage_invariants ();
          }
          BOOST_UBLAS_INLINE
          friend void swap (compressed_matrix &m1, compressed_matrix &m2) {
              m1.swap (m2);
          }
  
          // Back element insertion and erasure
          BOOST_UBLAS_INLINE
          void push_back (size_type i, size_type j, const_reference t) {
              if (filled2_ >= capacity_)
                  reserve (2 * filled2_, true);
              BOOST_UBLAS_CHECK (filled2_ < capacity_, internal_logic ());
              size_type element1 = layout_type::index_M (i, j);
              size_type element2 = layout_type::index_m (i, j);
              while (filled1_ < element1 + 2) {
                  index1_data_ [filled1_] = k_based (filled2_);
                  ++ filled1_;
              }
              // must maintain sort order
              BOOST_UBLAS_CHECK ((filled1_ == element1 + 2 &&
                                  (filled2_ == zero_based (index1_data_ [filled1_ - 2]) ||
                                  index2_data_ [filled2_ - 1] < k_based (element2))), external_logic ());
              ++ filled2_;
              index1_data_ [filled1_ - 1] = k_based (filled2_);
              index2_data_ [filled2_ - 1] = k_based (element2);
              value_data_ [filled2_ - 1] = t;
              storage_invariants ();
          }
          BOOST_UBLAS_INLINE
          void pop_back () {
              BOOST_UBLAS_CHECK (filled1_ > 0 && filled2_ > 0, external_logic ());
              -- filled2_;
              while (index1_data_ [filled1_ - 2] > k_based (filled2_)) {
                  index1_data_ [filled1_ - 1] = 0;
                  -- filled1_;
              }
              -- index1_data_ [filled1_ - 1];
              storage_invariants ();
          }
  
          // Iterator types
      private:
          // Use index array iterator
          typedef typename IA::const_iterator vector_const_subiterator_type;
          typedef typename IA::iterator vector_subiterator_type;
          typedef typename IA::const_iterator const_subiterator_type;
          typedef typename IA::iterator subiterator_type;
  
          BOOST_UBLAS_INLINE
          true_reference at_element (size_type i, size_type j) {
              pointer p = find_element (i, j);
              BOOST_UBLAS_CHECK (p, bad_index ());
              return *p;
          }
  
      public:
          class const_iterator1;
          class iterator1;
          class const_iterator2;
          class iterator2;
          typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
          typedef reverse_iterator_base1<iterator1> reverse_iterator1;
          typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
          typedef reverse_iterator_base2<iterator2> reverse_iterator2;
  
          // Element lookup
          // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
          const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const {
              for (;;) {
                  array_size_type address1 (layout_type::index_M (i, j));
                  array_size_type address2 (layout_type::index_m (i, j));
                  vector_const_subiterator_type itv (index1_data_.begin () + (std::min) (filled1_ - 1, address1));
                  if (filled1_ <= address1 + 1)
                      return const_iterator1 (*this, rank, i, j, itv, index2_data_.begin () + filled2_);
  
                  const_subiterator_type it_begin (index2_data_.begin () + zero_based (*itv));
                  const_subiterator_type it_end (index2_data_.begin () + zero_based (*(itv + 1)));
  
                  const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ()));
                  if (rank == 0)
                      return const_iterator1 (*this, rank, i, j, itv, it);
                  if (it != it_end && zero_based (*it) == address2)
                      return const_iterator1 (*this, rank, i, j, itv, it);
                  if (direction > 0) {
                      if (layout_type::fast_i ()) {
                          if (it == it_end)
                              return const_iterator1 (*this, rank, i, j, itv, it);
                          i = zero_based (*it);
                      } else {
                          if (i >= size1_)
                              return const_iterator1 (*this, rank, i, j, itv, it);
                          ++ i;
                      }
                  } else /* if (direction < 0)  */ {
                      if (layout_type::fast_i ()) {
                          if (it == index2_data_.begin () + zero_based (*itv))
                              return const_iterator1 (*this, rank, i, j, itv, it);
                          i = zero_based (*(it - 1));
                      } else {
                          if (i == 0)
                              return const_iterator1 (*this, rank, i, j, itv, it);
                          -- i;
                      }
                  }
              }
          }
          // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
          iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) {
              for (;;) {
                  array_size_type address1 (layout_type::index_M (i, j));
                  array_size_type address2 (layout_type::index_m (i, j));
                  vector_subiterator_type itv (index1_data_.begin () + (std::min) (filled1_ - 1, address1));
                  if (filled1_ <= address1 + 1)
                      return iterator1 (*this, rank, i, j, itv, index2_data_.begin () + filled2_);
  
                  subiterator_type it_begin (index2_data_.begin () + zero_based (*itv));
                  subiterator_type it_end (index2_data_.begin () + zero_based (*(itv + 1)));
  
                  subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ()));
                  if (rank == 0)
                      return iterator1 (*this, rank, i, j, itv, it);
                  if (it != it_end && zero_based (*it) == address2)
                      return iterator1 (*this, rank, i, j, itv, it);
                  if (direction > 0) {
                      if (layout_type::fast_i ()) {
                          if (it == it_end)
                              return iterator1 (*this, rank, i, j, itv, it);
                          i = zero_based (*it);
                      } else {
                          if (i >= size1_)
                              return iterator1 (*this, rank, i, j, itv, it);
                          ++ i;
                      }
                  } else /* if (direction < 0)  */ {
                      if (layout_type::fast_i ()) {
                          if (it == index2_data_.begin () + zero_based (*itv))
                              return iterator1 (*this, rank, i, j, itv, it);
                          i = zero_based (*(it - 1));
                      } else {
                          if (i == 0)
                              return iterator1 (*this, rank, i, j, itv, it);
                          -- i;
                      }
                  }
              }
          }
          // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
          const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const {
              for (;;) {
                  array_size_type address1 (layout_type::index_M (i, j));
                  array_size_type address2 (layout_type::index_m (i, j));
                  vector_const_subiterator_type itv (index1_data_.begin () + (std::min) (filled1_ - 1, address1));
                  if (filled1_ <= address1 + 1)
                      return const_iterator2 (*this, rank, i, j, itv, index2_data_.begin () + filled2_);
  
                  const_subiterator_type it_begin (index2_data_.begin () + zero_based (*itv));
                  const_subiterator_type it_end (index2_data_.begin () + zero_based (*(itv + 1)));
  
                  const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ()));
                  if (rank == 0)
                      return const_iterator2 (*this, rank, i, j, itv, it);
                  if (it != it_end && zero_based (*it) == address2)
                      return const_iterator2 (*this, rank, i, j, itv, it);
                  if (direction > 0) {
                      if (layout_type::fast_j ()) {
                          if (it == it_end)
                              return const_iterator2 (*this, rank, i, j, itv, it);
                          j = zero_based (*it);
                      } else {
                          if (j >= size2_)
                              return const_iterator2 (*this, rank, i, j, itv, it);
                          ++ j;
                      }
                  } else /* if (direction < 0)  */ {
                      if (layout_type::fast_j ()) {
                          if (it == index2_data_.begin () + zero_based (*itv))
                              return const_iterator2 (*this, rank, i, j, itv, it);
                          j = zero_based (*(it - 1));
                      } else {
                          if (j == 0)
                              return const_iterator2 (*this, rank, i, j, itv, it);
                          -- j;
                      }
                  }
              }
          }
          // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
          iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) {
              for (;;) {
                  array_size_type address1 (layout_type::index_M (i, j));
                  array_size_type address2 (layout_type::index_m (i, j));
                  vector_subiterator_type itv (index1_data_.begin () + (std::min) (filled1_ - 1, address1));
                  if (filled1_ <= address1 + 1)
                      return iterator2 (*this, rank, i, j, itv, index2_data_.begin () + filled2_);
  
                  subiterator_type it_begin (index2_data_.begin () + zero_based (*itv));
                  subiterator_type it_end (index2_data_.begin () + zero_based (*(itv + 1)));
  
                  subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ()));
                  if (rank == 0)
                      return iterator2 (*this, rank, i, j, itv, it);
                  if (it != it_end && zero_based (*it) == address2)
                      return iterator2 (*this, rank, i, j, itv, it);
                  if (direction > 0) {
                      if (layout_type::fast_j ()) {
                          if (it == it_end)
                              return iterator2 (*this, rank, i, j, itv, it);
                          j = zero_based (*it);
                      } else {
                          if (j >= size2_)
                              return iterator2 (*this, rank, i, j, itv, it);
                          ++ j;
                      }
                  } else /* if (direction < 0)  */ {
                      if (layout_type::fast_j ()) {
                          if (it == index2_data_.begin () + zero_based (*itv))
                              return iterator2 (*this, rank, i, j, itv, it);
                          j = zero_based (*(it - 1));
                      } else {
                          if (j == 0)
                              return iterator2 (*this, rank, i, j, itv, it);
                          -- j;
                      }
                  }
              }
          }
  
  
          class const_iterator1:
              public container_const_reference<compressed_matrix>,
              public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
                                                 const_iterator1, value_type> {
          public:
              typedef typename compressed_matrix::value_type value_type;
              typedef typename compressed_matrix::difference_type difference_type;
              typedef typename compressed_matrix::const_reference reference;
              typedef const typename compressed_matrix::pointer pointer;
  
              typedef const_iterator2 dual_iterator_type;
              typedef const_reverse_iterator2 dual_reverse_iterator_type;
  
              // Construction and destruction
              BOOST_UBLAS_INLINE
              const_iterator1 ():
                  container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
              BOOST_UBLAS_INLINE
              const_iterator1 (const self_type &m, int rank, size_type i, size_type j, const vector_const_subiterator_type &itv, const const_subiterator_type &it):
                  container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
              BOOST_UBLAS_INLINE
              const_iterator1 (const iterator1 &it):
                  container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {}
  
              // Arithmetic
              BOOST_UBLAS_INLINE
              const_iterator1 &operator ++ () {
                  if (rank_ == 1 && layout_type::fast_i ())
                      ++ it_;
                  else {
                      i_ = index1 () + 1;
                      if (rank_ == 1)
                          *this = (*this) ().find1 (rank_, i_, j_, 1);
                  }
                  return *this;
              }
              BOOST_UBLAS_INLINE
              const_iterator1 &operator -- () {
                  if (rank_ == 1 && layout_type::fast_i ())
                      -- it_;
                  else {
                      --i_;
                      if (rank_ == 1)
                          *this = (*this) ().find1 (rank_, i_, j_, -1);
                  }
                  return *this;
              }
  
              // Dereference
              BOOST_UBLAS_INLINE
              const_reference operator * () const {
                  BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                  BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
                  if (rank_ == 1) {
                      return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
                  } else {
                      return (*this) () (i_, j_);
                  }
              }
  
  #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_iterator2 begin () const {
                  const self_type &m = (*this) ();
                  return m.find2 (1, index1 (), 0);
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_iterator2 end () const {
                  const self_type &m = (*this) ();
                  return m.find2 (1, index1 (), m.size2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_reverse_iterator2 rbegin () const {
                  return const_reverse_iterator2 (end ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_reverse_iterator2 rend () const {
                  return const_reverse_iterator2 (begin ());
              }
  #endif
  
              // Indices
              BOOST_UBLAS_INLINE
              size_type index1 () const {
                  BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
                      return layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
                  } else {
                      return i_;
                  }
              }
              BOOST_UBLAS_INLINE
              size_type index2 () const {
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
                      return layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
                  } else {
                      return j_;
                  }
              }
  
              // Assignment
              BOOST_UBLAS_INLINE
              const_iterator1 &operator = (const const_iterator1 &it) {
                  container_const_reference<self_type>::assign (&it ());
                  rank_ = it.rank_;
                  i_ = it.i_;
                  j_ = it.j_;
                  itv_ = it.itv_;
                  it_ = it.it_;
                  return *this;
              }
  
              // Comparison
              BOOST_UBLAS_INLINE
              bool operator == (const const_iterator1 &it) const {
                  BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
                  // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
                  if (rank_ == 1 || it.rank_ == 1) {
                      return it_ == it.it_;
                  } else {
                      return i_ == it.i_ && j_ == it.j_;
                  }
              }
  
          private:
              int rank_;
              size_type i_;
              size_type j_;
              vector_const_subiterator_type itv_;
              const_subiterator_type it_;
          };
  
          BOOST_UBLAS_INLINE
          const_iterator1 begin1 () const {
              return find1 (0, 0, 0);
          }
          BOOST_UBLAS_INLINE
          const_iterator1 end1 () const {
              return find1 (0, size1_, 0);
          }
  
          class iterator1:
              public container_reference<compressed_matrix>,
              public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
                                                 iterator1, value_type> {
          public:
              typedef typename compressed_matrix::value_type value_type;
              typedef typename compressed_matrix::difference_type difference_type;
              typedef typename compressed_matrix::true_reference reference;
              typedef typename compressed_matrix::pointer pointer;
  
              typedef iterator2 dual_iterator_type;
              typedef reverse_iterator2 dual_reverse_iterator_type;
  
              // Construction and destruction
              BOOST_UBLAS_INLINE
              iterator1 ():
                  container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
              BOOST_UBLAS_INLINE
              iterator1 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it):
                  container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
  
              // Arithmetic
              BOOST_UBLAS_INLINE
              iterator1 &operator ++ () {
                  if (rank_ == 1 && layout_type::fast_i ())
                      ++ it_;
                  else {
                      i_ = index1 () + 1;
                      if (rank_ == 1)
                          *this = (*this) ().find1 (rank_, i_, j_, 1);
                  }
                  return *this;
              }
              BOOST_UBLAS_INLINE
              iterator1 &operator -- () {
                  if (rank_ == 1 && layout_type::fast_i ())
                      -- it_;
                  else {
                      --i_;
                      if (rank_ == 1)
                          *this = (*this) ().find1 (rank_, i_, j_, -1);
                  }
                  return *this;
              }
  
              // Dereference
              BOOST_UBLAS_INLINE
              reference operator * () const {
                  BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                  BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
                  if (rank_ == 1) {
                      return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
                  } else {
                      return (*this) ().at_element (i_, j_);
                  }
              }
  
  #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              iterator2 begin () const {
                  self_type &m = (*this) ();
                  return m.find2 (1, index1 (), 0);
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              iterator2 end () const {
                  self_type &m = (*this) ();
                  return m.find2 (1, index1 (), m.size2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              reverse_iterator2 rbegin () const {
                  return reverse_iterator2 (end ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              reverse_iterator2 rend () const {
                  return reverse_iterator2 (begin ());
              }
  #endif
  
              // Indices
              BOOST_UBLAS_INLINE
              size_type index1 () const {
                  BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
                      return layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
                  } else {
                      return i_;
                  }
              }
              BOOST_UBLAS_INLINE
              size_type index2 () const {
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
                      return layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
                  } else {
                      return j_;
                  }
              }
  
              // Assignment
              BOOST_UBLAS_INLINE
              iterator1 &operator = (const iterator1 &it) {
                  container_reference<self_type>::assign (&it ());
                  rank_ = it.rank_;
                  i_ = it.i_;
                  j_ = it.j_;
                  itv_ = it.itv_;
                  it_ = it.it_;
                  return *this;
              }
  
              // Comparison
              BOOST_UBLAS_INLINE
              bool operator == (const iterator1 &it) const {
                  BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
                  // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
                  if (rank_ == 1 || it.rank_ == 1) {
                      return it_ == it.it_;
                  } else {
                      return i_ == it.i_ && j_ == it.j_;
                  }
              }
  
          private:
              int rank_;
              size_type i_;
              size_type j_;
              vector_subiterator_type itv_;
              subiterator_type it_;
  
              friend class const_iterator1;
          };
  
          BOOST_UBLAS_INLINE
          iterator1 begin1 () {
              return find1 (0, 0, 0);
          }
          BOOST_UBLAS_INLINE
          iterator1 end1 () {
              return find1 (0, size1_, 0);
          }
  
          class const_iterator2:
              public container_const_reference<compressed_matrix>,
              public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
                                                 const_iterator2, value_type> {
          public:
              typedef typename compressed_matrix::value_type value_type;
              typedef typename compressed_matrix::difference_type difference_type;
              typedef typename compressed_matrix::const_reference reference;
              typedef const typename compressed_matrix::pointer pointer;
  
              typedef const_iterator1 dual_iterator_type;
              typedef const_reverse_iterator1 dual_reverse_iterator_type;
  
              // Construction and destruction
              BOOST_UBLAS_INLINE
              const_iterator2 ():
                  container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
              BOOST_UBLAS_INLINE
              const_iterator2 (const self_type &m, int rank, size_type i, size_type j, const vector_const_subiterator_type itv, const const_subiterator_type &it):
                  container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
              BOOST_UBLAS_INLINE
              const_iterator2 (const iterator2 &it):
                  container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {}
  
              // Arithmetic
              BOOST_UBLAS_INLINE
              const_iterator2 &operator ++ () {
                  if (rank_ == 1 && layout_type::fast_j ())
                      ++ it_;
                  else {
                      j_ = index2 () + 1;
                      if (rank_ == 1)
                          *this = (*this) ().find2 (rank_, i_, j_, 1);
                  }
                  return *this;
              }
              BOOST_UBLAS_INLINE
              const_iterator2 &operator -- () {
                  if (rank_ == 1 && layout_type::fast_j ())
                      -- it_;
                  else {
                      --j_;
                      if (rank_ == 1)
                          *this = (*this) ().find2 (rank_, i_, j_, -1);
                  }
                  return *this;
              }
  
              // Dereference
              BOOST_UBLAS_INLINE
              const_reference operator * () const {
                  BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                  BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
                  if (rank_ == 1) {
                      return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
                  } else {
                      return (*this) () (i_, j_);
                  }
              }
  
  #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_iterator1 begin () const {
                  const self_type &m = (*this) ();
                  return m.find1 (1, 0, index2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_iterator1 end () const {
                  const self_type &m = (*this) ();
                  return m.find1 (1, m.size1 (), index2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_reverse_iterator1 rbegin () const {
                  return const_reverse_iterator1 (end ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_reverse_iterator1 rend () const {
                  return const_reverse_iterator1 (begin ());
              }
  #endif
  
              // Indices
              BOOST_UBLAS_INLINE
              size_type index1 () const {
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
                      return layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
                  } else {
                      return i_;
                  }
              }
              BOOST_UBLAS_INLINE
              size_type index2 () const {
                  BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
                      return layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
                  } else {
                      return j_;
                  }
              }
  
              // Assignment
              BOOST_UBLAS_INLINE
              const_iterator2 &operator = (const const_iterator2 &it) {
                  container_const_reference<self_type>::assign (&it ());
                  rank_ = it.rank_;
                  i_ = it.i_;
                  j_ = it.j_;
                  itv_ = it.itv_;
                  it_ = it.it_;
                  return *this;
              }
  
              // Comparison
              BOOST_UBLAS_INLINE
              bool operator == (const const_iterator2 &it) const {
                  BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
                  // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
                  if (rank_ == 1 || it.rank_ == 1) {
                      return it_ == it.it_;
                  } else {
                      return i_ == it.i_ && j_ == it.j_;
                  }
              }
  
          private:
              int rank_;
              size_type i_;
              size_type j_;
              vector_const_subiterator_type itv_;
              const_subiterator_type it_;
          };
  
          BOOST_UBLAS_INLINE
          const_iterator2 begin2 () const {
              return find2 (0, 0, 0);
          }
          BOOST_UBLAS_INLINE
          const_iterator2 end2 () const {
              return find2 (0, 0, size2_);
          }
  
          class iterator2:
              public container_reference<compressed_matrix>,
              public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
                                                 iterator2, value_type> {
          public:
              typedef typename compressed_matrix::value_type value_type;
              typedef typename compressed_matrix::difference_type difference_type;
              typedef typename compressed_matrix::true_reference reference;
              typedef typename compressed_matrix::pointer pointer;
  
              typedef iterator1 dual_iterator_type;
              typedef reverse_iterator1 dual_reverse_iterator_type;
  
              // Construction and destruction
              BOOST_UBLAS_INLINE
              iterator2 ():
                  container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
              BOOST_UBLAS_INLINE
              iterator2 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it):
                  container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
  
              // Arithmetic
              BOOST_UBLAS_INLINE
              iterator2 &operator ++ () {
                  if (rank_ == 1 && layout_type::fast_j ())
                      ++ it_;
                  else {
                      j_ = index2 () + 1;
                      if (rank_ == 1)
                          *this = (*this) ().find2 (rank_, i_, j_, 1);
                  }
                  return *this;
              }
              BOOST_UBLAS_INLINE
              iterator2 &operator -- () {
                  if (rank_ == 1 && layout_type::fast_j ())
                      -- it_;
                  else {
                      --j_;
                      if (rank_ == 1)
                          *this = (*this) ().find2 (rank_, i_, j_, -1);
                  }
                  return *this;
              }
  
              // Dereference
              BOOST_UBLAS_INLINE
              reference operator * () const {
                  BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                  BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
                  if (rank_ == 1) {
                      return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
                  } else {
                      return (*this) ().at_element (i_, j_);
                  }
              }
  
  #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              iterator1 begin () const {
                  self_type &m = (*this) ();
                  return m.find1 (1, 0, index2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              iterator1 end () const {
                  self_type &m = (*this) ();
                  return m.find1 (1, m.size1 (), index2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              reverse_iterator1 rbegin () const {
                  return reverse_iterator1 (end ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              reverse_iterator1 rend () const {
                  return reverse_iterator1 (begin ());
              }
  #endif
  
              // Indices
              BOOST_UBLAS_INLINE
              size_type index1 () const {
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
                      return layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
                  } else {
                      return i_;
                  }
              }
              BOOST_UBLAS_INLINE
              size_type index2 () const {
                  BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
                      return layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
                  } else {
                      return j_;
                  }
              }
  
              // Assignment
              BOOST_UBLAS_INLINE
              iterator2 &operator = (const iterator2 &it) {
                  container_reference<self_type>::assign (&it ());
                  rank_ = it.rank_;
                  i_ = it.i_;
                  j_ = it.j_;
                  itv_ = it.itv_;
                  it_ = it.it_;
                  return *this;
              }
  
              // Comparison
              BOOST_UBLAS_INLINE
              bool operator == (const iterator2 &it) const {
                  BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
                  // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
                  if (rank_ == 1 || it.rank_ == 1) {
                      return it_ == it.it_;
                  } else {
                      return i_ == it.i_ && j_ == it.j_;
                  }
              }
  
          private:
              int rank_;
              size_type i_;
              size_type j_;
              vector_subiterator_type itv_;
              subiterator_type it_;
  
              friend class const_iterator2;
          };
  
          BOOST_UBLAS_INLINE
          iterator2 begin2 () {
              return find2 (0, 0, 0);
          }
          BOOST_UBLAS_INLINE
          iterator2 end2 () {
              return find2 (0, 0, size2_);
          }
  
          // Reverse iterators
  
          BOOST_UBLAS_INLINE
          const_reverse_iterator1 rbegin1 () const {
              return const_reverse_iterator1 (end1 ());
          }
          BOOST_UBLAS_INLINE
          const_reverse_iterator1 rend1 () const {
              return const_reverse_iterator1 (begin1 ());
          }
  
          BOOST_UBLAS_INLINE
          reverse_iterator1 rbegin1 () {
              return reverse_iterator1 (end1 ());
          }
          BOOST_UBLAS_INLINE
          reverse_iterator1 rend1 () {
              return reverse_iterator1 (begin1 ());
          }
  
          BOOST_UBLAS_INLINE
          const_reverse_iterator2 rbegin2 () const {
              return const_reverse_iterator2 (end2 ());
          }
          BOOST_UBLAS_INLINE
          const_reverse_iterator2 rend2 () const {
              return const_reverse_iterator2 (begin2 ());
          }
  
          BOOST_UBLAS_INLINE
          reverse_iterator2 rbegin2 () {
              return reverse_iterator2 (end2 ());
          }
          BOOST_UBLAS_INLINE
          reverse_iterator2 rend2 () {
              return reverse_iterator2 (begin2 ());
          }
  
           // Serialization
          template<class Archive>
          void serialize(Archive & ar, const unsigned int /* file_version */){
              serialization::collection_size_type s1 (size1_);
              serialization::collection_size_type s2 (size2_);
              ar & serialization::make_nvp("size1",s1);
              ar & serialization::make_nvp("size2",s2);
              if (Archive::is_loading::value) {
                  size1_ = s1;
                  size2_ = s2;
              }
              ar & serialization::make_nvp("capacity", capacity_);
              ar & serialization::make_nvp("filled1", filled1_);
              ar & serialization::make_nvp("filled2", filled2_);
              ar & serialization::make_nvp("index1_data", index1_data_);
              ar & serialization::make_nvp("index2_data", index2_data_);
              ar & serialization::make_nvp("value_data", value_data_);
              storage_invariants();
          }
  
      private:
          void storage_invariants () const {
              BOOST_UBLAS_CHECK (layout_type::size_M (size1_, size2_) + 1 == index1_data_.size (), internal_logic ());
              BOOST_UBLAS_CHECK (capacity_ == index2_data_.size (), internal_logic ());
              BOOST_UBLAS_CHECK (capacity_ == value_data_.size (), internal_logic ());
              BOOST_UBLAS_CHECK (filled1_ > 0 && filled1_ <= layout_type::size_M (size1_, size2_) + 1, internal_logic ());
              BOOST_UBLAS_CHECK (filled2_ <= capacity_, internal_logic ());
              BOOST_UBLAS_CHECK (index1_data_ [filled1_ - 1] == k_based (filled2_), internal_logic ());
          }
          
          size_type size1_;
          size_type size2_;
          array_size_type capacity_;
          array_size_type filled1_;
          array_size_type filled2_;
          index_array_type index1_data_;
          index_array_type index2_data_;
          value_array_type value_data_;
          static const value_type zero_;
  
          BOOST_UBLAS_INLINE
          static size_type zero_based (size_type k_based_index) {
              return k_based_index - IB;
          }
          BOOST_UBLAS_INLINE
          static size_type k_based (size_type zero_based_index) {
              return zero_based_index + IB;
          }
  
          friend class iterator1;
          friend class iterator2;
          friend class const_iterator1;
          friend class const_iterator2;
      };
  
      template<class T, class L, std::size_t IB, class IA, class TA>
      const typename compressed_matrix<T, L, IB, IA, TA>::value_type compressed_matrix<T, L, IB, IA, TA>::zero_ = value_type/*zero*/();
  
  
      // Coordinate array based sparse matrix class
      // Thanks to Kresimir Fresl for extending this to cover different index bases.
      template<class T, class L, std::size_t IB, class IA, class TA>
      class coordinate_matrix:
          public matrix_container<coordinate_matrix<T, L, IB, IA, TA> > {
  
          typedef T &true_reference;
          typedef T *pointer;
          typedef const T *const_pointer;
          typedef L layout_type;
          typedef coordinate_matrix<T, L, IB, IA, TA> self_type;
      public:
  #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
          using matrix_container<self_type>::operator ();
  #endif
          // ISSUE require type consistency check, is_convertable (IA::size_type, TA::size_type)
          typedef typename IA::value_type size_type;
          // ISSUE difference_type cannot be deduced for sparse indices, we only know the value_type
          typedef std::ptrdiff_t difference_type;
          // size_type for the data arrays.
          typedef typename IA::size_type array_size_type;
          typedef T value_type;
          typedef const T &const_reference;
  #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
          typedef T &reference;
  #else
          typedef sparse_matrix_element<self_type> reference;
  #endif
          typedef IA index_array_type;
          typedef TA value_array_type;
          typedef const matrix_reference<const self_type> const_closure_type;
          typedef matrix_reference<self_type> closure_type;
          typedef coordinate_vector<T, IB, IA, TA> vector_temporary_type;
          typedef self_type matrix_temporary_type;
          typedef sparse_tag storage_category;
          typedef typename L::orientation_category orientation_category;
  
          // Construction and destruction
          BOOST_UBLAS_INLINE
          coordinate_matrix ():
              matrix_container<self_type> (),
              size1_ (0), size2_ (0), capacity_ (restrict_capacity (0)),
              filled_ (0), sorted_filled_ (filled_), sorted_ (true),
              index1_data_ (capacity_), index2_data_ (capacity_), value_data_ (capacity_) {
              storage_invariants ();
          }
          BOOST_UBLAS_INLINE
          coordinate_matrix (size_type size1, size_type size2, array_size_type non_zeros = 0):
              matrix_container<self_type> (),
              size1_ (size1), size2_ (size2), capacity_ (restrict_capacity (non_zeros)),
              filled_ (0), sorted_filled_ (filled_), sorted_ (true),
              index1_data_ (capacity_), index2_data_ (capacity_), value_data_ (capacity_) {
              storage_invariants ();
          }
          BOOST_UBLAS_INLINE
          coordinate_matrix (const coordinate_matrix &m):
              matrix_container<self_type> (),
              size1_ (m.size1_), size2_ (m.size2_), capacity_ (m.capacity_),
              filled_ (m.filled_), sorted_filled_ (m.sorted_filled_), sorted_ (m.sorted_),
              index1_data_ (m.index1_data_), index2_data_ (m.index2_data_), value_data_ (m.value_data_) {
              storage_invariants ();
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          coordinate_matrix (const matrix_expression<AE> &ae, array_size_type non_zeros = 0):
              matrix_container<self_type> (),
              size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), capacity_ (restrict_capacity (non_zeros)),
              filled_ (0), sorted_filled_ (filled_), sorted_ (true),
              index1_data_ (capacity_), index2_data_ (capacity_), value_data_ (capacity_) {
              storage_invariants ();
              matrix_assign<scalar_assign> (*this, ae);
          }
  
          // Accessors
          BOOST_UBLAS_INLINE
          size_type size1 () const {
              return size1_;
          }
          BOOST_UBLAS_INLINE
          size_type size2 () const {
              return size2_;
          }
          BOOST_UBLAS_INLINE
          size_type nnz_capacity () const {
              return capacity_;
          }
          BOOST_UBLAS_INLINE
          size_type nnz () const {
              return filled_;
          }
  
          // Storage accessors
          BOOST_UBLAS_INLINE
          static size_type index_base () {
              return IB;
          }
          BOOST_UBLAS_INLINE
          array_size_type filled () const {
              return filled_;
          }
          BOOST_UBLAS_INLINE
          const index_array_type &index1_data () const {
              return index1_data_;
          }
          BOOST_UBLAS_INLINE
          const index_array_type &index2_data () const {
              return index2_data_;
          }
          BOOST_UBLAS_INLINE
          const value_array_type &value_data () const {
              return value_data_;
          }
          BOOST_UBLAS_INLINE
          void set_filled (const array_size_type &filled) {
              // Make sure that storage_invariants() succeeds
              if (sorted_ && filled < filled_)
                  sorted_filled_ = filled;
              else
                  sorted_ = (sorted_filled_ == filled);
              filled_ = filled;
              storage_invariants ();
          }
          BOOST_UBLAS_INLINE
          index_array_type &index1_data () {
              return index1_data_;
          }
          BOOST_UBLAS_INLINE
          index_array_type &index2_data () {
              return index2_data_;
          }
          BOOST_UBLAS_INLINE
          value_array_type &value_data () {
              return value_data_;
          }
  
          // Resizing
      private:
          BOOST_UBLAS_INLINE
          array_size_type restrict_capacity (array_size_type non_zeros) const {
              // minimum non_zeros
              non_zeros = (std::max) (non_zeros, array_size_type((std::min) (size1_, size2_)));
              // ISSUE no maximum as coordinate may contain inserted duplicates
              return non_zeros;
          }
      public:
          BOOST_UBLAS_INLINE
          void resize (size_type size1, size_type size2, bool preserve = true) {
              // FIXME preserve unimplemented
              BOOST_UBLAS_CHECK (!preserve, internal_logic ());
              size1_ = size1;
              size2_ = size2;
              capacity_ = restrict_capacity (capacity_);
              index1_data_.resize (capacity_);
              index2_data_.resize (capacity_);
              value_data_.resize (capacity_);
              filled_ = 0;
              sorted_filled_ = filled_;
              sorted_ = true;
              storage_invariants ();
          }
  
          // Reserving
          BOOST_UBLAS_INLINE
          void reserve (array_size_type non_zeros, bool preserve = true) {
              sort ();    // remove duplicate elements
              capacity_ = restrict_capacity (non_zeros);
              if (preserve) {
                  index1_data_.resize (capacity_, size_type ());
                  index2_data_.resize (capacity_, size_type ());
                  value_data_.resize (capacity_, value_type ());
                  filled_ = (std::min) (capacity_, filled_);
              }
              else {
                  index1_data_.resize (capacity_);
                  index2_data_.resize (capacity_);
                  value_data_.resize (capacity_);
                  filled_ = 0;
              }
              sorted_filled_ = filled_;
              storage_invariants ();
          }
  
          // Element support
          BOOST_UBLAS_INLINE
          pointer find_element (size_type i, size_type j) {
              return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i, j));
          }
          BOOST_UBLAS_INLINE
          const_pointer find_element (size_type i, size_type j) const {
              sort ();
              size_type element1 (layout_type::index_M (i, j));
              size_type element2 (layout_type::index_m (i, j));
              vector_const_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (element1), std::less<size_type> ()));
              vector_const_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (element1), std::less<size_type> ()));
              if (itv_begin == itv_end)
                  return 0;
              const_subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ()));
              const_subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ()));
              const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (element2), std::less<size_type> ()));
              if (it == it_end || *it != k_based (element2))
                  return 0;
              return &value_data_ [it - index2_data_.begin ()];
          }
  
          // Element access
          BOOST_UBLAS_INLINE
          const_reference operator () (size_type i, size_type j) const {
              const_pointer p = find_element (i, j);
              if (p)
                  return *p;
              else 
                  return zero_;
          }
          BOOST_UBLAS_INLINE
          reference operator () (size_type i, size_type j) {
  #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
              pointer p = find_element (i, j);
              if (p)
                  return *p;
              else
                  return insert_element (i, j, value_type/*zero*/());
  #else
              return reference (*this, i, j);
  #endif
          }
  
          // Element assignment
          BOOST_UBLAS_INLINE
          void append_element (size_type i, size_type j, const_reference t) {
              if (filled_ >= capacity_)
                  reserve (2 * filled_, true);
              BOOST_UBLAS_CHECK (filled_ < capacity_, internal_logic ());
              size_type element1 = layout_type::index_M (i, j);
              size_type element2 = layout_type::index_m (i, j);
              index1_data_ [filled_] = k_based (element1);
              index2_data_ [filled_] = k_based (element2);
              value_data_ [filled_] = t;
              ++ filled_;
              sorted_ = false;
              storage_invariants ();
          }
          BOOST_UBLAS_INLINE
          true_reference insert_element (size_type i, size_type j, const_reference t) {
              BOOST_UBLAS_CHECK (!find_element (i, j), bad_index ());        // duplicate element
              append_element (i, j, t);
              return value_data_ [filled_ - 1];
          }
          BOOST_UBLAS_INLINE
          void erase_element (size_type i, size_type j) {
              size_type element1 = layout_type::index_M (i, j);
              size_type element2 = layout_type::index_m (i, j);
              sort ();
              vector_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (element1), std::less<size_type> ()));
              vector_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (element1), std::less<size_type> ()));
              subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ()));
              subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ()));
              subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (element2), std::less<size_type> ()));
              if (it != it_end && *it == k_based (element2)) {
                  typename std::iterator_traits<subiterator_type>::difference_type n = it - index2_data_.begin ();
                  vector_subiterator_type itv (index1_data_.begin () + n);
                  std::copy (itv + 1, index1_data_.begin () + filled_, itv);
                  std::copy (it + 1, index2_data_.begin () + filled_, it);
                  typename value_array_type::iterator itt (value_data_.begin () + n);
                  std::copy (itt + 1, value_data_.begin () + filled_, itt);
                  -- filled_;
                  sorted_filled_ = filled_;
              }
              storage_invariants ();
          }
          
          // Zeroing
          BOOST_UBLAS_INLINE
          void clear () {
              filled_ = 0;
              sorted_filled_ = filled_;
              sorted_ = true;
              storage_invariants ();
          }
  
          // Assignment
          BOOST_UBLAS_INLINE
          coordinate_matrix &operator = (const coordinate_matrix &m) {
              if (this != &m) {
                  size1_ = m.size1_;
                  size2_ = m.size2_;
                  capacity_ = m.capacity_;
                  filled_ = m.filled_;
                  sorted_filled_ = m.sorted_filled_;
                  sorted_ = m.sorted_;
                  index1_data_ = m.index1_data_;
                  index2_data_ = m.index2_data_;
                  value_data_ = m.value_data_;
                  BOOST_UBLAS_CHECK (capacity_ == index1_data_.size (), internal_logic ());
                  BOOST_UBLAS_CHECK (capacity_ == index2_data_.size (), internal_logic ());
                  BOOST_UBLAS_CHECK (capacity_ == value_data_.size (), internal_logic ());
              }
              storage_invariants ();
              return *this;
          }
          template<class C>          // Container assignment without temporary
          BOOST_UBLAS_INLINE
          coordinate_matrix &operator = (const matrix_container<C> &m) {
              resize (m ().size1 (), m ().size2 (), false);
              assign (m);
              return *this;
          }
          BOOST_UBLAS_INLINE
          coordinate_matrix &assign_temporary (coordinate_matrix &m) {
              swap (m);
              return *this;
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          coordinate_matrix &operator = (const matrix_expression<AE> &ae) {
              self_type temporary (ae, capacity_);
              return assign_temporary (temporary);
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          coordinate_matrix &assign (const matrix_expression<AE> &ae) {
              matrix_assign<scalar_assign> (*this, ae);
              return *this;
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          coordinate_matrix& operator += (const matrix_expression<AE> &ae) {
              self_type temporary (*this + ae, capacity_);
              return assign_temporary (temporary);
          }
          template<class C>          // Container assignment without temporary
          BOOST_UBLAS_INLINE
          coordinate_matrix &operator += (const matrix_container<C> &m) {
              plus_assign (m);
              return *this;
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          coordinate_matrix &plus_assign (const matrix_expression<AE> &ae) {
              matrix_assign<scalar_plus_assign> (*this, ae);
              return *this;
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          coordinate_matrix& operator -= (const matrix_expression<AE> &ae) {
              self_type temporary (*this - ae, capacity_);
              return assign_temporary (temporary);
          }
          template<class C>          // Container assignment without temporary
          BOOST_UBLAS_INLINE
          coordinate_matrix &operator -= (const matrix_container<C> &m) {
              minus_assign (m);
              return *this;
          }
          template<class AE>
          BOOST_UBLAS_INLINE
          coordinate_matrix &minus_assign (const matrix_expression<AE> &ae) {
              matrix_assign<scalar_minus_assign> (*this, ae);
              return *this;
          }
          template<class AT>
          BOOST_UBLAS_INLINE
          coordinate_matrix& operator *= (const AT &at) {
              matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
              return *this;
          }
          template<class AT>
          BOOST_UBLAS_INLINE
          coordinate_matrix& operator /= (const AT &at) {
              matrix_assign_scalar<scalar_divides_assign> (*this, at);
              return *this;
          }
  
          // Swapping
          BOOST_UBLAS_INLINE
          void swap (coordinate_matrix &m) {
              if (this != &m) {
                  std::swap (size1_, m.size1_);
                  std::swap (size2_, m.size2_);
                  std::swap (capacity_, m.capacity_);
                  std::swap (filled_, m.filled_);
                  std::swap (sorted_filled_, m.sorted_filled_);
                  std::swap (sorted_, m.sorted_);
                  index1_data_.swap (m.index1_data_);
                  index2_data_.swap (m.index2_data_);
                  value_data_.swap (m.value_data_);
              }
              storage_invariants ();
          }
          BOOST_UBLAS_INLINE
          friend void swap (coordinate_matrix &m1, coordinate_matrix &m2) {
              m1.swap (m2);
          }
  
          // Sorting and summation of duplicates
          BOOST_UBLAS_INLINE
          void sort () const {
              if (! sorted_ && filled_ > 0) {
                  typedef index_triple_array<index_array_type, index_array_type, value_array_type> array_triple;
                  array_triple ita (filled_, index1_data_, index2_data_, value_data_);
                  const typename array_triple::iterator iunsorted = ita.begin () + sorted_filled_;
                  // sort new elements and merge
                  std::sort (iunsorted, ita.end ());
                  std::inplace_merge (ita.begin (), iunsorted, ita.end ());
                  
                  // sum duplicates with += and remove
                  array_size_type filled = 0;
                  for (array_size_type i = 1; i < filled_; ++ i) {
                      if (index1_data_ [filled] != index1_data_ [i] ||
                          index2_data_ [filled] != index2_data_ [i]) {
                          ++ filled;
                          if (filled != i) {
                              index1_data_ [filled] = index1_data_ [i];
                              index2_data_ [filled] = index2_data_ [i];
                              value_data_ [filled] = value_data_ [i];
                          }
                      } else {
                          value_data_ [filled] += value_data_ [i];
                      }
                  }
                  filled_ = filled + 1;
                  sorted_filled_ = filled_;
                  sorted_ = true;
                  storage_invariants ();
              }
          }
  
          // Back element insertion and erasure
          BOOST_UBLAS_INLINE
          void push_back (size_type i, size_type j, const_reference t) {
              size_type element1 = layout_type::index_M (i, j);
              size_type element2 = layout_type::index_m (i, j);
              // must maintain sort order
              BOOST_UBLAS_CHECK (sorted_ && 
                      (filled_ == 0 ||
                      index1_data_ [filled_ - 1] < k_based (element1) ||
                      (index1_data_ [filled_ - 1] == k_based (element1) && index2_data_ [filled_ - 1] < k_based (element2)))
                      , external_logic ());
              if (filled_ >= capacity_)
                  reserve (2 * filled_, true);
              BOOST_UBLAS_CHECK (filled_ < capacity_, internal_logic ());
              index1_data_ [filled_] = k_based (element1);
              index2_data_ [filled_] = k_based (element2);
              value_data_ [filled_] = t;
              ++ filled_;
              sorted_filled_ = filled_;
              storage_invariants ();
          }
          BOOST_UBLAS_INLINE
          void pop_back () {
              // ISSUE invariants could be simpilfied if sorted required as precondition
              BOOST_UBLAS_CHECK (filled_ > 0, external_logic ());
              -- filled_;
              sorted_filled_ = (std::min) (sorted_filled_, filled_);
              sorted_ = sorted_filled_ = filled_;
              storage_invariants ();
          }
  
          // Iterator types
      private:
          // Use index array iterator
          typedef typename IA::const_iterator vector_const_subiterator_type;
          typedef typename IA::iterator vector_subiterator_type;
          typedef typename IA::const_iterator const_subiterator_type;
          typedef typename IA::iterator subiterator_type;
  
          BOOST_UBLAS_INLINE
          true_reference at_element (size_type i, size_type j) {
              pointer p = find_element (i, j);
              BOOST_UBLAS_CHECK (p, bad_index ());
              return *p;
          }
  
      public:
          class const_iterator1;
          class iterator1;
          class const_iterator2;
          class iterator2;
          typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
          typedef reverse_iterator_base1<iterator1> reverse_iterator1;
          typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
          typedef reverse_iterator_base2<iterator2> reverse_iterator2;
  
          // Element lookup
          // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
          const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const {
              sort ();
              for (;;) {
                  size_type address1 (layout_type::index_M (i, j));
                  size_type address2 (layout_type::index_m (i, j));
                  vector_const_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ()));
                  vector_const_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ()));
  
                  const_subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ()));
                  const_subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ()));
  
                  const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ()));
                  vector_const_subiterator_type itv (index1_data_.begin () + (it - index2_data_.begin ()));
                  if (rank == 0)
                      return const_iterator1 (*this, rank, i, j, itv, it);
                  if (it != it_end && zero_based (*it) == address2)
                      return const_iterator1 (*this, rank, i, j, itv, it);
                  if (direction > 0) {
                      if (layout_type::fast_i ()) {
                          if (it == it_end)
                              return const_iterator1 (*this, rank, i, j, itv, it);
                          i = zero_based (*it);
                      } else {
                          if (i >= size1_)
                              return const_iterator1 (*this, rank, i, j, itv, it);
                          ++ i;
                      }
                  } else /* if (direction < 0)  */ {
                      if (layout_type::fast_i ()) {
                          if (it == index2_data_.begin () + array_size_type (zero_based (*itv)))
                              return const_iterator1 (*this, rank, i, j, itv, it);
                          i = zero_based (*(it - 1));
                      } else {
                          if (i == 0)
                              return const_iterator1 (*this, rank, i, j, itv, it);
                          -- i;
                      }
                  }
              }
          }
          // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
          iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) {
              sort ();
              for (;;) {
                  size_type address1 (layout_type::index_M (i, j));
                  size_type address2 (layout_type::index_m (i, j));
                  vector_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ()));
                  vector_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ()));
  
                  subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ()));
                  subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ()));
  
                  subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ()));
                  vector_subiterator_type itv (index1_data_.begin () + (it - index2_data_.begin ()));
                  if (rank == 0)
                      return iterator1 (*this, rank, i, j, itv, it);
                  if (it != it_end && zero_based (*it) == address2)
                      return iterator1 (*this, rank, i, j, itv, it);
                  if (direction > 0) {
                      if (layout_type::fast_i ()) {
                          if (it == it_end)
                              return iterator1 (*this, rank, i, j, itv, it);
                          i = zero_based (*it);
                      } else {
                          if (i >= size1_)
                              return iterator1 (*this, rank, i, j, itv, it);
                          ++ i;
                      }
                  } else /* if (direction < 0)  */ {
                      if (layout_type::fast_i ()) {
                          if (it == index2_data_.begin () + array_size_type (zero_based (*itv)))
                              return iterator1 (*this, rank, i, j, itv, it);
                          i = zero_based (*(it - 1));
                      } else {
                          if (i == 0)
                              return iterator1 (*this, rank, i, j, itv, it);
                          -- i;
                      }
                  }
              }
          }
          // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
          const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const {
              sort ();
              for (;;) {
                  size_type address1 (layout_type::index_M (i, j));
                  size_type address2 (layout_type::index_m (i, j));
                  vector_const_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ()));
                  vector_const_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ()));
  
                  const_subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ()));
                  const_subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ()));
  
                  const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ()));
                  vector_const_subiterator_type itv (index1_data_.begin () + (it - index2_data_.begin ()));
                  if (rank == 0)
                      return const_iterator2 (*this, rank, i, j, itv, it);
                  if (it != it_end && zero_based (*it) == address2)
                      return const_iterator2 (*this, rank, i, j, itv, it);
                  if (direction > 0) {
                      if (layout_type::fast_j ()) {
                          if (it == it_end)
                              return const_iterator2 (*this, rank, i, j, itv, it);
                          j = zero_based (*it);
                      } else {
                          if (j >= size2_)
                              return const_iterator2 (*this, rank, i, j, itv, it);
                          ++ j;
                      }
                  } else /* if (direction < 0)  */ {
                      if (layout_type::fast_j ()) {
                          if (it == index2_data_.begin () + array_size_type (zero_based (*itv)))
                              return const_iterator2 (*this, rank, i, j, itv, it);
                          j = zero_based (*(it - 1));
                      } else {
                          if (j == 0)
                              return const_iterator2 (*this, rank, i, j, itv, it);
                          -- j;
                      }
                  }
              }
          }
          // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
          iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) {
              sort ();
              for (;;) {
                  size_type address1 (layout_type::index_M (i, j));
                  size_type address2 (layout_type::index_m (i, j));
                  vector_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ()));
                  vector_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ()));
  
                  subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ()));
                  subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ()));
  
                  subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ()));
                  vector_subiterator_type itv (index1_data_.begin () + (it - index2_data_.begin ()));
                  if (rank == 0)
                      return iterator2 (*this, rank, i, j, itv, it);
                  if (it != it_end && zero_based (*it) == address2)
                      return iterator2 (*this, rank, i, j, itv, it);
                  if (direction > 0) {
                      if (layout_type::fast_j ()) {
                          if (it == it_end)
                              return iterator2 (*this, rank, i, j, itv, it);
                          j = zero_based (*it);
                      } else {
                          if (j >= size2_)
                              return iterator2 (*this, rank, i, j, itv, it);
                          ++ j;
                      }
                  } else /* if (direction < 0)  */ {
                      if (layout_type::fast_j ()) {
                          if (it == index2_data_.begin () + array_size_type (zero_based (*itv)))
                              return iterator2 (*this, rank, i, j, itv, it);
                          j = zero_based (*(it - 1));
                      } else {
                          if (j == 0)
                              return iterator2 (*this, rank, i, j, itv, it);
                          -- j;
                      }
                  }
              }
          }
  
  
          class const_iterator1:
              public container_const_reference<coordinate_matrix>,
              public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
                                                 const_iterator1, value_type> {
          public:
              typedef typename coordinate_matrix::value_type value_type;
              typedef typename coordinate_matrix::difference_type difference_type;
              typedef typename coordinate_matrix::const_reference reference;
              typedef const typename coordinate_matrix::pointer pointer;
  
              typedef const_iterator2 dual_iterator_type;
              typedef const_reverse_iterator2 dual_reverse_iterator_type;
  
              // Construction and destruction
              BOOST_UBLAS_INLINE
              const_iterator1 ():
                  container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
              BOOST_UBLAS_INLINE
              const_iterator1 (const self_type &m, int rank, size_type i, size_type j, const vector_const_subiterator_type &itv, const const_subiterator_type &it):
                  container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
              BOOST_UBLAS_INLINE
              const_iterator1 (const iterator1 &it):
                  container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {}
  
              // Arithmetic
              BOOST_UBLAS_INLINE
              const_iterator1 &operator ++ () {
                  if (rank_ == 1 && layout_type::fast_i ())
                      ++ it_;
                  else {
                      i_ = index1 () + 1;
                      if (rank_ == 1)
                          *this = (*this) ().find1 (rank_, i_, j_, 1);
                  }
                  return *this;
              }
              BOOST_UBLAS_INLINE
              const_iterator1 &operator -- () {
                  if (rank_ == 1 && layout_type::fast_i ())
                      -- it_;
                  else {
                      i_ = index1 () - 1;
                      if (rank_ == 1)
                          *this = (*this) ().find1 (rank_, i_, j_, -1);
                  }
                  return *this;
              }
  
              // Dereference
              BOOST_UBLAS_INLINE
              const_reference operator * () const {
                  BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                  BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
                  if (rank_ == 1) {
                      return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
                  } else {
                      return (*this) () (i_, j_);
                  }
              }
  
  #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_iterator2 begin () const {
                  const self_type &m = (*this) ();
                  return m.find2 (1, index1 (), 0);
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_iterator2 end () const {
                  const self_type &m = (*this) ();
                  return m.find2 (1, index1 (), m.size2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_reverse_iterator2 rbegin () const {
                  return const_reverse_iterator2 (end ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_reverse_iterator2 rend () const {
                  return const_reverse_iterator2 (begin ());
              }
  #endif
  
              // Indices
              BOOST_UBLAS_INLINE
              size_type index1 () const {
                  BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
                      return layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
                  } else {
                      return i_;
                  }
              }
              BOOST_UBLAS_INLINE
              size_type index2 () const {
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
                      return layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
                  } else {
                      return j_;
                  }
              }
  
              // Assignment
              BOOST_UBLAS_INLINE
              const_iterator1 &operator = (const const_iterator1 &it) {
                  container_const_reference<self_type>::assign (&it ());
                  rank_ = it.rank_;
                  i_ = it.i_;
                  j_ = it.j_;
                  itv_ = it.itv_;
                  it_ = it.it_;
                  return *this;
              }
  
              // Comparison
              BOOST_UBLAS_INLINE
              bool operator == (const const_iterator1 &it) const {
                  BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
                  // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
                  if (rank_ == 1 || it.rank_ == 1) {
                      return it_ == it.it_;
                  } else {
                      return i_ == it.i_ && j_ == it.j_;
                  }
              }
  
          private:
              int rank_;
              size_type i_;
              size_type j_;
              vector_const_subiterator_type itv_;
              const_subiterator_type it_;
          };
  
          BOOST_UBLAS_INLINE
          const_iterator1 begin1 () const {
              return find1 (0, 0, 0);
          }
          BOOST_UBLAS_INLINE
          const_iterator1 end1 () const {
              return find1 (0, size1_, 0);
          }
  
          class iterator1:
              public container_reference<coordinate_matrix>,
              public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
                                                 iterator1, value_type> {
          public:
              typedef typename coordinate_matrix::value_type value_type;
              typedef typename coordinate_matrix::difference_type difference_type;
              typedef typename coordinate_matrix::true_reference reference;
              typedef typename coordinate_matrix::pointer pointer;
  
              typedef iterator2 dual_iterator_type;
              typedef reverse_iterator2 dual_reverse_iterator_type;
  
              // Construction and destruction
              BOOST_UBLAS_INLINE
              iterator1 ():
                  container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
              BOOST_UBLAS_INLINE
              iterator1 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it):
                  container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
  
              // Arithmetic
              BOOST_UBLAS_INLINE
              iterator1 &operator ++ () {
                  if (rank_ == 1 && layout_type::fast_i ())
                      ++ it_;
                  else {
                      i_ = index1 () + 1;
                      if (rank_ == 1)
                          *this = (*this) ().find1 (rank_, i_, j_, 1);
                  }
                  return *this;
              }
              BOOST_UBLAS_INLINE
              iterator1 &operator -- () {
                  if (rank_ == 1 && layout_type::fast_i ())
                      -- it_;
                  else {
                      i_ = index1 () - 1;
                      if (rank_ == 1)
                          *this = (*this) ().find1 (rank_, i_, j_, -1);
                  }
                  return *this;
              }
  
              // Dereference
              BOOST_UBLAS_INLINE
              reference operator * () const {
                  BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                  BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
                  if (rank_ == 1) {
                      return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
                  } else {
                      return (*this) ().at_element (i_, j_);
                  }
              }
  
  #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              iterator2 begin () const {
                  self_type &m = (*this) ();
                  return m.find2 (1, index1 (), 0);
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              iterator2 end () const {
                  self_type &m = (*this) ();
                  return m.find2 (1, index1 (), m.size2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              reverse_iterator2 rbegin () const {
                  return reverse_iterator2 (end ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              reverse_iterator2 rend () const {
                  return reverse_iterator2 (begin ());
              }
  #endif
  
              // Indices
              BOOST_UBLAS_INLINE
              size_type index1 () const {
                  BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
                      return layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
                  } else {
                      return i_;
                  }
              }
              BOOST_UBLAS_INLINE
              size_type index2 () const {
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
                      return layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
                  } else {
                      return j_;
                  }
              }
  
              // Assignment
              BOOST_UBLAS_INLINE
              iterator1 &operator = (const iterator1 &it) {
                  container_reference<self_type>::assign (&it ());
                  rank_ = it.rank_;
                  i_ = it.i_;
                  j_ = it.j_;
                  itv_ = it.itv_;
                  it_ = it.it_;
                  return *this;
              }
  
              // Comparison
              BOOST_UBLAS_INLINE
              bool operator == (const iterator1 &it) const {
                  BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
                  // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
                  if (rank_ == 1 || it.rank_ == 1) {
                      return it_ == it.it_;
                  } else {
                      return i_ == it.i_ && j_ == it.j_;
                  }
              }
  
          private:
              int rank_;
              size_type i_;
              size_type j_;
              vector_subiterator_type itv_;
              subiterator_type it_;
  
              friend class const_iterator1;
          };
  
          BOOST_UBLAS_INLINE
          iterator1 begin1 () {
              return find1 (0, 0, 0);
          }
          BOOST_UBLAS_INLINE
          iterator1 end1 () {
              return find1 (0, size1_, 0);
          }
  
          class const_iterator2:
              public container_const_reference<coordinate_matrix>,
              public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
                                                 const_iterator2, value_type> {
          public:
              typedef typename coordinate_matrix::value_type value_type;
              typedef typename coordinate_matrix::difference_type difference_type;
              typedef typename coordinate_matrix::const_reference reference;
              typedef const typename coordinate_matrix::pointer pointer;
  
              typedef const_iterator1 dual_iterator_type;
              typedef const_reverse_iterator1 dual_reverse_iterator_type;
  
              // Construction and destruction
              BOOST_UBLAS_INLINE
              const_iterator2 ():
                  container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
              BOOST_UBLAS_INLINE
              const_iterator2 (const self_type &m, int rank, size_type i, size_type j, const vector_const_subiterator_type itv, const const_subiterator_type &it):
                  container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
              BOOST_UBLAS_INLINE
              const_iterator2 (const iterator2 &it):
                  container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {}
  
              // Arithmetic
              BOOST_UBLAS_INLINE
              const_iterator2 &operator ++ () {
                  if (rank_ == 1 && layout_type::fast_j ())
                      ++ it_;
                  else {
                      j_ = index2 () + 1;
                      if (rank_ == 1)
                          *this = (*this) ().find2 (rank_, i_, j_, 1);
                  }
                  return *this;
              }
              BOOST_UBLAS_INLINE
              const_iterator2 &operator -- () {
                  if (rank_ == 1 && layout_type::fast_j ())
                      -- it_;
                  else {
                      j_ = index2 () - 1;
                      if (rank_ == 1)
                          *this = (*this) ().find2 (rank_, i_, j_, -1);
                  }
                  return *this;
              }
  
              // Dereference
              BOOST_UBLAS_INLINE
              const_reference operator * () const {
                  BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                  BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
                  if (rank_ == 1) {
                      return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
                  } else {
                      return (*this) () (i_, j_);
                  }
              }
  
  #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_iterator1 begin () const {
                  const self_type &m = (*this) ();
                  return m.find1 (1, 0, index2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_iterator1 end () const {
                  const self_type &m = (*this) ();
                  return m.find1 (1, m.size1 (), index2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_reverse_iterator1 rbegin () const {
                  return const_reverse_iterator1 (end ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              const_reverse_iterator1 rend () const {
                  return const_reverse_iterator1 (begin ());
              }
  #endif
  
              // Indices
              BOOST_UBLAS_INLINE
              size_type index1 () const {
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
                      return layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
                  } else {
                      return i_;
                  }
              }
              BOOST_UBLAS_INLINE
              size_type index2 () const {
                  BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
                      return layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
                  } else {
                      return j_;
                  }
              }
  
              // Assignment
              BOOST_UBLAS_INLINE
              const_iterator2 &operator = (const const_iterator2 &it) {
                  container_const_reference<self_type>::assign (&it ());
                  rank_ = it.rank_;
                  i_ = it.i_;
                  j_ = it.j_;
                  itv_ = it.itv_;
                  it_ = it.it_;
                  return *this;
              }
  
              // Comparison
              BOOST_UBLAS_INLINE
              bool operator == (const const_iterator2 &it) const {
                  BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
                  // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
                  if (rank_ == 1 || it.rank_ == 1) {
                      return it_ == it.it_;
                  } else {
                      return i_ == it.i_ && j_ == it.j_;
                  }
              }
  
          private:
              int rank_;
              size_type i_;
              size_type j_;
              vector_const_subiterator_type itv_;
              const_subiterator_type it_;
          };
  
          BOOST_UBLAS_INLINE
          const_iterator2 begin2 () const {
              return find2 (0, 0, 0);
          }
          BOOST_UBLAS_INLINE
          const_iterator2 end2 () const {
              return find2 (0, 0, size2_);
          }
  
          class iterator2:
              public container_reference<coordinate_matrix>,
              public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
                                                 iterator2, value_type> {
          public:
              typedef typename coordinate_matrix::value_type value_type;
              typedef typename coordinate_matrix::difference_type difference_type;
              typedef typename coordinate_matrix::true_reference reference;
              typedef typename coordinate_matrix::pointer pointer;
  
              typedef iterator1 dual_iterator_type;
              typedef reverse_iterator1 dual_reverse_iterator_type;
  
              // Construction and destruction
              BOOST_UBLAS_INLINE
              iterator2 ():
                  container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
              BOOST_UBLAS_INLINE
              iterator2 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it):
                  container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
  
              // Arithmetic
              BOOST_UBLAS_INLINE
              iterator2 &operator ++ () {
                  if (rank_ == 1 && layout_type::fast_j ())
                      ++ it_;
                  else {
                      j_ = index2 () + 1;
                      if (rank_ == 1)
                          *this = (*this) ().find2 (rank_, i_, j_, 1);
                  }
                  return *this;
              }
              BOOST_UBLAS_INLINE
              iterator2 &operator -- () {
                  if (rank_ == 1 && layout_type::fast_j ())
                      -- it_;
                  else {
                      j_ = index2 ();
                      if (rank_ == 1)
                          *this = (*this) ().find2 (rank_, i_, j_, -1);
                  }
                  return *this;
              }
  
              // Dereference
              BOOST_UBLAS_INLINE
              reference operator * () const {
                  BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                  BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
                  if (rank_ == 1) {
                      return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
                  } else {
                      return (*this) ().at_element (i_, j_);
                  }
              }
  
  #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              iterator1 begin () const {
                  self_type &m = (*this) ();
                  return m.find1 (1, 0, index2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              iterator1 end () const {
                  self_type &m = (*this) ();
                  return m.find1 (1, m.size1 (), index2 ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              reverse_iterator1 rbegin () const {
                  return reverse_iterator1 (end ());
              }
              BOOST_UBLAS_INLINE
  #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
              typename self_type::
  #endif
              reverse_iterator1 rend () const {
                  return reverse_iterator1 (begin ());
              }
  #endif
  
              // Indices
              BOOST_UBLAS_INLINE
              size_type index1 () const {
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
                      return layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
                  } else {
                      return i_;
                  }
              }
              BOOST_UBLAS_INLINE
              size_type index2 () const {
                  BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
                  if (rank_ == 1) {
                      BOOST_UBLAS_CHECK (layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
                      return layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
                  } else {
                      return j_;
                  }
              }
  
              // Assignment
              BOOST_UBLAS_INLINE
              iterator2 &operator = (const iterator2 &it) {
                  container_reference<self_type>::assign (&it ());
                  rank_ = it.rank_;
                  i_ = it.i_;
                  j_ = it.j_;
                  itv_ = it.itv_;
                  it_ = it.it_;
                  return *this;
              }
  
              // Comparison
              BOOST_UBLAS_INLINE
              bool operator == (const iterator2 &it) const {
                  BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
                  // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
                  if (rank_ == 1 || it.rank_ == 1) {
                      return it_ == it.it_;
                  } else {
                      return i_ == it.i_ && j_ == it.j_;
                  }
              }
  
          private:
              int rank_;
              size_type i_;
              size_type j_;
              vector_subiterator_type itv_;
              subiterator_type it_;
  
              friend class const_iterator2;
          };
  
          BOOST_UBLAS_INLINE
          iterator2 begin2 () {
              return find2 (0, 0, 0);
          }
          BOOST_UBLAS_INLINE
          iterator2 end2 () {
              return find2 (0, 0, size2_);
          }
  
          // Reverse iterators
  
          BOOST_UBLAS_INLINE
          const_reverse_iterator1 rbegin1 () const {
              return const_reverse_iterator1 (end1 ());
          }
          BOOST_UBLAS_INLINE
          const_reverse_iterator1 rend1 () const {
              return const_reverse_iterator1 (begin1 ());
          }
  
          BOOST_UBLAS_INLINE
          reverse_iterator1 rbegin1 () {
              return reverse_iterator1 (end1 ());
          }
          BOOST_UBLAS_INLINE
          reverse_iterator1 rend1 () {
              return reverse_iterator1 (begin1 ());
          }
  
          BOOST_UBLAS_INLINE
          const_reverse_iterator2 rbegin2 () const {
              return const_reverse_iterator2 (end2 ());
          }
          BOOST_UBLAS_INLINE
          const_reverse_iterator2 rend2 () const {
              return const_reverse_iterator2 (begin2 ());
          }
  
          BOOST_UBLAS_INLINE
          reverse_iterator2 rbegin2 () {
              return reverse_iterator2 (end2 ());
          }
          BOOST_UBLAS_INLINE
          reverse_iterator2 rend2 () {
              return reverse_iterator2 (begin2 ());
          }
  
           // Serialization
          template<class Archive>
          void serialize(Archive & ar, const unsigned int /* file_version */){
              serialization::collection_size_type s1 (size1_);
              serialization::collection_size_type s2 (size2_);
              ar & serialization::make_nvp("size1",s1);
              ar & serialization::make_nvp("size2",s2);
              if (Archive::is_loading::value) {
                  size1_ = s1;
                  size2_ = s2;
              }
              ar & serialization::make_nvp("capacity", capacity_);
              ar & serialization::make_nvp("filled", filled_);
              ar & serialization::make_nvp("sorted_filled", sorted_filled_);
              ar & serialization::make_nvp("sorted", sorted_);
              ar & serialization::make_nvp("index1_data", index1_data_);
              ar & serialization::make_nvp("index2_data", index2_data_);
              ar & serialization::make_nvp("value_data", value_data_);
              storage_invariants();
          }
  
      private:
          void storage_invariants () const
          {
              BOOST_UBLAS_CHECK (capacity_ == index1_data_.size (), internal_logic ());
              BOOST_UBLAS_CHECK (capacity_ == index2_data_.size (), internal_logic ());
              BOOST_UBLAS_CHECK (capacity_ == value_data_.size (), internal_logic ());
              BOOST_UBLAS_CHECK (filled_ <= capacity_, internal_logic ());
              BOOST_UBLAS_CHECK (sorted_filled_ <= filled_, internal_logic ());
              BOOST_UBLAS_CHECK (sorted_ == (sorted_filled_ == filled_), internal_logic ());
          }
  
          size_type size1_;
          size_type size2_;
          array_size_type capacity_;
          mutable array_size_type filled_;
          mutable array_size_type sorted_filled_;
          mutable bool sorted_;
          mutable index_array_type index1_data_;
          mutable index_array_type index2_data_;
          mutable value_array_type value_data_;
          static const value_type zero_;
  
          BOOST_UBLAS_INLINE
          static size_type zero_based (size_type k_based_index) {
              return k_based_index - IB;
          }
          BOOST_UBLAS_INLINE
          static size_type k_based (size_type zero_based_index) {
              return zero_based_index + IB;
          }
  
          friend class iterator1;
          friend class iterator2;
          friend class const_iterator1;
          friend class const_iterator2;
      };
  
      template<class T, class L, std::size_t IB, class IA, class TA>
      const typename coordinate_matrix<T, L, IB, IA, TA>::value_type coordinate_matrix<T, L, IB, IA, TA>::zero_ = value_type/*zero*/();
  
  }}}
  
  #endif
  