AMDiS/lib/mtl4/boost/numeric/mtl/matrix/banded_view.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_MATRIX_BANDED_VIEW_INCLUDE
+#define MTL_MATRIX_BANDED_VIEW_INCLUDE
+
+#include <utility>
+#include <boost/shared_ptr.hpp>
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/mtl/utility/range_generator.hpp>
+#include <boost/numeric/mtl/utility/property_map.hpp>
+#include <boost/numeric/mtl/utility/parameters.hpp>
+#include <boost/numeric/mtl/matrix/crtp_base_matrix.hpp>
+#include <boost/numeric/mtl/matrix/base_matrix.hpp>
+#include <boost/numeric/mtl/operation/sfunctor.hpp>
+#include <boost/numeric/mtl/matrix/mat_expr.hpp>
+#include <boost/numeric/mtl/matrix/map_view.hpp>
+#include <boost/numeric/linear_algebra/identity.hpp>
+
+
+// Is not mutable because masking out some values forbids returning references
+//
+// Arbitrary combinations with other views (using shared_ptr) is planned
+
+namespace mtl { namespace matrix {
+
+// Forward
+namespace detail {
+    template <typename> struct banded_value;
+    template <typename, typename> struct mapped_value;
+}
+
+
+template <typename Matrix>
+struct banded_view
+  : public const_crtp_base_matrix< banded_view<Matrix>,
+				   typename Matrix::value_type, typename Matrix::size_type >,
+    public mat_expr< banded_view<Matrix> >,
+    public base_matrix<typename Matrix::value_type,
+		       typename mtl::traits::parameters<Matrix>::type>
+{
+    typedef banded_view                                self;
+    typedef mat_expr< self >                           expr_base;
+    typedef typename mtl::traits::parameters<Matrix>::type parameters;
+
+    typedef base_matrix<typename Matrix::value_type, parameters> base;
+
+    typedef Matrix                                     other;
+    typedef typename Matrix::orientation               orientation;
+    typedef typename Matrix::index_type                index_type;
+    // typedef typename Matrix::parameters                parameters;
+
+    typedef typename Matrix::value_type                value_type;
+    typedef typename Matrix::const_reference           const_reference;
+
+    typedef typename Matrix::key_type                  key_type;
+    typedef typename Matrix::size_type                 size_type;
+    typedef typename Matrix::dim_type                  dim_type;
+
+    typedef long int                                   bsize_type;
+
+    banded_view(const other& ref, bsize_type begin, bsize_type end)
+      : base(dim_type(mtl::matrix::num_rows(ref), mtl::matrix::num_cols(ref)), ref.nnz()),
+	ref(ref), begin(begin), end(end)
+    {}
+
+    banded_view(const boost::shared_ptr<Matrix>& p, bsize_type begin, bsize_type end)
+	: base(dim_type(mtl::matrix::num_rows(*p), mtl::matrix::num_cols(*p)), p->nnz()),
+	  my_copy(p), ref(*p), begin(begin), end(end)
+    {}
+
+#ifdef MTL_WITH_CPP11_MOVE
+    banded_view (self&& that) : my_copy(std::move(that.my_copy)), ref(that.ref), begin(that.begin), end(that.end) {}
+  banded_view (const self& that) : ref(that.ref), begin(that.begin), end(that.end) { assert(that.my_copy.use_count() == 0); }
+#endif
+
+    value_type operator() (size_type r, size_type c) const
+    {
+	using math::zero;
+	bsize_type bc= static_cast<bsize_type>(c), br= static_cast<bsize_type>(r),
+	           band= bc - br;
+	// Need value to return correct zero as well (i.e. matrices itself)
+	value_type v= ref(r, c);
+	return begin <= band && band < end ? v : zero(v);
+    }
+
+    // need const functions
+    bsize_type get_begin() const { return begin; }
+    bsize_type get_end() const { return end; }
+
+    template <typename> friend struct detail::banded_value;
+    template <typename, typename> friend struct detail::map_value;
+    //template <typename> friend struct ::mtl::sub_matrix_t<self>;
+
+    friend size_type inline num_rows(const self& A)
+    { 	using mtl::matrix::num_rows; return num_rows(A.ref);     }
+    friend size_type inline num_cols(const self& A)
+    { 	using mtl::matrix::num_cols; return num_cols(A.ref);     }
+
+  protected:
+    boost::shared_ptr<Matrix>           my_copy;
+  public:
+    const other&      ref;
+    bsize_type        begin, end;
+};
+
+template <typename Matrix>
+inline std::size_t size(const banded_view<Matrix>& A)
+{
+    return num_rows(A) * num_rows(A);
+}
+
+// ==========
+// Sub matrix
+// ==========
+
+template <typename Matrix>
+struct sub_matrix_t< mtl::matrix::banded_view<Matrix> >
+{
+    typedef mtl::matrix::banded_view<Matrix>                                           view_type;
+
+    // Mapping of sub-matrix type
+    typedef typename sub_matrix_t<Matrix>::sub_matrix_type                        ref_sub_type;
+    typedef mtl::matrix::banded_view<ref_sub_type>                                     const_sub_matrix_type;
+    typedef mtl::matrix::banded_view<ref_sub_type>                                     sub_matrix_type;
+    typedef typename view_type::size_type                                         size_type;
+
+    sub_matrix_type operator()(view_type const& view, size_type begin_r, size_type end_r,
+				     size_type begin_c, size_type end_c)
+    {
+	typedef boost::shared_ptr<ref_sub_type>                        pointer_type;
+
+	// Submatrix of referred matrix (or view)
+	// Create a submatrix, whos address will be kept by banded_view
+	pointer_type p(new ref_sub_type(sub_matrix(view.ref, begin_r, end_r, begin_c, end_c)));
+	return sub_matrix_type(p, view.begin, view.end);
+    }
+};
+
+
+}} // namespace mtl::matrix
+
+
+
+
+namespace mtl { namespace traits {
+
+    using mtl::matrix::banded_view;
+
+    template <typename Matrix>
+    struct row<banded_view<Matrix> >
+    {
+	// from map_view
+	typedef detail::mapped_row<sfunctor::identity<typename Matrix::value_type>, Matrix>   type;
+    };
+
+    template <typename Matrix>
+    struct col<banded_view<Matrix> >
+    {
+	// from map_view
+	typedef detail::mapped_col<sfunctor::identity<typename Matrix::value_type>, Matrix>   type;
+    };
+
+    namespace detail {
+
+	template <typename Matrix>
+	struct banded_value
+	{
+	    typedef typename Matrix::key_type              key_type;
+	    typedef typename Matrix::value_type            value_type;
+	    typedef banded_view<Matrix>                    view_type;
+
+	    banded_value(view_type const& view)
+		: view(view), its_row(view.ref), its_col(view.ref), its_value(view.ref)
+	    {}
+
+	    value_type operator() (key_type const& key) const
+	    {
+		using math::zero;
+		typedef typename view_type::bsize_type   bsize_type;
+
+		bsize_type br= static_cast<bsize_type>(its_row(key)),
+                           bc= static_cast<bsize_type>(its_col(key)),
+		           band= bc - br;
+		// Need value to return correct zero as well (i.e. matrices itself)
+		const value_type v= its_value(key);
+
+		return view.get_begin() <= band && band < view.get_end() ? v : zero(v);
+	    }
+
+	  protected:
+	    view_type const&                                view;
+	    typename row<Matrix>::type         its_row;
+	    typename col<Matrix>::type         its_col;
+	    typename const_value<Matrix>::type its_value;
+        };
+
+    } // detail
+
+    template <typename Matrix>
+    struct const_value<banded_view<Matrix> >
+    {
+	typedef detail::banded_value<Matrix>  type;
+    };
+
+    // ================
+    // Range generators
+    // ================
+
+    // Use range_generator of original matrix
+    template <typename Tag, typename Matrix>
+    struct range_generator<Tag, banded_view<Matrix> >
+	: public detail::referred_range_generator<banded_view<Matrix>, range_generator<Tag, Matrix> >
+    {};
+
+#if 0 // It is more complicated than this because referred_range_generator returns Matrix's cursor and we
+      // cannot dispatch on this anymore
+    template <typename Matrix>
+    struct range_generator<glas::tag::nz,
+			   typename detail::referred_range_generator<banded_view<Matrix>, range_generator<Tag, Matrix> >::type>
+
+			   detail::sub_matrix_cursor<banded_view<Matrix>, glas::tag::row, 2> >
+    {
+	typedef range_generator<glas::tag::row, banded_view<Matrix> >                    collection_type;
+	typedef range_generator<glas::tag::nz, range_generator<glas::tag::row, Matrix> > other_generator;
+	static int const                            level = other_generator::level;
+	typedef typename other_generator::type       type;
+
+	type begin(const collection_type& c)
+	{
+	    return
+
+	//
+	typedef typename range_generator<glas::tag::nz, detail::sub_matrix_cursor<Matrix>, glas::tag::row, 2>::type type;
+#endif
+
+
+    // To disambiguate
+    template <typename Matrix>
+    struct range_generator<tag::major, banded_view<Matrix> >
+	: public detail::referred_range_generator<banded_view<Matrix>, range_generator<tag::major, Matrix> >
+    {};
+
+
+
+
+
+
+}} // mtl::traits
+
+
+#endif // MTL_MATRIX_BANDED_VIEW_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/bands.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_MATRIX_BANDS_INCLUDE
+#define MTL_MATRIX_BANDS_INCLUDE
+
+#include <boost/numeric/mtl/matrix/banded_view.hpp>
+
+namespace mtl { namespace matrix {
+
+namespace traits {
+
+    template <typename Matrix>
+    struct bands
+    {
+	typedef banded_view<Matrix> type;
+    };
+}
+
+/// Returns a view of a matrix \p A from diagonal \p begin to \p end
+/** The main diagonal is numbered 0; the off-diagonal below the main one is -1.
+    Accordingly, the off-diagonal above the main is 1.
+    The parameters \p begin and \p end specify a right-open interval.
+    For, instance bands(A, -1, 2) yields a tridiagonal matrix. **/
+template <typename Matrix>
+typename traits::bands<Matrix>::type
+inline bands(const Matrix& A, long begin, long end)
+{
+    typedef typename traits::bands<Matrix>::type result;
+    return result(A, begin, end);
+}
+
+} // namespace matrix
+
+    using matrix::bands;
+
+} // namespace mtl
+
+#endif // MTL_MATRIX_BANDS_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/base_matrix.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_BASE_MATRIX_INCLUDE
+#define MTL_BASE_MATRIX_INCLUDE
+
+#include <algorithm>
+#include <boost/static_assert.hpp>
+#include <boost/mpl/bool.hpp>
+#include <boost/numeric/mtl/matrix/dimension.hpp>
+#include <boost/numeric/mtl/detail/index.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/is_static.hpp>
+
+namespace mtl { namespace matrix {
+
+/// Base class for other matrices, contains only very simple functionality that is used in all matrices.
+template <class Elt, class Parameters>
+struct base_matrix
+{
+    typedef base_matrix                       self;
+    typedef Elt                               value_type;
+    typedef typename Parameters::orientation  orientation;
+    typedef typename Parameters::index        index_type;
+    typedef typename Parameters::dimensions   dim_type;
+    static bool const                         on_stack= Parameters::on_stack;
+    typedef typename Parameters::size_type    size_type;
+  protected:
+    dim_type                        dim;       ///< # of rows and columns
+    size_type                       my_nnz;    ///< # of non-zeros, to be set by derived matrix
+    typedef mtl::traits::is_static<dim_type>       static_bool;
+
+  public:
+    base_matrix(size_type n= 0) :  my_nnz(n) {}
+
+    /// Setting dimension
+    explicit base_matrix(mtl::non_fixed::dimensions d, size_type n= 0) : dim(d), my_nnz(n) {}
+
+	///	Swap base matrix
+    friend void swap(self& x, self& y)
+    {
+	using std::swap;
+	swap(x.my_nnz, y.my_nnz);
+	swap(x.dim, y.dim);
+    }
+
+    /// Either matrix to be changed is uninitialized (i.e. 0x0) or dimensions are equal
+    /** The matrices with dimension 0 x 0 are considered like stem cells: they can still
+	change into an arbitrary dimension and are compatible with any other matrix.  Once a matrix has a non-trivial dimension
+	it can be only changed explicitly and is only compatible with matrices of the same dimensionality. **/
+    void check_dim(size_type MTL_DEBUG_ARG(num_rows), size_type MTL_DEBUG_ARG(num_cols)) const
+    {
+	MTL_DEBUG_THROW_IF(this->num_rows() * this->num_cols() != 0
+			   && (this->num_rows() != num_rows || this->num_cols() != num_cols),
+			   incompatible_size());
+    }
+
+#if 0
+    /** Will fail for fixed::dimension **/
+    void change_dim(mtl::non_fixed::dimensions d) { dim= d; }
+
+    template <std::size_t Rows, std::size_t Cols>
+    void change_dim(mtl::fixed::dimensions<Rows, Cols> d) {}
+#endif
+
+    void change_dim(size_type r, size_type c, boost::mpl::false_) { dim= dim_type(r, c); }
+    void change_dim(size_type r, size_type c, boost::mpl::true_) { check_dim(r, c); }
+
+    void change_dim(size_type r, size_type c) {	change_dim(r, c, static_bool()); }
+
+public:
+    /// Number of rows
+    size_type num_rows() const
+    {
+      return dim.num_rows();
+    }
+    /// First row taking indexing into account
+    size_type begin_row() const
+    {
+      return index::change_to(index_type(), 0);
+    }
+    /// Past-end row taking indexing into account
+    size_type end_row() const
+    {
+      return index::change_to(index_type(), num_rows());
+    }
+
+    /// number of colums
+    size_type num_cols() const
+    {
+      return dim.num_cols();
+    }
+    /// First column taking indexing into account
+    size_type begin_col() const
+    {
+      return index::change_to(index_type(), 0);
+    }
+    /// Past-end column taking indexing into account
+    size_type end_col() const
+    {
+      return index::change_to(index_type(), num_cols());
+    }
+
+    /// Number of non-zeros
+    size_type nnz() const
+    {
+      return my_nnz;
+    }
+
+  protected:
+    // dispatched functions for major dimension
+    size_type dim1(row_major) const
+    {
+      return num_rows();
+    }
+    size_type dim1(col_major) const
+    {
+      return num_cols();
+    }
+
+    // dispatched functions for minor dimension
+    size_type dim2(row_major) const
+    {
+      return num_cols();
+    }
+    size_type dim2(col_major) const
+    {
+      return num_rows();
+    }
+
+    // Dispatched functions for major
+    // Trailing _ due to conflicts with macro major
+    size_type major_(size_type r, size_type, row_major) const
+    {
+	return r;
+    }
+    size_type major_(size_type, size_type c, col_major) const
+    {
+	return c;
+    }
+
+  public:
+    /// Major dimension
+    size_type dim1() const
+    {
+      return dim1(orientation());
+    }
+
+    /// Minor dimension
+    size_type dim2() const
+    {
+      return dim2(orientation());
+    }
+
+    // Returns the row for row_major otherwise the column
+    // Trailing _ due to conflicts with macro major
+    size_type major_(size_type r, size_type c) const
+    {
+	return major_(r, c, orientation());
+    }
+
+    // Returns the row for col_major otherwise the column
+    // Trailing _ for consistency with major
+    size_type minor_(size_type r, size_type c) const
+    {
+	return major_(c, r, orientation());
+    }
+
+    // returns copy of dim
+    dim_type get_dimensions() const
+    {
+      return dim;
+    }
+};
+
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_BASE_MATRIX_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/base_sub_matrix.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_BASE_SUB_MATRIX_INCLUDE
+#define MTL_BASE_SUB_MATRIX_INCLUDE
+
+#include <algorithm>
+#include <boost/static_assert.hpp>
+#include <boost/numeric/mtl/matrix/dimension.hpp>
+#include <boost/numeric/mtl/detail/index.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+
+namespace mtl { namespace matrix {
+
+// Base class for sub-matrices
+// Contains only very simple functionality that is used in all sub-matrices
+// But also used in some complete matrices
+template <class Elt, class Parameters>
+struct base_sub_matrix
+{
+    typedef Elt                               value_type;
+    typedef typename Parameters::orientation  orientation;
+    typedef typename Parameters::index        index_type;
+    typedef typename Parameters::dimensions   dim_type;
+    static bool const                         on_stack= Parameters::on_stack;
+    typedef std::size_t                       size_type;
+    typedef base_sub_matrix                   self;
+
+  protected:
+    size_type                       my_nnz,       // # of non-zeros, to be set by derived matrix (drop maybe?)
+                                    my_begin_row, my_end_row,
+                                    my_begin_col, my_end_col;
+
+    void constructor_helper(const dim_type& dim)
+    {
+	my_begin_row= index::change_to(index_type(), 0);
+	my_end_row=   index::change_to(index_type(), dim.num_rows());
+	my_begin_col= index::change_to(index_type(), 0);
+	my_end_col=   index::change_to(index_type(), dim.num_cols());
+	my_nnz= 0;
+    }
+
+  public:
+    // base_sub_matrix() :  my_nnz(0), my_begin_row(0), my_end_row(0), my_begin_col(0), my_end_col(0) {}
+
+    base_sub_matrix()
+    {
+	// With no static dimension information, it is by default 0
+	constructor_helper(dim_type());
+    }
+
+    explicit base_sub_matrix(const dim_type& d)
+    //explicit base_sub_matrix(mtl::non_fixed::dimensions d)
+    {
+	constructor_helper(d);
+    }
+
+    friend void swap(self& x, self& y)
+    {
+	std::swap(x.my_nnz, y.my_nnz);
+	std::swap(x.my_begin_row, y.my_begin_row);
+	std::swap(x.my_end_row, y.my_end_row);
+	std::swap(x.my_begin_col, y.my_begin_col);
+	std::swap(x.my_end_col, y.my_end_col);
+    }
+
+    // Either changed matrix is uninitialized (i.e. 0x0) or dimensions are equal
+    void check_dim(size_type MTL_DEBUG_ARG(num_rows), size_type MTL_DEBUG_ARG(num_cols) ) const
+    {
+	MTL_DEBUG_THROW_IF(this->num_rows() * this->num_cols() != 0
+			   && (this->num_rows() != num_rows || this->num_cols() != num_cols),
+			   incompatible_size());
+    }
+
+protected:
+    void change_dim(non_fixed::dimensions d) { constructor_helper(d); }
+    template <std::size_t Rows, std::size_t Cols>
+    void change_dim(fixed::dimensions<Rows, Cols> d) { check_dim(d.num_rows(), d.num_cols()); }
+
+
+
+    void change_dim(size_type r, size_type c) {	change_dim(dim_type(r, c));  }
+
+    void set_ranges(size_type br, size_type er, size_type bc, size_type ec)
+    {
+	MTL_DEBUG_THROW_IF(br > er, range_error("begin row > end row"));
+	MTL_DEBUG_THROW_IF(bc > ec, range_error("begin column > end column"));
+	my_begin_row= br; my_end_row= er; my_begin_col= bc; my_end_col= ec;
+    }
+
+public:
+    void check_ranges(size_type MTL_DEBUG_ARG(begin_r), size_type MTL_DEBUG_ARG(end_r),
+		      size_type MTL_DEBUG_ARG(begin_c), size_type MTL_DEBUG_ARG(end_c) ) const
+    {
+	MTL_DEBUG_THROW_IF(begin_r < begin_row(), range_error("begin_row out of range"));
+	// if (end_r > end_row()) std::cout << "end_row out of range\n";
+	MTL_DEBUG_THROW_IF(end_r > end_row(), range_error("end_row out of range"));
+
+	MTL_DEBUG_THROW_IF(begin_c < begin_col(), range_error("begin_col out of range"));
+	MTL_DEBUG_THROW_IF(end_c > end_col(), range_error("end_col out of range"));
+    }
+
+    explicit base_sub_matrix(size_type br, size_type er, size_type bc, size_type ec) : my_nnz(0)
+    {
+	set_ranges(br, er, bc, ec);
+    }
+
+    // Number of rows
+    size_type num_rows() const
+    {
+	return my_end_row - my_begin_row;
+    }
+
+    // First row taking indexing into account (already stored as such)
+    size_type begin_row() const
+    {
+	return my_begin_row;
+    }
+
+    // Past-end row taking indexing into account (already stored as such)
+    size_type end_row() const
+    {
+	return my_end_row;
+    }
+
+    // Number of columns
+    size_type num_cols() const
+    {
+	return my_end_col - my_begin_col;
+    }
+
+    // First column taking indexing into account (already stored as such)
+    size_type begin_col() const
+    {
+	return my_begin_col;
+    }
+
+    // Past-end column taking indexing into account (already stored as such)
+    size_type end_col() const
+    {
+	return my_end_col;
+    }
+
+    // Number of non-zeros
+    size_type nnz() const
+    {
+      return my_nnz;
+    }
+
+  protected:
+    // dispatched functions for major dimension
+    size_type dim1(row_major) const
+    {
+	return num_rows();
+    }
+
+    size_type dim1(col_major) const
+    {
+	return num_cols();
+    }
+
+    // dispatched functions for minor dimension
+    size_type dim2(row_major) const
+    {
+	return num_cols();
+    }
+
+    size_type dim2(col_major) const
+    {
+	return num_rows();
+    }
+
+    // Dispatched functions for major
+    // Trailing _ due to conflicts with macro major
+    size_type major_(size_type r, size_type, row_major) const
+    {
+	return r;
+    }
+
+    size_type major_(size_type, size_type c, col_major) const
+    {
+	return c;
+    }
+
+  public:
+    // return major dimension
+    size_type dim1() const
+    {
+	return dim1(orientation());
+    }
+
+    // return minor dimension
+    size_type dim2() const
+    {
+	return dim2(orientation());
+    }
+
+    // Returns the row for row_major otherwise the column
+    // Trailing _ due to conflicts with macro major
+    size_type major_(size_type r, size_type c) const
+    {
+	return major_(r, c, orientation());
+    }
+
+    // Returns the row for col_major otherwise the column
+    // Trailing _ for consistency with major
+    size_type minor_(size_type r, size_type c) const
+    {
+	return major_(c, r, orientation());
+    }
+
+    dim_type get_dimensions() const
+    {
+	return dim_type(num_rows(), num_cols());
+    }
+
+};
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_BASE_SUB_MATRIX_INCLUDE
+
+
+/*
+   Question:
+   - Shall we keep the position in the original matrix?
+*/

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/block_diagonal2D.hpp

+
+#ifndef MTL_DIST_BLOCK_DIAGONAL2D_INCLUDE
+#define MTL_DIST_BLOCK_DIAGONAL2D_INCLUDE
+
+
+
+#include <vector>
+#include <cassert>
+#include <limits>
+#include <boost/type_traits/is_same.hpp>
+#include <boost/mpl/bool.hpp>
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/irange.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+#include <boost/numeric/mtl/matrix/mat_expr.hpp>
+
+
+namespace mtl {
+
+  namespace matrix {
+
+
+/// Block diagonal matrix structure
+/** Blocks can be any existing matrix_type in mtl4. **/
+template <typename Matrix>
+class block_diagonal2D
+  : public mat_expr< block_diagonal2D<Matrix> >
+{
+  public:
+
+    typedef unsigned int	      size_type;
+    typedef Matrix		      block_type;
+    typedef std::vector<Matrix>       block_matrix_type;
+    typedef std::vector<size_type>    index_vector_type;
+    typedef block_diagonal2D<Matrix>  self;
+    typedef typename Collection<Matrix>::value_type value_type;
+
+    /// Constructor: number or rows and columns and optionally estimated number of blocks
+    explicit block_diagonal2D(size_type rows, size_type cols, size_type init_size= 0)
+      : nrows(rows), ncols(cols), nb_blocks(0), my_nnz(0),
+	min_ind(std::numeric_limits<size_type>::max()), max_ind(std::numeric_limits<size_type>::min()),
+	compact_heap(0)
+    {
+	start_block.reserve(init_size);
+	end_block.reserve(init_size);
+	blocks.reserve(init_size);
+    }
+
+    ~block_diagonal2D() { delete[] compact_heap; }
+
+    /// Returns the global number of rows
+    size_type num_rows() const {
+	return nrows ;
+    }
+    /// Returns the global number of coluums
+    size_type num_cols() const {
+	return ncols ;
+    }
+    /// Returns the global number of blocks
+    size_type num_blocks() const {
+	return  nb_blocks;
+    }
+
+    /// Minimal index
+    size_type min_index() const { return min_ind; }
+
+    /// Maximal index
+    size_type max_index() const { return max_ind; }
+
+
+    block_type const& block(size_type i) const
+    {
+	MTL_DEBUG_THROW_IF(is_negative(i) || i >= nb_blocks, index_out_of_range());
+	return blocks[i];
+    }
+
+    /// Insert a block from start x start to end x end
+    void insert(size_type start, size_type end, const block_type& A)
+    {
+	MTL_DEBUG_THROW_IF(start > end, logic_error());
+	MTL_DEBUG_THROW_IF(is_negative(start) || end > nrows || end > ncols, index_out_of_range());
+	assert(compact_heap == 0); // insertion after make_compact; might be relaxed later
+
+	start_block.push_back(start);
+	end_block.push_back(end);
+	blocks.push_back(A);
+	my_nnz+= A.nnz();
+	nb_blocks++;
+
+	if (start < min_ind)
+	    min_ind= start;
+	if (end > max_ind)
+	    max_ind= end;
+#if 0
+	boost::mpi::communicator world;
+	if (world.rank() == 0) {
+	    std::cout << "block_diag.insert " << nb_blocks-1 << "th block: start == "
+		      << start_block << ", end == " << end_block << "block is:\n"
+		      << blocks.back();
+	}
+#endif
+    }
+
+    /// Element A[i][j] by summing over all blocks, use only for debugging because it is very slow
+    value_type operator()(size_type i, size_type j) const
+    {
+	value_type s= value_type(0);
+	for (std::size_t b= 0; b < blocks.size(); ++b) {
+	    size_type st= start_block[b], e= end_block[b];
+	    if (st <= i && st <= j && i < e && j < e)
+		s+= blocks[b][i - st][j - st];
+	}
+	return s;
+    }
+
+    /// Memory of inserted
+    void make_compact(boost::mpl::true_)
+    {
+	assert(compact_heap == 0); // might be relaxed later
+
+	size_type entries= 0;
+	for (size_type i= 0; i < nb_blocks; i++)
+	    entries+= size(blocks[i]);
+	compact_heap= new value_type[entries];
+
+	size_type pos= 0;
+	for (size_type i= 0; i < nb_blocks; i++) {
+	    size_type s= end_block[i] - start_block[i];
+	    block_type tmp(s, s, compact_heap + pos);
+	    tmp= blocks[i];
+	    pos+= size(blocks[i]);
+	    swap(blocks[i], tmp);
+	}
+	assert(pos == entries);
+    }
+    void make_compact(boost::mpl::false_) {}
+
+    void make_compact() { make_compact(boost::is_same<typename mtl::traits::category<block_type>::type, tag::dense2D>()); }
+
+    /// Number of non-zeros (accumulated over blocks)
+    size_type nnz() const { return my_nnz; }
+
+    template <typename VectorIn, typename VectorOut>
+    void add_mult(const VectorIn& x, VectorOut& y) const
+    {
+	mtl::vampir_trace<3062> tracer;
+	MTL_DEBUG_THROW_IF(ncols != size(x) || nrows != size(y), incompatible_size());
+
+	// set_to_zero(y);
+	for(size_type i= 0; i < nb_blocks; i++) {
+	    mtl::irange r(start_block[i], end_block[i]);
+	    y[r]+= blocks[i] * x[r];
+	}
+    }
+
+    ///block_diagonal-matrix  times cvec
+    template <typename VectorIn, typename VectorOut>
+    void mult(const VectorIn& x, VectorOut& y) const
+    {
+	mtl::vampir_trace<3062> tracer;
+	MTL_DEBUG_THROW_IF(ncols != size(x) || nrows != size(y), incompatible_size());
+
+	set_to_zero(y);
+	for(size_type i= 0; i < nb_blocks; i++) {
+	    mtl::irange r(start_block[i], end_block[i]);
+	    y[r]= blocks[i] * x[r];
+	}
+    }
+
+    template <typename VectorIn>
+    struct multiplier
+      : mtl::vector::assigner<multiplier<VectorIn> >
+    {
+	explicit multiplier(const self& P, const VectorIn& x) : P(P), x(x) {}
+
+	template <typename VectorOut>
+	void assign_to(VectorOut& y) const
+	{   P.mult(x, y); }
+
+	const self& P;
+	const VectorIn& x;
+    };
+
+
+    template <typename VectorIn>
+    multiplier<VectorIn> operator*(const VectorIn& x) // const
+    {  return multiplier<VectorIn>(*this, x); }
+
+
+  private:
+    size_type           nrows, ncols, nb_blocks, my_nnz, min_ind, max_ind;
+    index_vector_type  	start_block, end_block;
+    block_matrix_type   blocks;
+    value_type*         compact_heap;
+};
+
+// ================
+// Free functions
+// ================
+
+
+/// Number of rows
+template <typename Value>
+typename block_diagonal2D<Value>::size_type
+inline num_rows(const block_diagonal2D<Value>& matrix)
+{
+    return matrix.num_rows();
+}
+
+/// Number of columns
+template <typename Value>
+typename block_diagonal2D<Value>::size_type
+inline num_cols(const block_diagonal2D<Value>& matrix)
+{
+    return matrix.num_cols();
+}
+
+/// Size of the matrix, i.e. the number of row times columns
+template <typename Value>
+typename block_diagonal2D<Value>::size_type
+inline size(const block_diagonal2D<Value>& matrix)
+{
+    return matrix.num_cols() * matrix.num_rows();
+}
+
+// /// Distributed Vector of start block entrys in the matrix
+// template <typename Value>
+// typename block_diagonal2D<Value>::index_type
+// inline start(const block_diagonal2D<Value>& matrix)
+// {
+//     return matrix.start();
+// }
+
+// /// Distributed Vector of en block entrys in the matrix
+// template <typename Value>
+// typename block_diagonal2D<Value>::index_type
+// inline end(const block_diagonal2D<Value>& matrix)
+// {
+//     return matrix.end();
+// }
+
+
+
+
+} // namespace matrix
+} // namespace mtl
+
+#endif // MTL_DIST_BLOCK_DIAGONAL2D_INCLUDE
+

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/compressed2D.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_COMPRESSED2D_INCLUDE
+#define MTL_COMPRESSED2D_INCLUDE
+
+#include <algorithm>
+#include <vector>
+#include <map>
+#include <cmath>
+#include <boost/tuple/tuple.hpp>
+#include <boost/type_traits/is_same.hpp>
+#include <boost/mpl/if.hpp>
+
+#include <boost/numeric/linear_algebra/identity.hpp>
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/config.hpp>
+#include <boost/numeric/mtl/utility/common_include.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/maybe.hpp>
+#include <boost/numeric/mtl/utility/shrink_stl_vector.hpp>
+#include <boost/numeric/mtl/utility/zipped_sort.hpp>
+#include <boost/numeric/mtl/detail/base_cursor.hpp>
+#include <boost/numeric/mtl/operation/is_negative.hpp>
+#include <boost/numeric/mtl/operation/update.hpp>
+#include <boost/numeric/mtl/operation/shift_block.hpp>
+#include <boost/numeric/mtl/matrix/parameter.hpp>
+#include <boost/numeric/mtl/matrix/base_matrix.hpp>
+#include <boost/numeric/mtl/matrix/crtp_base_matrix.hpp>
+#include <boost/numeric/mtl/matrix/mat_expr.hpp>
+#include <boost/numeric/mtl/matrix/element_matrix.hpp>
+#include <boost/numeric/mtl/matrix/element_array.hpp>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+#include <boost/numeric/mtl/operation/compute_factors.hpp>
+#include <boost/numeric/mtl/operation/num_rows.hpp>
+#include <boost/numeric/mtl/operation/num_cols.hpp>
+#include <boost/numeric/mtl/operation/size.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+#ifdef MTL_WITH_INITLIST
+# include <initializer_list>
+#endif
+
+#undef major // fight namespace pollution
+
+namespace mtl { namespace matrix {
+
+
+struct compressed_key
+{
+    typedef std::size_t                               size_t;
+    typedef compressed_key                            self;
+
+    template <typename Elt, typename Parameters>
+    explicit compressed_key(compressed2D<Elt, Parameters> const& matrix, size_t offset) : offset(offset)
+    {
+	std::size_t my_major= matrix.indexer.find_major(matrix, offset);
+	major= my_major;
+    }
+
+    template <typename Elt, typename Parameters>
+    explicit compressed_key(compressed2D<Elt, Parameters> const& matrix, size_t r, size_t c)
+    {
+	offset= matrix.indexer(matrix, r, c).value();
+	major= matrix.indexer.major_minor_c(matrix, r, c).first;
+    }
+
+    compressed_key(compressed_key const& other) { offset= other.offset; major= other.major; }
+
+    self& operator= (self const& other)
+    {
+	offset= other.offset; major= other.major;
+	return *this;
+    }
+
+    bool operator== (compressed_key const& other) const
+    {
+	//if (offset == other.offset && major != other.major)
+	//    std::cout << offset << " " << other.offset << " " << major << " " << other.major << '\n';
+	// The following tests doesn't hold everywhere (anymore)
+	// MTL_DEBUG_THROW_IF(offset == other.offset && major != other.major, logic_error("equal offsets imply equal major"));
+	return offset == other.offset;
+    }
+
+    bool operator!= (compressed_key const& other) const { return !(*this == other); }
+
+    size_t       major;
+    size_t       offset;
+};
+
+
+// Cursor over every element
+template <typename Elt, typename Parameters>
+struct compressed_el_cursor
+    : public compressed_key
+{
+    typedef Elt                           value_type;
+    typedef compressed_key                base;
+    typedef compressed_el_cursor          self;
+    typedef std::size_t                   size_t;
+
+    explicit compressed_el_cursor(compressed2D<Elt, Parameters> const& matrix, size_t r, size_t c)
+	: base(matrix, r, c), matrix(matrix)    {}
+
+    explicit compressed_el_cursor(compressed2D<Elt, Parameters> const& matrix, size_t offset)
+	: base(matrix, offset), matrix(matrix)    {}
+
+    compressed_el_cursor(const compressed_el_cursor<Elt, Parameters>& other)
+	: base(other), matrix(other.matrix)     {}
+
+    self& operator= (self const& other) { base::operator=(other); return *this; }
+
+    self& operator++ ()
+    {
+	++offset;
+	MTL_DEBUG_THROW_IF(matrix.starts[major+1] < offset, runtime_error("Inconsistent incrementation!"));
+	while (major < matrix.starts.size()-1 && matrix.starts[major+1] == offset)
+	    ++major;
+	return *this;
+    }
+
+    base& operator* () { return *this; }
+    const base& operator* () const { return *this; }
+
+    compressed2D<Elt, Parameters> const& matrix;
+};
+
+
+// Cursor over every element
+template <typename Elt, typename Parameters>
+struct compressed_minor_cursor
+    : public compressed_key
+{
+    typedef Elt                           value_type;
+    typedef compressed_key                base;
+    typedef compressed_minor_cursor       self;
+    typedef std::size_t                   size_t;
+
+    explicit compressed_minor_cursor(mtl::matrix::compressed2D<Elt, Parameters> const& matrix, size_t r, size_t c)
+	: base(matrix, r, c), matrix(matrix)
+    {}
+
+    explicit compressed_minor_cursor(mtl::matrix::compressed2D<Elt, Parameters> const& matrix, size_t offset)
+	: base(matrix, offset), matrix(matrix)
+    {}
+
+    compressed_minor_cursor(self const& other) : base(other), matrix(other.matrix) {}
+
+    self& operator= (self const& other)
+    {	base::operator=(other);	return *this; }
+
+    self& operator++() { ++offset; return *this;  }
+    self& operator+=(size_t inc) { offset+= inc; return *this; }
+    self operator+(size_t inc) const  {	self tmp(*this); tmp+= inc; return tmp; }
+
+    self& operator--() { --offset; return *this; }
+    base& operator* () { return *this; }
+    const base& operator* () const { return *this; }
+
+    mtl::matrix::compressed2D<Elt, Parameters> const& matrix;
+};
+
+
+
+
+// Indexing for compressed matrices
+template <typename SizeType>
+struct compressed2D_indexer
+{
+    typedef SizeType                          size_type;
+    typedef std::pair<size_type, size_type>   size_pair;
+  private:
+    // helpers for public functions
+    template <class Matrix>
+    utilities::maybe<size_type> offset(const Matrix& ma, size_type major, size_type minor) const
+    {
+	typedef utilities::maybe<size_type>      result_type;
+	assert(ma.starts[major] <= ma.starts[major+1]); // Check sortedness
+	assert(ma.starts[major+1] <= ma.my_nnz);        // Check bounds of indices
+	// Now we are save to use past-end addresses as iterators
+
+	// Empty matrices are special cases
+	if (ma.indices.empty())
+	    return result_type(0, false);
+
+	const size_type *first = &ma.indices[0] + ma.starts[major],
+	             *last = &ma.indices[0] + ma.starts[major+1];
+	// if empty row (or column) return start of next one
+	if (first == last)
+	    return result_type(first - &ma.indices[0], false);
+
+	const size_type *index= first;
+	if (last - index <= int(compressed_linear_search_limit))
+	    while (index != last && *index < minor) ++index;
+	else
+	    index = std::lower_bound(first, last, minor);
+	return result_type(index - &ma.indices[0], index != last && *index == minor);
+    }
+
+  public:
+    // Returns major and minor index in C style (starting with 0)
+    template <class Matrix>
+    size_pair major_minor_c(const Matrix& ma, size_type row, size_type col) const
+    {
+	using std::make_pair;
+	// convert into c indices
+	typename Matrix::index_type my_index;
+	size_type my_row= index::change_from(my_index, row),
+	       my_col= index::change_from(my_index, col);
+	return make_pair(ma.major_(my_row, my_col), ma.minor_(my_row, my_col));
+    }
+
+    // Returns the offset if found
+    // If not found it returns the position where it would be inserted
+    template <class Matrix>
+    utilities::maybe<size_type> operator() (const Matrix& ma, size_type row, size_type col) const
+    {
+	size_type major, minor;
+	boost::tie(major, minor) = major_minor_c(ma, row, col);
+	return offset(ma, major, minor);
+    }
+
+    // Same as above if internal representation is already known
+    template <class Matrix>
+    utilities::maybe<size_type> operator() (const Matrix& ma, size_pair major_minor) const
+    {
+	return offset(ma, major_minor.first, major_minor.second);
+    }
+
+    // For a given offset the minor can be accessed directly, the major dim has to be searched
+    // Returned in internal (c) representation
+    template <class Matrix>
+    size_type find_major(const Matrix& ma, size_type offset) const
+    {
+	MTL_DEBUG_THROW_IF(ma.starts.empty(), logic_error("Major vector can't be empty"));
+	size_type my_major= std::upper_bound(ma.starts.begin(), ma.starts.end(), offset) - ma.starts.begin();
+	return --my_major;
+    }
+
+    template <class Matrix>
+    size_type minor_from_offset(const Matrix& ma, size_type offset) const
+    {
+	typedef typename Matrix::index_type my_index;
+	return index::change_to(my_index(), ma.indices[offset]);
+    }
+
+}; // compressed2D_indexer
+
+
+/// Compressed 2D matrix type
+// For now no external data
+template <typename Elt, typename Parameters = matrix::parameters<> >
+class compressed2D
+  : public base_matrix<Elt, Parameters>,
+    public const_crtp_base_matrix< compressed2D<Elt, Parameters>, Elt, typename Parameters::size_type >,
+    public crtp_matrix_assign< compressed2D<Elt, Parameters>, Elt, typename Parameters::size_type >,
+    public mat_expr< compressed2D<Elt, Parameters> >
+{
+    typedef std::size_t                              size_t;
+    typedef base_matrix<Elt, Parameters>             super;
+    typedef compressed2D                             self;
+    typedef mat_expr< compressed2D<Elt, Parameters> >          expr_base;
+
+    // Only allocation of new data, doesn't copy if already existent
+    void allocate(size_t new_nnz)
+    {
+	if (new_nnz) {
+	    this->my_nnz = new_nnz;
+	    data.resize(this->my_nnz);
+	    indices.resize(this->my_nnz, 0);
+	}
+    }
+
+  public:
+    typedef Parameters                               parameters;
+    typedef typename Parameters::orientation         orientation;
+    typedef typename Parameters::index               index_type;
+    typedef typename Parameters::dimensions          dimensions;
+    typedef Elt                                      value_type;
+    typedef compressed_key                           key_type;
+    // const value_type& isn't defined everywhere
+    typedef value_type                               const_reference;
+
+    typedef typename Parameters::size_type           size_type;
+    typedef crtp_matrix_assign<self, Elt, size_type>  assign_base;
+    typedef compressed2D_indexer<size_type>          indexer_type;
+
+    void check() const { MTL_DEBUG_THROW_IF(inserting, access_during_insertion()); }
+
+    /// Removes all values; e.g. for set_to_zero
+    void make_empty()
+    {
+	check();
+	this->my_nnz = 0;
+	data.resize(0);
+	indices.resize(0);
+	std::fill(starts.begin(), starts.end(), 0);
+    }
+
+    /// Change dimension of the matrix; data get lost.
+    void change_dim(size_type r, size_type c)
+    {
+	check();
+	if (this->num_rows() != r || this->num_cols() != c) {
+	    super::change_dim(r, c);
+	    starts.resize(this->dim1()+1);
+	    make_empty();
+	}
+    }
+
+    /// Sets nnz and resizes indices and data
+    void set_nnz(size_type n)
+    {
+	check();
+	indices.resize(n);
+	data.resize(n);
+	this->my_nnz= n;
+    }
+
+    // if compile time matrix size, we can set the start vector
+    /// Default constructor
+    explicit compressed2D () : inserting(false)
+    {
+	if (super::dim_type::is_static) starts.resize(super::dim1() + 1);
+    }
+
+
+    /// Setting dimension and allocate starting vector
+    explicit compressed2D (mtl::non_fixed::dimensions d, size_t nnz = 0)
+      : super(d), inserting(false)
+    {
+	starts.resize(super::dim1() + 1, 0);
+	allocate(nnz);
+    }
+
+    /// Setting dimension and allocate starting vector
+    explicit compressed2D (size_type num_rows, size_type num_cols, size_t nnz = 0)
+      : super(non_fixed::dimensions(num_rows, num_cols)), inserting(false)
+    {
+	starts.resize(super::dim1() + 1, 0);
+	allocate(nnz);
+    }
+
+    /// Initializing matrix from 3 arrays: \p Major, \p Minor,\p Entries.
+    /** For a row-major matrix Major shall contain the ranges of rows and Minor the column-indices.
+	This function is for advanced users only and should be used with care:
+	errors can lead to inconsistent matrices and segmentation faults.
+     **/
+    compressed2D (size_type m, size_type n, size_type nnz,
+		  size_type* Major, size_type* Minor, value_type* Entries)
+      : super(non_fixed::dimensions(m, n)), inserting(false)
+    {
+      allocate(nnz);
+
+      data.assign(Entries, Entries + nnz);
+      starts.assign(Major, Major + super::dim1() + 1);
+      indices.assign(Minor, Minor + nnz);
+    }
+
+
+#if 0 // Default is faster !!!
+    /// Copy constructor (just in case that is not generated by all compilers (generic matrix copy slower))
+    compressed2D(const self& src)
+      : super(non_fixed::dimensions(src.num_rows(), src.num_cols())), data(src.data),
+	starts(src.starts), indices(src.indices), inserting(false)
+    {}
+#endif
+
+    /// Copy from other types
+    template <typename MatrixSrc>
+    explicit compressed2D (const MatrixSrc& src) : inserting(false)
+    {
+	if (super::dim_type::is_static) starts.resize(super::dim1() + 1);
+	*this= src;
+    }
+
+
+#if defined(MTL_WITH_INITLIST) && defined(MTL_WITH_AUTO) && defined(MTL_WITH_RANGEDFOR)
+    /// Constructor for initializer list \p values
+    template <typename Value2>
+    compressed2D(std::initializer_list<std::initializer_list<Value2> > values)
+      : super(mtl::non_fixed::dimensions(values.size(), values.size()? values.begin()->size() : 0)),
+	inserting(false)
+    {
+	starts.resize(super::dim1() + 1, 0);
+	allocate(super::dim1() * super::dim2());
+	*this= values;
+    }
+#endif
+
+#ifdef MTL_WITH_DEFAULTIMPL
+    compressed2D(const compressed2D&) = default;
+#endif
+
+#ifdef MTL_WITH_MOVE
+    self& operator=(self&& src)
+    {
+	std::cout << "I am moving !!!\n";
+	check();
+	this->checked_change_dim(src.num_rows(), src.num_cols());
+	swap(*this, src);
+	return *this;
+    }
+
+    self& operator=(const self& src)
+    {
+	if (this == &src)
+	    return *this;
+	check();
+	this->checked_change_dim(src.num_rows(), src.num_cols());
+	set_nnz(src.nnz());
+
+	starts= src.starts;
+	indices= src.indices;
+	data= src.data;
+	return *this;
+    }
+#else
+    /// Consuming assignment operator (without move semantics)
+    self& operator=(self src)
+    {
+	assert(this != &src);  // Self-copy would be an indication of an error
+
+	check();
+	// TODO Check the following line
+	this->checked_change_dim(src.num_rows(), src.num_cols());
+	swap(*this, src);
+	return *this;
+    }
+#endif
+
+    using assign_base::operator=;
+
+#if 0
+    void self_copy(const self& src)
+    {
+	if (this == &src)
+	    return;
+	change_dim(num_rows(src), num_cols(src));
+	set_nnz(src.nnz());
+
+	#ifdef MTL_WITH_OPENMP
+        #  pragma omp parallel for schedule(static, 16)
+	#endif
+	for (int i= 0; i < int(src.dim1()); i++) {
+	    const size_type cj0= src.starts[i], cj1= src.starts[i+1];
+	    starts[i]= cj0;
+	    for (size_type j0= cj0; j0 != cj1; ++j0) {
+		indices[j0]= src.indices[j0];
+		data[j0]= src.data[j0];
+	    }
+	}
+	starts[src.dim1()]= src.starts[src.dim1()];
+    }
+#endif
+
+    /// Copy the pattern of matrix \p src
+    /** Advanced feature. Not generic, i.e. does not exist for all matrix types.
+	value_type of \p src and target can be different.
+	Changes dimension of target to that of \p src.
+	Existing entries are deleted. **/
+    template <typename Src>
+    void copy_pattern(const Src& src)
+    {
+	using std::copy;
+	change_dim(num_rows(src), num_cols(src));
+	copy(src.ref_major().begin(), src.ref_major().end(), starts.begin());
+
+	set_nnz(src.ref_minor().size());
+	copy(src.ref_minor().begin(), src.ref_minor().end(), indices.begin());
+    }
+
+    void make_symmetric_pattern()
+    {
+	MTL_THROW_IF(this->num_cols() != this->num_rows(), matrix_not_square());
+	self A(*this);
+	{
+	    using math::zero;
+	    inserter<self> ins(A, 2 * this->nnz() / this->dim1());
+
+	    typename traits::row<self>::type                                 row(*this);
+	    typename traits::col<self>::type                                 col(*this);
+	    typename traits::const_value<self>::type                         value(*this);
+	    typedef typename traits::range_generator<mtl::tag::major, self>::type        cursor_type;
+	    for (cursor_type cursor = mtl::begin<mtl::tag::major>(*this), cend = mtl::end<mtl::tag::major>(*this);
+		 cursor != cend; ++cursor) {
+		typedef mtl::tag::nz     inner_tag;
+		typedef typename traits::range_generator<inner_tag, cursor_type>::type icursor_type;
+		for (icursor_type icursor = mtl::begin<inner_tag>(cursor), icend = mtl::end<inner_tag>(cursor); icursor != icend; ++icursor) {
+		    size_type r= row(*icursor), c= col(*icursor);
+		    if (!indexer(*this, c, r))
+			ins[c][r] << zero(value(*icursor));
+	      	}
+	    }
+	}
+	swap(*this, A);
+    }
+
+
+    // Copies range of values and their coordinates into compressed matrix
+    // For brute force initialization, should be used with uttermost care
+    // Won't be suitable for distributed matrices, take care of this to this later
+    template <typename ValueIterator, typename StartIterator, typename IndexIterator>
+    void raw_copy(ValueIterator first_value, ValueIterator last_value,
+		  StartIterator first_start, IndexIterator first_index)
+    {
+	using std::copy;
+
+	// check if starts has right size
+	allocate(last_value - first_value); // ????
+	// check if nnz and indices has right size
+
+	copy(first_value, last_value, data.begin());
+	copy(first_start, first_start + this->dim1() + 1, starts.begin());
+	copy(first_index, first_index + this->nnz(), indices.begin());
+    }
+
+    /// Value of matrix entry
+    const_reference operator() (size_type row, size_type col) const
+    {
+	using math::zero;
+	check(); MTL_DEBUG_THROW_IF(is_negative(row) || row >= this->num_rows() || is_negative(col) || col >= this->num_cols(), index_out_of_range());
+	utilities::maybe<size_type> pos = indexer(*this, row, col);
+	return pos ? data[pos.value()] : zero(value_type());
+    }
+
+    /// L-value reference of stored matrix entry
+    /** To be used with care; in debug mode, exception is thrown if entry is not found **/
+    value_type& lvalue(size_type row, size_type col)
+    {
+	utilities::maybe<size_type> pos = indexer(*this, row, col);
+	check(); MTL_DEBUG_THROW_IF(!pos, logic_error("This entry does not exist in the matrix"));
+	return data[pos.value()];
+    }
+
+    // For internal use
+    const value_type& value_from_offset(size_type offset) const
+    {
+	check(); MTL_DEBUG_THROW_IF(offset >= this->my_nnz, index_out_of_range("Offset larger than matrix"));
+	return data[offset];
+    }
+
+    value_type& value_from_offset(size_type offset)
+    {
+	check(); MTL_DEBUG_THROW_IF(offset >= this->my_nnz, index_out_of_range("Offset larger than matrix"));
+	return data[offset];
+    }
+
+    /// Swap matrices
+    friend void swap(self& matrix1, self& matrix2)
+    {
+	using std::swap;
+	swap(static_cast<super&>(matrix1), static_cast<super&>(matrix2));
+
+	swap(matrix1.data, matrix2.data);
+	swap(matrix1.starts, matrix2.starts);
+	swap(matrix1.indices, matrix2.indices);
+	swap(matrix1.inserting, matrix2.inserting);
+    }
+
+    /// Remove zero entries
+    void crop()
+    {
+	check();
+	if (data.empty()) return;
+
+	using math::zero;
+	value_type z= zero(data[0]);
+	size_type nzi= 0; // Where to copy next non-zero
+
+	std::vector<size_type>  new_starts(this->dim1() + 1);
+	new_starts[0] = 0;
+
+	for (size_type i = 0; i < this->dim1(); i++) {
+	    for (size_type j= starts[i], end= starts[i+1]; j != end; ++j)
+		if (data[j] != z)
+		    indices[nzi]= indices[j], data[nzi++]= data[j];
+	    new_starts[i+1]= nzi;
+	}
+	this->my_nnz= nzi;
+	data.resize(nzi);
+	indices.resize(nzi);
+	swap(starts, new_starts);
+    }
+
+
+    /// Address of first major index; to be used with care. [advanced]
+    size_type* address_major() { check(); return &starts[0]; }
+    /// Address of first major index; to be used with care. [advanced]
+    const size_type* address_major() const { check(); return &starts[0]; }
+    /// Address of first minor index; to be used with care. [advanced]
+    size_type* address_minor() { check(); return &indices[0]; }
+    /// Address of first minor index; to be used with care. [advanced]
+    const size_type* address_minor() const { check(); return &indices[0]; }
+    /// Address of first data entry; to be used with care. [advanced]
+    value_type* address_data() { check(); return &data[0]; }
+    /// Address of first data entry; to be used with care. [advanced]
+    const value_type* address_data() const { check(); return &data[0]; }
+
+    const std::vector<size_type>& ref_major() const { return starts; } ///< Refer start vector [advanced]
+          std::vector<size_type>& ref_major()       { return starts; } ///< Refer start vector [advanced]
+    const std::vector<size_type>& ref_minor() const { return indices; } ///< Refer index vector [advanced]
+          std::vector<size_type>& ref_minor()       { return indices; } ///< Refer index vector [advanced]
+
+    /// Release unused space in STL vectors
+    void shrink()
+    {
+	shrink_stl_vector(data);
+	shrink_stl_vector(starts);
+	shrink_stl_vector(indices);
+    }
+
+    /// Number of non-zeros in row/column \p r_or_c when matrix is row-/column-major
+    size_type nnz_local(size_type r_or_c) const
+    {
+	MTL_DEBUG_THROW_IF(r_or_c >= this->dim1(), index_out_of_range());
+	return starts[r_or_c+1] - starts[r_or_c];
+    }
+
+    template <typename> friend struct compressed2D_indexer;
+    template <typename, typename, typename> friend struct compressed2D_inserter;
+    template <typename, typename> friend struct compressed_el_cursor;
+    template <typename, typename> friend struct compressed_minor_cursor;
+
+    indexer_type            indexer;
+    std::vector<value_type> data;
+  protected:
+    std::vector<size_type>  starts;
+    std::vector<size_type>  indices;
+    bool                    inserting;
+};
+
+
+
+// ========
+// Inserter
+// ========
+
+/// Additional data structure to insert entries into a compressed2D matrix
+template <typename Elt, typename Parameters, typename Updater = mtl::operations::update_store<Elt> >
+struct compressed2D_inserter
+{
+    typedef compressed2D_inserter             self;
+    typedef compressed2D<Elt, Parameters>     matrix_type;
+    typedef typename matrix_type::size_type   size_type;
+    typedef typename matrix_type::value_type  value_type;
+    typedef std::pair<size_type, size_type>   size_pair;
+    typedef std::map<size_pair, value_type>   map_type;
+    typedef operations::update_proxy<self, size_type>   proxy_type;
+
+  private:
+    // stretch matrix rows or columns to slot size (or leave it if equal or greater)
+    void stretch();
+
+    struct bracket_proxy
+    {
+	bracket_proxy(self& ref, size_type row) : ref(ref), row(row) {}
+
+	proxy_type operator[](size_type col) { return proxy_type(ref, row, col); }
+
+	self&      ref;
+	size_type  row;
+    };
+
+  protected:
+    void finish()
+    {
+	// std::cout << "~compressed2D_inserter: " << matrix.my_nnz << " entries in matrix already (reserved for " <<  elements.size()
+	// 	  << ") and " << spare.size() << " entries in map.\n";
+	vampir_trace<3051> tracer;
+	if (num_rows(matrix) > 0 && num_cols(matrix) > 0) {
+	    final_place();
+	    insert_spare();
+	}
+	matrix.inserting = false;
+    }
+
+  public:
+    /// Construction with matrix reference and optional slot size, see \ref matrix_insertion
+    explicit compressed2D_inserter(matrix_type& matrix, size_type slot_size = 5)
+	: matrix(matrix), elements(matrix.data), starts(matrix.starts), indices(matrix.indices),
+	  slot_size(std::min(slot_size, matrix.dim2())), slot_ends(matrix.dim1()+1)
+    {
+	vampir_trace<3050> tracer;
+	MTL_THROW_IF(matrix.inserting, runtime_error("Two inserters on same matrix"));
+	matrix.inserting = true;
+	if (size(matrix) > 0)
+	    stretch();
+    }
+
+    ~compressed2D_inserter()
+    {
+	// Check if finish wasn't called explicitly before
+	if (matrix.inserting)
+	    finish();
+    }
+
+    /// Proxy to insert into A[row][col]
+    bracket_proxy operator[] (size_type row)
+    {
+	return bracket_proxy(*this, row);
+    }
+
+    /// Proxy to insert into A[row][col]
+    proxy_type operator() (size_type row, size_type col)
+    {
+	return proxy_type(*this, row, col);
+    }
+
+    /// Modify A[row][col] with \p val using \p Modifier
+    template <typename Modifier>
+    void modify(size_type row, size_type col, value_type val);
+
+    /// Modify A[row][col] with \p val using the class' updater
+    void update(size_type row, size_type col, value_type val)
+    {
+	using math::zero;
+	modify<Updater>(row, col, val);
+    }
+
+    /// For debugging only; print entries in all slots; ignores entries in spare map [advanced]
+    void print() const
+    {
+	for (size_type j= 0; j < matrix.dim1(); j++)
+	    print(j);
+    }
+
+    /// For debugging only; print entries in slot; ignores entries in spare map [advanced]
+    void print(size_type i) const
+    {
+	std::cout << "in slot " << i << ": ";
+	for (size_type j= starts[i]; j < slot_ends[i]; ++j)
+	    std::cout << "[" << indices[j] << ": " << elements[j] << "]";
+	std::cout << "\n";
+    }
+
+    /// Empties slot i (row or column according to orientation); for experts only [advanced]
+    /** Does not work if entries are in spare map !!!! **/
+    void make_empty(size_type i)
+    {	slot_ends[i]= starts[i];    }
+
+    /// Value in A[r][c]
+    value_type value(size_type r, size_type c) const
+    {
+	size_pair                      mm= matrix.indexer.major_minor_c(matrix, r, c);
+	utilities::maybe<size_type>    offset= matrix_offset(mm);
+	return offset ? matrix.data[offset.value()] : value_type(0);
+    }
+
+    /// Insert \p elements into %matrix with pre-sorting
+    /** It should be faster for sufficiently large element matrices but our benchmarks showed the opposite. **/
+    template <typename Matrix, typename Rows, typename Cols>
+    self& sorted_block_insertion(const element_matrix_t<Matrix, Rows, Cols>& elements);
+
+    /// Insert \p elements into %matrix
+    template <typename Matrix, typename Rows, typename Cols>
+    self& operator<< (const element_matrix_t<Matrix, Rows, Cols>& elements)
+    {
+	using mtl::size;
+#if 0 // shouldn't be slower now
+	if (size(elements.cols) > sorted_block_insertion_limit
+	    && boost::is_same<typename Parameters::orientation, row_major>::value)
+	    return sorted_block_insertion(elements);
+#endif
+
+	for (unsigned ri= 0; ri < size(elements.rows); ri++)
+	    for (unsigned ci= 0; ci < size(elements.cols); ci++)
+		update (elements.rows[ri], elements.cols[ci], elements.matrix(ri, ci));
+	return *this;
+    }
+
+    /// Insert \p elements into %matrix
+    template <typename Matrix, typename Rows, typename Cols>
+    self& operator<< (const element_array_t<Matrix, Rows, Cols>& elements)
+    {
+	using mtl::size;
+	for (unsigned ri= 0; ri < size(elements.rows); ri++)
+	    for (unsigned ci= 0; ci < size(elements.cols); ci++)
+		update (elements.rows[ri], elements.cols[ci], elements.array[ri][ci]);
+	return *this;
+    }
+
+    // not so nice functions needed for direct access, e.g. in factorizations
+    std::vector<size_type> const& ref_major() const { return starts; } ///< Refer start vector [advanced]
+    std::vector<size_type> const& ref_minor() const { return indices; } ///< Refer index vector [advanced]
+    std::vector<size_type> const& ref_slot_ends() const { return slot_ends; } ///< Refer slot-end vector [advanced]
+    std::vector<value_type> const& ref_elements() const { return elements; } ///< Refer element vector [advanced]
+
+  private:
+    utilities::maybe<typename self::size_type> matrix_offset(size_pair) const;
+    void final_place();
+    void insert_spare();
+
+    // Show ith row/column
+    void display(size_type i)
+    {
+	std::cout << "slot " << i << " is [";
+	for (size_type j= starts[i]; j < slot_ends[i]; ++j)
+	    std::cout << '(' << indices[j] << ": " <<  elements[j] << ")"
+		      << (j < slot_ends[i] - 1 ? ", " : "");
+	std::cout << "]\n";
+    }
+
+    // eliminate double entries in a slot by using the updater
+    // indices are supposed to be sorted
+    void reduce_slot(size_type i)
+    {
+	size_type tgt= starts[i], src= tgt + 1, end= slot_ends[i];
+	for (; src < end;) {
+	    while (src < end && indices[tgt] == indices[src])
+		Updater()(elements[tgt], elements[src++]);
+	    ++tgt;
+	    if (tgt == src)
+		++src;
+	    else if (src < end && indices[tgt] != indices[src]) {
+		indices[tgt]= indices[src];
+		elements[tgt]= elements[src++];
+		if (src >= end) ++tgt; // correction to not remain on the last entry
+	    }
+	}
+	slot_ends[i]= tgt;
+    }
+
+  protected:
+    compressed2D<Elt, Parameters>&      matrix;
+    std::vector<value_type>&            elements;
+    std::vector<size_type>&             starts;
+    std::vector<size_type>&             indices;
+    size_type                           slot_size;
+    std::vector<size_type>              slot_ends;
+    map_type                            spare;
+};
+
+template <typename Elt, typename Parameters, typename Updater>
+void compressed2D_inserter<Elt, Parameters, Updater>::stretch()
+{
+    using std::copy;
+    using std::copy_backward;
+    using std::swap;
+
+    vampir_trace<3052> tracer;
+    // Stretching is much simpler for empty matrices
+    if (elements.empty()) {
+	for (size_type i= 0, s= 0; i <= matrix.dim1(); i++, s+= slot_size)
+	    slot_ends[i]= starts[i]= s;
+	size_type new_total= (slot_ends[matrix.dim1()]= starts[matrix.dim1()]) + slot_size;
+	elements.reserve(new_total); indices.reserve(new_total);
+	elements.resize(new_total); indices.resize(new_total);
+	return;
+    }
+
+
+
+    // If there are enough existing entries then skip the stretching (expensive)
+    if (elements.size() + matrix.dim1()/2 > std::size_t(slot_size * matrix.dim1())) {
+	// Use start of next row/col as slot_ends
+	copy(starts.begin() + 1, starts.end(), slot_ends.begin());
+	slot_ends[matrix.dim1()]= starts[matrix.dim1()];
+	return;
+    }
+
+    std::vector<size_type>  new_starts(matrix.dim1() + 1);
+    new_starts[0] = 0;
+    for (size_type i = 0; i < matrix.dim1(); i++) {
+	size_type entries = starts[i+1] - starts[i];
+	slot_ends[i] = new_starts[i] + entries;
+	new_starts[i+1] = new_starts[i] + std::max(entries, slot_size);
+    }
+    // Add an additional slot for temporaries
+    size_type new_total= (slot_ends[matrix.dim1()]= new_starts[matrix.dim1()]) + slot_size;
+    elements.resize(new_total);
+    indices.resize(new_total);
+    // for (int i= 0; i < matrix.dim1()+1; i++) std::cout << "Slot " << i << " is [" << new_starts[i] << ", " << slot_ends[i] << ")\n";
+
+    // copy normally if not overlapping and backward if overlapping
+    // i goes down to 1 (not to 0) because i >= 0 never stops for unsigned ;-)
+	// &v[i] is replaced by &v[0]+i to enable past-end addresses for STL copy
+    for (size_type i = matrix.dim1(); i > 0; i--)
+	if (starts[i] <= new_starts[i-1]) {
+	    // std::cout << "Kopiere vorwaerts von " << starts[i-1] << " bis " << starts[i] << " nach " << new_starts[i-1] << "\n";
+	    copy(&elements[0] + starts[i-1], &elements[0] + starts[i], &elements[0] + new_starts[i-1]);
+	    copy(&indices[0] + starts[i-1], &indices[0] + starts[i], &indices[0] + new_starts[i-1]);
+	} else {
+	    // std::cout << "Kopiere rueckwaerts von " << starts[i-1] << " bis " << starts[i] << " nach " << slot_ends[i-1] << "\n";
+	    copy_backward(&elements[0] + starts[i-1], &elements[0] + starts[i], &elements[0] + slot_ends[i-1]);
+	    copy_backward(&indices[0] + starts[i-1], &indices[0] + starts[i], &indices[0] + slot_ends[i-1]);
+	}
+    swap(starts, new_starts);
+}
+
+template <typename Elt, typename Parameters, typename Updater>
+inline utilities::maybe<typename compressed2D_inserter<Elt, Parameters, Updater>::size_type>
+compressed2D_inserter<Elt, Parameters, Updater>::matrix_offset(size_pair mm) const
+{
+    size_type major, minor;
+    boost::tie(major, minor) = mm;
+
+    if (indices.empty())
+	return utilities::maybe<size_type> (0, false);
+
+    // &v[i] isn't liked by all libs -> &v[0]+i circumvents complaints
+    const size_type *first = &indices[0] + starts[major],
+  	            *last =  &indices[0] + slot_ends[major];
+    if (first == last)
+	return utilities::maybe<size_type> (first - &indices[0], false);
+
+    const size_type *index= first;
+    if (last - index < 10)
+	while (index != last && *index < minor) ++index;
+    else
+	index = std::lower_bound(first, last, minor);
+    return utilities::maybe<size_type>(index - &indices[0], index != last && *index == minor);
+}
+
+
+template <typename Elt, typename Parameters, typename Updater>
+template <typename Modifier>
+inline void compressed2D_inserter<Elt, Parameters, Updater>::modify(size_type row, size_type col, value_type val)
+{
+    using std::copy_backward;
+    MTL_DEBUG_THROW_IF(is_negative(row) || row >= num_rows(matrix) || is_negative(col) || col >= num_cols(matrix), index_out_of_range());
+
+    Modifier                          modifier;
+    compressed2D_indexer<size_type>   indexer;
+    size_pair                         mm = indexer.major_minor_c(matrix, row, col);
+    size_type                         major, minor;
+    boost::tie(major, minor) = mm;
+
+    utilities::maybe<size_type>       pos = matrix_offset(mm);
+    // Check if already in matrix and update it
+    if (pos)
+	modifier (elements[pos.value()], val);
+    else {
+	size_type& my_end = slot_ends[major];
+	// Check if place in matrix to insert there
+	if (my_end != starts[major+1]) {
+	    if (pos.value() != my_end) {
+		copy_backward(&elements[0] + pos.value(), &elements[0] + my_end, &elements[0] + (my_end+1));
+		copy_backward(&indices[0] + pos.value(), &indices[0] + my_end, &indices[0] + (my_end+1));
+	    }
+	    elements[pos.value()] = modifier.init(val); indices[pos.value()] = minor;
+	    my_end++;
+	    matrix.my_nnz++;      // new entry
+	} else {
+	    typename map_type::iterator it = spare.find(mm);
+	    // If not in map insert it, otherwise update the value
+	    if (it == spare.end()) {
+		spare.insert(std::make_pair(mm, modifier.init(val)));
+		matrix.my_nnz++;      // new entry
+	    } else
+		modifier(it->second, val);
+	}
+    }
+    // std::cout << "inserter update: " << matrix.my_nnz << " non-zero elements, new value is " << elements[pos] << "\n";
+}
+
+namespace detail {
+
+    struct cmp_first
+    {
+	template <typename F, typename S>
+	bool operator()(const std::pair<F, S>& x, const std::pair<F, S>& y) const
+	{
+	    return x.first < y.first;
+	}
+    };
+}
+
+template <typename Elt, typename Parameters, typename Updater>
+template <typename Matrix, typename Rows, typename Cols>
+compressed2D_inserter<Elt, Parameters, Updater>&
+compressed2D_inserter<Elt, Parameters, Updater>::sorted_block_insertion(const element_matrix_t<Matrix, Rows, Cols>& iblock)
+{
+    using std::copy; using std::copy_backward; using std::min;
+    using mtl::size; using mtl::num_rows; using mtl::num_cols;
+    using namespace mtl::utility;
+    typedef zip_it<size_type, value_type> it_type;
+
+    size_type m= size(iblock.rows), n= size(iblock.cols);
+    MTL_THROW_IF(m != num_rows(iblock.matrix) || n != num_cols(iblock.matrix), incompatible_size());
+    MTL_THROW_IF(&iblock.matrix[0][1] - &iblock.matrix[0][0] != 1, logic_error("Rows must be consecutive"));
+
+    size_type rmax= matrix.dim1(), &index_max0= indices[starts[rmax]];
+    value_type& value_max0= elements[starts[rmax]];
+    for (size_type i= 0; i < m; i++) {
+	size_type r= iblock.rows[i], &index_0= indices[starts[r]];
+	value_type& value_0= elements[starts[r]];
+	if (slot_ends[r] == starts[r]) {
+	    size_type to_copy= min(starts[r+1] - starts[r], n);
+	    copy(&iblock.cols[0], &iblock.cols[0]+to_copy, &index_0);
+	    copy(&iblock.matrix[i][0], &iblock.matrix[i][0]+to_copy, &value_0);
+	    slot_ends[r]= starts[r] + to_copy;
+	    std::sort(it_type(&index_0, &value_0, 0), it_type(&index_0, &value_0, to_copy), less_0());
+	    reduce_slot(r);
+	    // display(r);
+	    matrix.my_nnz+= slot_ends[r] - starts[r];
+	    for (size_type j= to_copy; j < n; ++j)
+		update(r, iblock.cols[j], iblock.matrix[i][j]);
+	} else {
+	    size_type to_copy= min(slot_size, n);
+	    copy(&iblock.cols[0], &iblock.cols[0]+to_copy, &index_max0);
+	    copy(&iblock.matrix[i][0], &iblock.matrix[i][0]+to_copy, &value_max0);
+	    slot_ends[rmax]= starts[rmax] + to_copy;
+	    std::sort(it_type(&index_max0, &value_max0, 0), it_type(&index_max0, &value_max0, to_copy), less_0());
+	    size_type tgt= starts[r], tend= slot_ends[r], later= starts[rmax], src= later, end= slot_ends[rmax];
+	    for (; src < end; ) {
+		// search for next equal entry in slot r
+		while (tgt < tend && indices[src] > indices[tgt]) ++tgt;
+		// reduce equal entries (they are consecutive)
+		if (tgt < tend)
+		    while (src < end && indices[src] == indices[tgt])
+			Updater()(elements[tgt], elements[src++]);
+		else { // all new entries are larger than the largest existing
+		    for (; src < end && tgt < starts[r+1]; ++src) { // copy at the end of slot
+			indices[tgt]= indices[src];
+			elements[tgt]= elements[src];
+			slot_ends[r]= ++tgt;
+			++matrix.my_nnz;
+		    }
+		    for (; src < end; ++src) // remainder goes into spare
+			update(r, indices[src], elements[src]);
+		    later= starts[rmax]; // we're done, avoid postponed loop
+		}
+		// collect indices not found in r to deal with later
+		while (src < end && indices[src] < indices[tgt]) {
+		    if (later != src) {
+			indices[later]= indices[src];
+			elements[later]= elements[src];
+		    }
+		    ++src; ++later;
+		}
+	    }
+	    for (size_type j= starts[rmax]; j < later; ++j)
+		update(r, indices[j], elements[j]);
+	}
+    }
+    return *this;
+}
+
+
+template <typename Elt, typename Parameters, typename Updater>
+void compressed2D_inserter<Elt, Parameters, Updater>::final_place()
+{
+    using std::swap;
+    vampir_trace<3053> tracer;
+
+    size_type          dim1 = matrix.dim1();
+    std::vector<size_type>  new_starts(dim1 + 1);
+    new_starts[0] = 0;
+
+    if (spare.empty()) {
+	// Check if everything is already in place
+	if (slot_ends[dim1-1] == matrix.my_nnz) {
+	    starts[dim1]= slot_ends[dim1-1];
+	    if (std::size_t(matrix.my_nnz) < elements.size())
+		elements.resize(matrix.my_nnz),
+		    indices.resize(matrix.my_nnz);
+	    return;
+	}
+	size_type pos= 0;
+	for (size_type i= 0; i < dim1; ++i) {
+	    new_starts[i]= pos;
+	    for (size_type j= starts[i], je= slot_ends[i]; j < je; ++j, ++pos) {
+		elements[pos]= elements[j];
+		indices[pos]= indices[j];
+	    }
+	}
+	new_starts[dim1]= pos;
+	swap(new_starts, starts);
+	if (std::size_t(matrix.my_nnz) < elements.size())
+	    elements.resize(matrix.my_nnz),
+		indices.resize(matrix.my_nnz);
+	return;
+    }
+
+    typename map_type::iterator it = spare.begin();
+    for (size_type i = 0; i < dim1; i++) {
+	size_type entries = slot_ends[i] - starts[i];
+	while (it != spare.end() && it->first.first == i)
+	    entries++, it++;
+	new_starts[i+1] = new_starts[i] + entries;
+    }
+
+    size_type new_total = new_starts[dim1], old_total = starts[dim1];
+    if (new_total > old_total) {
+	elements.resize(new_total);
+	indices.resize(new_total); }
+
+    operations::shift_blocks(dim1, starts, new_starts, slot_ends, elements);
+    operations::shift_blocks(dim1, starts, new_starts, slot_ends, indices);
+
+    if (std::size_t(new_total) < elements.size()) {
+	elements.resize(new_total);
+	indices.resize(new_total); }
+
+    for (size_type i = 0; i < dim1; i++)
+	slot_ends[i] = new_starts[i] + slot_ends[i] - starts[i];
+
+    swap(starts, new_starts);
+}
+
+template <typename Elt, typename Parameters, typename Updater>
+void compressed2D_inserter<Elt, Parameters, Updater>::insert_spare()
+{
+    vampir_trace<3054> tracer;
+    using std::copy_backward;
+
+    for (typename map_type::iterator it = spare.begin(); it != spare.end(); ++it) {
+	size_pair              mm = it->first;
+	size_type              major = mm.first, minor = mm.second;
+	utilities::maybe<size_type>       pos = matrix_offset(mm);
+	size_type&             my_end = slot_ends[major];
+
+	// &v[i] see above
+	copy_backward(&elements[0] + pos.value(), &elements[0] + my_end, &elements[0] + (my_end+1));
+	copy_backward(&indices[0] + pos.value(), &indices[0] + my_end, &indices[0] + (my_end+1));
+	elements[pos.value()] = it->second; indices[pos.value()] = minor;
+	my_end++;
+    }
+}
+
+// ================
+// Free functions
+// ================
+
+template <typename Value, typename Parameters>
+typename compressed2D<Value, Parameters>::size_type
+inline num_rows(const compressed2D<Value, Parameters>& matrix)
+{
+    return matrix.num_rows();
+}
+
+template <typename Value, typename Parameters>
+typename compressed2D<Value, Parameters>::size_type
+inline num_cols(const compressed2D<Value, Parameters>& matrix)
+{
+    return matrix.num_cols();
+}
+
+template <typename Value, typename Parameters>
+// typename compressed2D<Value, Parameters>::size_type risks overflow
+std::size_t
+inline size(const compressed2D<Value, Parameters>& matrix)
+{
+    return std::size_t(matrix.num_cols()) * std::size_t(matrix.num_rows());
+}
+
+
+}} // namespace mtl::matrix
+
+namespace mtl {
+	using matrix::compressed2D;
+}
+
+// ================
+// Range generators
+// ================
+
+namespace mtl { namespace traits {
+
+    // VC 8.0 finds ambiguity with mtl::tag::morton_dense (I wonder why, especially here)
+    using mtl::matrix::compressed2D;
+    using mtl::matrix::compressed_el_cursor;
+    using mtl::matrix::compressed_minor_cursor;
+
+    // ===========
+    // For cursors
+    // ===========
+
+    template <class Elt, class Parameters>
+    struct range_generator<glas::tag::nz, compressed2D<Elt, Parameters> >
+      : detail::all_offsets_range_generator<compressed2D<Elt, Parameters>,
+					    compressed_el_cursor<Elt, Parameters>,
+					    complexity_classes::linear_cached>
+    {};
+
+    // Cursor over all rows
+    // Supported if row major matrix
+    template <typename Elt, typename Parameters>
+    struct range_generator<glas::tag::row, compressed2D<Elt, Parameters> >
+      : boost::mpl::if_<
+	    boost::is_same<typename Parameters::orientation, row_major>
+ 	  , detail::all_rows_range_generator<compressed2D<Elt, Parameters>, complexity_classes::linear_cached>
+ 	  , range_generator<tag::unsupported, compressed2D<Elt, Parameters> >
+        >::type {};
+
+
+    template <class Elt, class Parameters>
+    struct range_generator<glas::tag::nz,
+			   detail::sub_matrix_cursor<compressed2D<Elt, Parameters>, glas::tag::row, 2> >
+    {
+	typedef typename Collection<compressed2D<Elt, Parameters> >::size_type  size_type;
+	typedef detail::sub_matrix_cursor<compressed2D<Elt, Parameters>, glas::tag::row, 2> cursor_type;
+	typedef complexity_classes::linear_cached         complexity;
+	typedef compressed_minor_cursor<Elt, Parameters>  type;
+	static int const                                  level = 1;
+
+	type begin(cursor_type const& cursor) const
+	{
+	    return type(cursor.ref, cursor.key, cursor.ref.begin_col());
+	}
+	type end(cursor_type const& cursor) const
+	{
+	    return type(cursor.ref, cursor.key, cursor.ref.end_col());
+	}
+	type lower_bound(cursor_type const& cursor, size_type position) const
+	{
+	    return type(cursor.ref, cursor.key, std::min(position, cursor.ref.end_col()));
+	}
+    };
+
+    // Cursor over all columns
+    // Supported if column major matrix
+    template <typename Elt, typename Parameters>
+    struct range_generator<glas::tag::col, compressed2D<Elt, Parameters> >
+      : boost::mpl::if_<
+	    boost::is_same<typename Parameters::orientation, col_major>
+ 	  , detail::all_cols_range_generator<compressed2D<Elt, Parameters>, complexity_classes::linear_cached>
+ 	  , range_generator<tag::unsupported, compressed2D<Elt, Parameters> >
+        >::type {};
+
+
+    template <class Elt, class Parameters>
+    struct range_generator<glas::tag::nz,
+			   detail::sub_matrix_cursor<compressed2D<Elt, Parameters>, glas::tag::col, 2> >
+    {
+	typedef typename Collection<compressed2D<Elt, Parameters> >::size_type  size_type;
+	typedef detail::sub_matrix_cursor<compressed2D<Elt, Parameters>, glas::tag::col, 2> cursor_type;
+	typedef complexity_classes::linear_cached         complexity;
+	typedef compressed_minor_cursor<Elt, Parameters>  type;
+	static int const                                  level = 1;
+
+	type begin(cursor_type const& cursor) const
+	{
+	    return type(cursor.ref, cursor.ref.begin_row(), cursor.key);
+	}
+	type end(cursor_type const& cursor) const
+	{
+	    return type(cursor.ref, cursor.ref.end_row(), cursor.key);
+	}
+	type lower_bound(cursor_type const& cursor, size_type position) const
+	{
+	    return type(cursor.ref, std::min(position, cursor.ref.end_row()), cursor.key);
+	}
+    };
+
+    // Cursor over all rows or columns, depending which one is major
+    template <typename Elt, typename Parameters>
+    struct range_generator<glas::tag::major, compressed2D<Elt, Parameters> >
+      : boost::mpl::if_<
+	    boost::is_same<typename Parameters::orientation, row_major>
+	  , range_generator<glas::tag::row, compressed2D<Elt, Parameters> >
+	  , range_generator<glas::tag::col, compressed2D<Elt, Parameters> >
+        >::type {};
+
+
+// =============
+// For iterators
+// =============
+
+
+    template <class Elt, class Parameters>
+    struct range_generator<tag::const_iter::nz,
+			   detail::sub_matrix_cursor<compressed2D<Elt, Parameters>, glas::tag::row, 2> >
+    {
+	typedef compressed2D<Elt, Parameters>                                         matrix_type;
+	typedef typename matrix_type::size_type                                       size_type;
+	typedef typename matrix_type::value_type                                      value_type;
+	typedef detail::sub_matrix_cursor<matrix_type, glas::tag::row, 2>             cursor;
+
+	typedef complexity_classes::linear_cached                                     complexity;
+	static int const                                                              level = 1;
+	typedef const value_type*                                                     type;
+
+	type begin(cursor const& c)
+	{
+	    const matrix_type& matrix= c.ref;
+	    size_type offset= matrix.indexer(matrix, c.key, matrix.begin_col());
+	    return &matrix.data[0] + offset;
+	}
+
+	// returned pointer can pass the end and must only be used for comparison
+	type end(cursor const& c)
+	{
+	    const matrix_type& matrix= c.ref;
+	    size_type offset= matrix.indexer(matrix, c.key, matrix.end_col());
+	    return &matrix.data[0] + offset;
+	}
+    };
+
+
+    template <class Elt, class Parameters>
+    struct range_generator<tag::const_iter::nz,
+			   detail::sub_matrix_cursor<compressed2D<Elt, Parameters>, glas::tag::col, 2> >
+    {
+	typedef compressed2D<Elt, Parameters>                                         matrix_type;
+	typedef typename matrix_type::size_type                                       size_type;
+	typedef typename matrix_type::value_type                                      value_type;
+	typedef detail::sub_matrix_cursor<matrix_type, glas::tag::col, 2>             cursor;
+
+	typedef complexity_classes::linear_cached                                     complexity;
+	static int const                                                              level = 1;
+	typedef const value_type*                                                     type;
+
+	type begin(cursor const& c)
+	{
+	    const matrix_type& matrix= c.ref;
+	    size_type offset= matrix.indexer(matrix, matrix.begin_row(), c.key);
+	    return &matrix.data[0] + offset;
+	}
+
+	// returned pointer can pass the end and must only be used for comparison
+	type end(cursor const& c)
+	{
+	    const matrix_type& matrix= c.ref;
+	    size_type offset= matrix.indexer(matrix, matrix.end_row(), c.key);
+	    return &matrix.data[0] + offset;
+	}
+    };
+
+
+}} // namespace mtl::traits
+
+
+
+#endif // MTL_COMPRESSED2D_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/coordinate2D.hpp

+
+#ifndef MTL_COORDINATE2D_INCLUDE
+#define MTL_COORDINATE2D_INCLUDE
+
+#include <vector>
+#include <cassert>
+#include <boost/mpl/bool.hpp>
+
+
+#include <boost/numeric/mtl/operation/is_negative.hpp>
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/operation/sort.hpp>
+#include <boost/numeric/mtl/operation/iota.hpp>
+#include <boost/numeric/mtl/matrix/parameter.hpp>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+#include <boost/numeric/mtl/utility/is_row_major.hpp>
+#include <boost/numeric/mtl/utility/static_assert.hpp>
+
+namespace mtl {  namespace matrix {
+
+
+/// Sparse matrix structure in coordinate format
+template <typename Value, typename Parameters = matrix::parameters<> >
+class coordinate2D
+  : public base_matrix<Value, Parameters>,
+    public const_crtp_base_matrix< coordinate2D<Value, Parameters>, Value, typename Parameters::size_type >,
+    public crtp_matrix_assign< coordinate2D<Value, Parameters>, Value, typename Parameters::size_type >,
+    public mat_expr< coordinate2D<Value, Parameters> >
+{
+  public:
+
+    typedef Value                                                       value_type ;
+    typedef Value&                                                      reference ;
+    typedef Value const&                                                const_reference ;
+    typedef typename Parameters::size_type				size_type ;
+    typedef typename Parameters::dimensions                             dim_type;
+    typedef typename Parameters::orientation                            orientation;
+
+    typedef base_matrix<Value, Parameters>                              super;
+
+    typedef std::vector< size_type >                                    row_index_array_type ;
+    typedef std::vector< size_type >                                    column_index_array_type ;
+    typedef std::vector< value_type >                                   value_array_type ;
+
+    /// Common constructor
+    explicit coordinate2D(size_type nrows, size_type ncols, size_type expected= 0)
+      : super(dim_type(nrows, ncols))
+    {
+	if (expected > 0) {
+	    rows.reserve(expected);
+	    cols.reserve(expected);
+	    values.reserve(expected);
+	}
+    }
+
+    size_type nnz() const { return rows.size(); } ///< Number of non-zeros
+
+    value_array_type const& value_array() const { return values; } ///< Array of values (const)
+    value_array_type& value_array() { return values; }             ///< Array of values (mutable)
+
+    row_index_array_type const& row_index_array() const { return rows; }       ///< Array of rows  (const)
+    column_index_array_type const& column_index_array() const {	return cols; } ///< Array of columns (const)
+
+    row_index_array_type& row_index_array() { return rows; }       ///< Array of rows   (mutable)
+    column_index_array_type& column_index_array() { return cols; } ///< Array of columns  (mutable)
+
+    /// Insert an entry at the end of the row-,col- and value-array
+    void push_back(size_type r, size_type c, const_reference v)
+    {
+	rows.push_back(r); cols.push_back(c); values.push_back(v); my_is_sorted= false;
+    }
+
+    /// Insert an entry at the end of the row-,col- and value-array, like push_back
+    void insert(size_type r, size_type c, const_reference v) {	push_back(r, c, v); }
+
+    /// Whether the entries are sorted
+    bool is_sorted() const { return my_is_sorted; }
+
+    /// sorting standard by rows
+    void sort()
+    {
+	if (nnz() > 0)
+	    sort(mtl::traits::is_row_major<Parameters>());
+	my_is_sorted= true;
+    }
+
+  private:
+    // sorting by rows
+    void sort(boost::mpl::true_)
+    {
+	mtl::vector::sort_xy(rows, cols, values);
+    }
+
+    // sorting by columns
+    void sort(boost::mpl::false_)
+    {
+	mtl::vector::sort_xy(cols, rows, values);
+    }
+
+    template <typename OStream, typename Vector>
+    void print_stl_vector(OStream& os, const Vector& v) const
+    {
+	os << "[";
+	for (unsigned i= 0; i < v.size(); i++)
+	    os << v[i] << (i+1 < v.size() ? "," : "");
+	os << "]\n";
+    }
+
+  public:
+    template <typename OStream>
+    void print_internal(OStream& os) const
+    {
+       	os << "rows   = "; print_stl_vector(os, rows);
+	os << "cols   = "; print_stl_vector(os, cols);
+	os << "values = "; print_stl_vector(os, values);
+    }
+
+    void print_internal() const { print(std::cout); }
+
+    ///operator * for  vector= coordinaten-matrix * vector
+    template <typename Vector >
+    Vector operator*(const Vector& x)
+    {
+
+	Vector res(this->num_rows());
+	res= 0;
+	for (size_type i= 0; i < nnz(); i++)
+	    res[rows[i]]+=  values[i] * x[cols[i]];
+	return res;
+    }
+
+    value_type operator() (const size_type r, const size_type c) const
+    {
+	MTL_DEBUG_THROW_IF(is_negative(r) || r >= this->num_rows()
+			   || is_negative(c) || c >= this->num_cols(), index_out_of_range());
+
+#if 0
+	if (my_is_sorted)
+	    return find(r, c, mtl::traits::is_row_major<Parameters>());
+#endif
+
+	for (size_type i= 0; i < nnz(); i++)
+	    if (rows[i] == r && cols[i] == c)
+		return values[i];
+	return value_type(0);
+    }
+
+    template <typename Updater>
+    void compress(Updater up)
+    {
+	if (!my_is_sorted)
+	    sort();
+
+	size_type i= 0, j= 1, end= rows.size();
+	for (; j < end; ++j)
+	    if (rows[i] == rows[j] && cols[i] == cols[j]) {
+		up(values[i], values[j]);
+	    } else {
+		i++;
+		if (i != j) {
+		    rows[i]= rows[j];
+		    cols[i]= cols[j];
+		    values[i]= values[j];
+		}
+	    }
+	if (end > 0) i++;
+	rows.resize(i);
+	cols.resize(i);
+	values.resize(i);
+   }
+
+    template <typename Matrix, typename Updater> friend struct coordinate2D_inserter;
+
+  private:
+    row_index_array_type      rows;
+    column_index_array_type   cols;
+    value_array_type          values;
+    bool                      my_is_sorted;
+};
+
+// ================
+// Free functions
+// ================
+
+
+/// Number of rows
+template <typename Value, typename Parameters>
+typename coordinate2D<Value, Parameters>::size_type
+inline num_rows(const coordinate2D<Value, Parameters>& matrix)
+{
+    return matrix.num_rows();
+}
+
+/// Number of columns
+template <typename Value, typename Parameters>
+typename coordinate2D<Value, Parameters>::size_type
+inline num_cols(const coordinate2D<Value, Parameters>& matrix)
+{
+    return matrix.num_cols();
+}
+
+/// Size of the matrix, i.e. the number of row times columns
+template <typename Value, typename Parameters>
+typename coordinate2D<Value, Parameters>::size_type
+inline size(const coordinate2D<Value, Parameters>& matrix)
+{
+    return matrix.num_cols() * matrix.num_rows();
+}
+
+/// Number of NoZeros of the matrix
+template <typename Value, typename Parameters>
+typename coordinate2D<Value, Parameters>::size_type
+inline nnz(const coordinate2D<Value, Parameters>& matrix)
+{
+    return matrix.nnz();
+}
+
+
+template <typename Matrix,
+	  typename Updater = mtl::operations::update_store<typename Matrix::value_type> >
+struct coordinate2D_inserter
+{
+    typedef coordinate2D_inserter                       self;
+    typedef Matrix                                      matrix_type;
+    typedef typename matrix_type::size_type             size_type;
+    typedef typename matrix_type::value_type            value_type;
+    typedef operations::update_proxy<self, size_type>   proxy_type;
+
+    // We only support storing so far !!!
+    // STATIC_ASSERT((boost::is_same<Updater, mtl::operations::update_store<value_type> >::value), "We only support storing so far");
+
+    coordinate2D_inserter(matrix_type& matrix, size_type slot_size= 1)
+      : matrix(matrix)
+    {
+	std::size_t ns= slot_size * matrix.dim1();
+	if (ns > matrix.nnz()) {
+	    matrix.rows.reserve(ns);
+	    matrix.cols.reserve(ns);
+	    matrix.values.reserve(ns);
+	}
+    }
+
+    ~coordinate2D_inserter() { matrix.compress(Updater()); }
+
+ private:
+
+    struct update_proxy
+    {
+	// self is type of inserter not update_proxy !!!
+	update_proxy(self& ref, size_type row, size_type col) : ref(ref), row(row), col(col) {}
+
+	template <typename Value>
+	update_proxy& operator<< (Value const& val)
+	{
+	    ref.matrix.push_back(row, col, val);
+	    return *this;
+	}
+	self& ref;
+	size_type row, col;
+    };
+
+    proxy_type operator() (size_type row, size_type col)
+    {
+	return proxy_type(*this, row, col);
+    }
+
+
+    struct bracket_proxy
+    {
+	bracket_proxy(self& ref, size_type row) : ref(ref), row(row) {}
+
+	proxy_type operator[](size_type col)
+	{
+	    return proxy_type(ref, row, col);
+	}
+
+	self&      ref;
+	size_type  row;
+    };
+
+  public:
+
+    bracket_proxy operator[] (size_type row)
+    {
+	return bracket_proxy(*this, row);
+    }
+
+    template <typename Value>
+    void update(size_type row, size_type col, Value val)
+    {
+	matrix.push_back(row, col, val);
+    }
+
+    template <typename Modifier, typename Value>
+    void modify(size_type row, size_type col, Value val)
+    {
+	matrix.push_back(row, col, val);
+    }
+
+    template <typename EMatrix, typename Rows, typename Cols>
+    self& operator<< (const matrix::element_matrix_t<EMatrix, Rows, Cols>& elements)
+    {
+	using mtl::size;
+	for (unsigned ri= 0; ri < size(elements.rows); ri++)
+	    for (unsigned ci= 0; ci < size(elements.cols); ci++)
+		update (elements.rows[ri], elements.cols[ci], elements.matrix(ri, ci));
+	return *this;
+    }
+
+    template <typename EMatrix, typename Rows, typename Cols>
+    self& operator<< (const matrix::element_array_t<EMatrix, Rows, Cols>& elements)
+    {
+	using mtl::size;
+	for (unsigned ri= 0; ri < size(elements.rows); ri++)
+	    for (unsigned ci= 0; ci < size(elements.cols); ci++)
+		update (elements.rows[ri], elements.cols[ci], elements.array[ri][ci]);
+	return *this;
+    }
+
+  protected:
+    matrix_type&         matrix;
+};
+
+struct coordinate_key
+{
+    typedef std::size_t                               size_t;
+
+    explicit coordinate_key(size_t offset) : offset(offset) {}
+
+    bool operator== (coordinate_key const& other) const { return offset == other.offset; }
+    bool operator!= (coordinate_key const& other) const { return offset != other.offset; }
+
+    size_t offset;
+};
+
+
+// Cursor over every element
+template <typename Value, typename Parameters>
+struct coordinate_minor_cursor
+ : public coordinate_key
+{
+    typedef coordinate_minor_cursor<Value, Parameters>   self;
+    typedef typename Parameters::size_type               size_type;
+    typedef const coordinate2D<Value, Parameters>&       matrix_ref_type;
+    static const int                                     level= 2;
+
+    coordinate_minor_cursor(matrix_ref_type ref, size_type offset)
+      : coordinate_key(offset), ref(ref)  {}
+
+    bool operator!=(const self& that) const
+    {
+	assert(&ref == &that.ref);
+	return this->offset != that.offset;
+    }
+
+    self& operator++() { this->offset++; return *this; }
+    coordinate_key operator*() const { return *this; }
+
+    matrix_ref_type ref;
+};
+
+
+template <typename Value, typename Parameters>
+struct coordinate_major_cursor
+{
+    typedef coordinate_major_cursor<Value, Parameters>   self;
+    typedef typename Parameters::size_type               size_type;
+    typedef const coordinate2D<Value, Parameters>&       matrix_ref_type;
+    typedef coordinate_minor_cursor<Value, Parameters>   inner_cursor;
+    static const int                                     level= 2;
+
+    void find_next_offset(boost::mpl::true_)
+    {
+	size_type i= offset;
+	for ( ; i < nnz(ref) && ref.row_index_array()[i] <= major; i++) ;
+	next_offset= i;
+    }
+
+    void find_next_offset(boost::mpl::false_)
+    {
+	size_type i= offset;
+	for ( ; i < nnz(ref) && ref.col_index_array()[i] <= major; i++) ;
+	next_offset= i;
+    }
+
+    void find_next_offset() { find_next_offset(mtl::traits::is_row_major<Parameters>()); }
+
+    coordinate_major_cursor(matrix_ref_type ref, size_type major, size_type offset)
+      : ref(ref), major(major), offset(offset)
+    {
+	find_next_offset();
+    }
+
+    bool operator!=(const self& that) const
+    {
+	assert(&ref == &that.ref);
+	return this->offset != that.offset;
+    }
+
+    self& operator++()
+    {
+	offset= next_offset;
+	major++;
+	find_next_offset();
+	return *this;
+    }
+
+    matrix_ref_type ref;
+    size_type       major, offset, next_offset;
+};
+
+template <typename Value, typename Parameters>
+struct coordinate_minor_range_generator
+{
+    typedef coordinate_major_cursor<Value, Parameters>   outer_cursor_type;
+    typedef coordinate_minor_cursor<Value, Parameters>   type;
+    static const int                                     level= 2;
+
+    type begin(outer_cursor_type c) const { return type(c.ref, c.offset); }
+    type end(outer_cursor_type c) const { return type(c.ref, c.next_offset); }
+};
+
+template <typename Value, typename Parameters>
+struct coordinate_row_range_generator
+{
+    typedef const coordinate2D<Value, Parameters>&       matrix_ref_type;
+    typedef coordinate_major_cursor<Value, Parameters>   type;
+    typedef complexity_classes::linear_cached            complexity;
+    static const int                                     level= 1;
+
+    type begin(matrix_ref_type A) const { return type(A, 0, 0); }
+    type end(matrix_ref_type A) const { return type(A, num_rows(A), nnz(A)); }
+};
+
+template <typename Value, typename Parameters>
+struct coordinate_col_range_generator
+{
+    typedef const coordinate2D<Value, Parameters>&       matrix_ref_type;
+    typedef coordinate_major_cursor<Value, Parameters>   type;
+    typedef complexity_classes::linear_cached            complexity;
+    static const int                                     level= 1;
+
+    type begin(matrix_ref_type A) const { return type(A, 0, 0); }
+    type end(matrix_ref_type A) const { return type(A, num_cols(A), nnz(A)); }
+};
+
+
+}} // namespace mtl::matrix
+
+
+namespace mtl { namespace traits {
+
+    // Cursor over all rows
+    // Supported if row major matrix
+    template <typename Value, typename Parameters>
+    struct range_generator<glas::tag::row, matrix::coordinate2D<Value, Parameters> >
+      : boost::mpl::if_<
+	    boost::is_same<typename Parameters::orientation, row_major>
+	  , matrix::coordinate_row_range_generator<Value, Parameters>
+ 	  , range_generator<tag::unsupported, matrix::coordinate2D<Value, Parameters> >
+        >::type {};
+
+    template <typename Value, typename Parameters>
+    struct range_generator<glas::tag::col, matrix::coordinate2D<Value, Parameters> >
+      : boost::mpl::if_<
+	    boost::is_same<typename Parameters::orientation, col_major>
+	  , matrix::coordinate_col_range_generator<Value, Parameters>
+ 	  , range_generator<tag::unsupported, matrix::coordinate2D<Value, Parameters> >
+        >::type {};
+
+    template <class Value, class Parameters>
+    struct range_generator<glas::tag::nz, matrix::coordinate_major_cursor<Value, Parameters> >
+      : matrix::coordinate_minor_range_generator<Value, Parameters>
+    {};
+
+
+}} // namespace mtl::traits
+
+namespace mtl {
+	using matrix::coordinate2D;
+}
+
+#endif // MTL_COORDINATE2D_INCLUDE
+

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/crtp_base_matrix.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_CRTP_BASE_MATRIX_INCLUDE
+#define MTL_CRTP_BASE_MATRIX_INCLUDE
+
+#include <iostream>
+#include <algorithm>
+#include <boost/mpl/bool.hpp>
+#include <boost/utility/enable_if.hpp>
+#include <boost/numeric/mtl/operation/print.hpp>
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/operation/matrix_bracket.hpp>
+#include <boost/numeric/mtl/operation/copy.hpp>
+#include <boost/numeric/mtl/operation/mult.hpp>
+#include <boost/numeric/mtl/operation/right_scale_inplace.hpp>
+#include <boost/numeric/mtl/operation/divide_by_inplace.hpp>
+#include <boost/numeric/mtl/matrix/all_mat_expr.hpp>
+#include <boost/numeric/mtl/matrix/diagonal_setup.hpp>
+#include <boost/numeric/mtl/matrix/inserter.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/ashape.hpp>
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/eval_dense.hpp>
+#include <boost/numeric/mtl/utility/irange.hpp>
+#include <boost/numeric/mtl/utility/iset.hpp>
+#include <boost/numeric/mtl/operation/mult_assign_mode.hpp>
+#include <boost/numeric/mtl/operation/compute_factors.hpp>
+#include <boost/numeric/mtl/operation/column_in_matrix.hpp>
+#include <boost/numeric/mtl/operation/row_in_matrix.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+#ifdef MTL_WITH_INITLIST
+# include <initializer_list>
+#endif
+
+namespace mtl { namespace matrix {
+
+template <typename Source, typename Matrix>
+struct crtp_assign
+{
+    Matrix& operator()(const Source& source, Matrix& matrix)
+    {
+	return assign(source, matrix, typename ashape::ashape<Source>::type());
+    }
+private:
+    /// Assign scalar to a matrix by setting the matrix to a multiple of unity matrix
+    /** Uses internally \sa diagonal_setup, for details see there. **/
+    Matrix& assign(const Source& source, Matrix& matrix, ashape::scal)
+    {
+	vampir_trace<3055> tracer;
+	MTL_DEBUG_THROW_IF(num_rows(matrix) * num_cols(matrix) == 0,
+			   range_error("Trying to initialize a 0 by 0 matrix with a value"));
+	diagonal_setup(matrix, source);
+	return matrix;
+    }
+
+    /// Assign matrix expressions by copying except for some special expressions
+    Matrix& assign(const Source& source, Matrix& matrix, typename ashape::ashape<Matrix>::type)
+    {
+	vampir_trace<3056> tracer;
+	// Self-assignment between different types shouldn't happen.
+	matrix.checked_change_resource(source);
+	matrix_copy(source, matrix);
+	return matrix;
+    }
+};
+
+
+
+/// Assign sum by assigning first argument and adding second
+/*  Note that this is more special then assigning arbitrary expressions including matrices itself
+    because mat_mat_plus_expr <E1, E2> is a derived class from mat_expr < MatrixSrc >. **/
+template <typename E1, typename E2, typename Matrix>
+struct crtp_assign<mat_mat_plus_expr<E1, E2>, Matrix>
+{
+    Matrix& operator()(const mat_mat_plus_expr<E1, E2>& src, Matrix& matrix)
+    {
+	vampir_trace<3056> tracer;
+	matrix.checked_change_resource(src.first);
+	matrix= src.first;
+	return matrix+= src.second;
+    }
+};
+
+/// Assign difference by assigning first argument and subtracting second
+/*  Note that this is more special then assigning arbitrary expressions including matrices itself
+    because mat_mat_minus_expr <E1, E2> is a derived class from mat_expr < MatrixSrc >. **/
+template <typename E1, typename E2, typename Matrix>
+struct crtp_assign<mat_mat_minus_expr<E1, E2>, Matrix>
+{
+    Matrix& operator()(const mat_mat_minus_expr<E1, E2>& src, Matrix& matrix)
+    {
+	vampir_trace<3057> tracer;
+	matrix.checked_change_resource(src.first);
+	matrix= src.first;
+	return matrix-= src.second;
+    }
+};
+
+/// Assign product by calling mult
+template <typename E1, typename E2, typename Matrix>
+struct crtp_assign<mat_mat_times_expr<E1, E2>, Matrix>
+{
+    Matrix& operator()(const mat_mat_times_expr<E1, E2>& src, Matrix& matrix)
+    {
+	vampir_trace<4012> tracer;
+	operation::compute_factors<Matrix, mat_mat_times_expr<E1, E2> > factors(src);
+	matrix.checked_change_resource(factors.first, factors.second);
+	mult(factors.first, factors.second, matrix);
+	return matrix;
+    }
+};
+
+
+/// Assign element-wise product
+template <typename E1, typename E2, typename Matrix>
+struct crtp_assign<mat_mat_ele_times_expr<E1, E2>, Matrix>
+{
+    Matrix& operator()(const mat_mat_ele_times_expr<E1, E2>& src, Matrix& matrix)
+    {
+	vampir_trace<3028> tracer;
+	operation::compute_factors<Matrix, mat_mat_ele_times_expr<E1, E2> > factors(src);
+	matrix.checked_change_resource(factors.first);
+	matrix= factors.first;
+	return matrix.ele_rscale(factors.second);
+    }
+};
+
+
+
+/// Assign c-style 2D-array, because it's easier to initialize.
+template <typename Value, unsigned Rows, unsigned Cols, typename Matrix>
+struct crtp_assign<Value[Rows][Cols], Matrix>
+{
+    Matrix& operator()(const Value src[Rows][Cols], Matrix& matrix)
+    {
+	vampir_trace<3059> tracer;
+	typedef typename Collection<Matrix>::size_type size_type;
+
+	matrix.checked_change_dim(Rows, Cols);
+	inserter<Matrix>  ins(matrix, matrix.dim2());
+
+	for (size_type r= 0; r < Rows; ++r)
+	    for (size_type c= 0; c < Cols; ++c)
+		ins(r, c) << src[r][c];
+	return matrix;
+    }
+};
+
+#if defined(MTL_WITH_INITLIST) && defined(MTL_WITH_AUTO) && defined(MTL_WITH_RANGEDFOR)
+    /// Constructor for initializer list \p values
+    template <typename Value2, typename Matrix>
+    struct crtp_assign<std::initializer_list<std::initializer_list<Value2> >, Matrix>
+    {
+	Matrix& operator()(std::initializer_list<std::initializer_list<Value2> > values, Matrix& matrix)
+	{
+	    typedef typename Collection<Matrix>::size_type size_type;
+	    size_type nr= values.size(), nc= nr > 0? values.begin()->size() : 0;
+	    matrix.checked_change_dim(nr, nc);
+	    inserter<Matrix>  ins(matrix, matrix.dim2());
+
+	    size_t r= 0;
+	    for (auto l : values) {
+		size_t c= 0;
+		MTL_THROW_IF(l.size() != nc, logic_error("All sub-lists must have same size!"));
+		for (auto v : l)
+		    ins(r, c++) << v;
+		r++;
+	    }
+	    return matrix;
+	}
+    };
+#endif
+
+
+template <typename Vector, typename Matrix>
+struct crtp_assign<multi_vector<Vector>, Matrix>
+{
+    Matrix& operator()(const multi_vector<Vector>& src, Matrix& matrix)
+    {
+	vampir_trace<3060> tracer;
+	typedef typename Collection<Matrix>::size_type size_type;
+
+	matrix.checked_change_resource(src);
+	// del checked_change_dim(num_rows(src), num_cols(src));
+	inserter<Matrix>  ins(matrix);
+
+	for (size_type r= 0; r < num_rows(src); ++r)
+	    for (size_type c= 0; c < num_cols(src); ++c)
+		ins(r, c) << src[r][c];
+	return matrix;
+    }
+};
+
+
+/// Assign content of a file to the matrix
+template <typename IFStream, typename OFStream, typename Matrix>
+struct crtp_assign<io::matrix_file<IFStream, OFStream>, Matrix>
+{
+    Matrix& operator()(const io::matrix_file<IFStream, OFStream>& file, Matrix& matrix)
+    {
+	vampir_trace<3029> tracer;
+	IFStream stream(file.file_name().c_str());
+	stream >> matrix;
+	return matrix;
+    }
+};
+
+/// Assign-add matrix expressions by incrementally copying except for some special expressions
+template <typename Source, typename Matrix>
+struct crtp_plus_assign
+{
+    Matrix& operator()(const Source& source, Matrix& matrix)
+    {
+	vampir_trace<3030> tracer;
+	return assign(source, matrix, typename ashape::ashape<Source>::type());
+    }
+  private:
+    Matrix& assign(const Source& source, Matrix& matrix, typename ashape::ashape<Matrix>::type)
+    {
+	matrix_copy_plus(source, matrix);
+	return matrix;
+    }
+};
+
+/// Assign-add sum by adding both arguments
+/** Note that this is more special then assigning arbitrary expressions including matrices itself
+	because mat_mat_plus_expr <E1, E2> is a derived class from
+	mat_expr < MatrixSrc >. **/
+template <typename E1, typename E2, typename Matrix>
+struct crtp_plus_assign<mat_mat_plus_expr<E1, E2>, Matrix>
+{
+    Matrix& operator()(const mat_mat_plus_expr<E1, E2>& src, Matrix& matrix)
+    {
+	vampir_trace<3030> tracer;
+	matrix+= src.first;
+	return matrix+= src.second;
+    }
+};
+
+template <typename E1, typename E2, typename Matrix>
+struct crtp_plus_assign<mat_mat_minus_expr<E1, E2>, Matrix>
+{
+    Matrix& operator()(const mat_mat_minus_expr<E1, E2>& src, Matrix& matrix)
+    {
+	vampir_trace<3030> tracer;
+	matrix+= src.first;
+	return matrix-= src.second;
+    }
+};
+
+template <typename E1, typename E2, typename Matrix>
+struct crtp_plus_assign<mat_mat_ele_times_expr<E1, E2>, Matrix>
+{
+    Matrix& operator()(const mat_mat_ele_times_expr<E1, E2>& src, Matrix& matrix)
+    {
+	vampir_trace<3030> tracer;
+	Matrix Prod(ele_prod(src.first, src.second));
+	return matrix+= Prod;
+    }
+};
+
+template <typename E1, typename E2, typename Matrix>
+struct crtp_plus_assign<mat_mat_times_expr<E1, E2>, Matrix>
+{
+    Matrix& operator()(const mat_mat_times_expr<E1, E2>& src, Matrix& matrix)
+    {
+	vampir_trace<3030> tracer;
+	operation::compute_factors<Matrix, mat_mat_times_expr<E1, E2> > factors(src);
+	gen_mult(factors.first, factors.second, matrix, assign::plus_sum(),
+		 tag::flat<tag::matrix>(), tag::flat<tag::matrix>(), tag::flat<tag::matrix>());
+	return matrix;
+    }
+};
+
+
+/// Assign-subtract matrix expressions by decrementally copying except for some special expressions
+template <typename Source, typename Matrix>
+struct crtp_minus_assign
+{
+    Matrix& operator()(const Source& source, Matrix& matrix)
+    {
+	vampir_trace<3031> tracer;
+	return assign(source, matrix, typename ashape::ashape<Source>::type());
+    }
+private:
+    Matrix& assign(const Source& source, Matrix& matrix, typename ashape::ashape<Matrix>::type)
+    {
+	matrix_copy_minus(source, matrix);
+	return matrix;
+    }
+};
+
+/// Assign-subtract sum by adding both arguments
+/** Note that this is more special then assigning arbitrary expressions including matrices itself
+	because mat_mat_plus_expr <E1, E2> is a derived class from
+	mat_expr < MatrixSrc >. **/
+template <typename E1, typename E2, typename Matrix>
+struct crtp_minus_assign<mat_mat_plus_expr<E1, E2>, Matrix>
+{
+    Matrix& operator()(const mat_mat_plus_expr<E1, E2>& src, Matrix& matrix)
+    {
+	vampir_trace<3031> tracer;
+	matrix-= src.first;
+	return matrix-= src.second;
+    }
+};
+
+/// Assign-subtracting difference by subtracting first argument and adding the second one
+/** Note that this is more special then assigning arbitrary expressions including matrices itself
+	because mat_mat_minus_expr <E1, E2> is a derived class from
+	mat_expr < MatrixSrc >. **/
+template <typename E1, typename E2, typename Matrix>
+struct crtp_minus_assign<mat_mat_minus_expr<E1, E2>, Matrix>
+{
+    Matrix& operator()(const mat_mat_minus_expr<E1, E2>& src, Matrix& matrix)
+    {
+	vampir_trace<3031> tracer;
+	matrix-= src.first;
+	return matrix+= src.second;
+    }
+};
+
+template <typename E1, typename E2, typename Matrix>
+struct crtp_minus_assign<mat_mat_ele_times_expr<E1, E2>, Matrix>
+{
+    Matrix& operator()(const mat_mat_ele_times_expr<E1, E2>& src, Matrix& matrix)
+    {
+	vampir_trace<3031> tracer;
+	Matrix Prod(ele_prod(src.first, src.second));
+	return matrix-= Prod;
+    }
+};
+
+/// Assign-subtract product by calling gen_mult
+/** Note that this does not work for arbitrary expressions. **/
+template <typename E1, typename E2, typename Matrix>
+struct crtp_minus_assign<mat_mat_times_expr<E1, E2>, Matrix>
+{
+    Matrix& operator()(const mat_mat_times_expr<E1, E2>& src, Matrix& matrix)
+    {
+	vampir_trace<3031> tracer;
+	operation::compute_factors<Matrix, mat_mat_times_expr<E1, E2> > factors(src);
+	gen_mult(factors.first, factors.second, matrix, assign::minus_sum(), tag::flat<tag::matrix>(), tag::flat<tag::matrix>(), tag::flat<tag::matrix>());
+	return matrix;
+    }
+};
+
+
+
+/// Base class to provide matrix assignment operators generically
+template <typename Matrix, typename ValueType, typename SizeType>
+struct crtp_matrix_assign
+{
+private:
+
+    // For (compatible) dense matrices do a loop over all entries
+    template <typename Source>
+    Matrix& density_assign(const Source& src, boost::mpl::true_)
+    {
+	vampir_trace<3032> tracer;
+	// typedef typename Collection<Source>::size_type size_type;
+	typedef unsigned size_type;
+
+	// std::cout << "Dense assignment\n";
+	checked_change_resource(src);
+
+	Matrix& matrix= static_cast<Matrix&>(*this);
+	for (size_type r= 0; r < num_rows(matrix); ++r)
+	    for (size_type c= 0; c < num_cols(matrix); ++c)
+		matrix[r][c]= src[r][c];
+	return matrix;
+    }
+
+    // If sparse matrices are involved evaluate step-wise (or assignment from scalar)
+    template <typename Source>
+    Matrix& density_assign(const Source& src, boost::mpl::false_)
+    {
+	// std::cout << "Sparse assignment\n";
+	return crtp_assign<Source, Matrix>()(src, static_cast<Matrix&>(*this));
+    }
+
+
+    // For (compatible) dense matrices do a loop over all entries
+    template <typename Source>
+    Matrix& density_plus_assign(const Source& src, boost::mpl::true_)
+    {
+	vampir_trace<3033> tracer;
+	// typedef typename Collection<Source>::size_type size_type;
+	typedef unsigned size_type;
+
+	// std::cout << "Dense assignment\n";
+	checked_change_resource(src);
+	// del checked_change_dim(num_rows(src), num_cols(src));
+
+	Matrix& matrix= static_cast<Matrix&>(*this);
+	for (size_type r= 0; r < num_rows(matrix); ++r)
+	    for (size_type c= 0; c < num_cols(matrix); ++c)
+		matrix[r][c]+= src[r][c];
+	return matrix;
+    }
+
+    // If sparse matrices are involved evaluate step-wise (or assignment from scalar)
+    template <typename Source>
+    Matrix& density_plus_assign(const Source& src, boost::mpl::false_)
+    {
+	// std::cout << "Sparse assignment\n";
+	return crtp_plus_assign<Source, Matrix>()(src, static_cast<Matrix&>(*this));
+    }
+
+    // For (compatible) dense matrices do a loop over all entries
+    template <typename Source>
+    Matrix& density_minus_assign(const Source& src, boost::mpl::true_)
+    {
+	vampir_trace<3034> tracer;
+	// typedef typename Collection<Source>::size_type size_type;
+	typedef unsigned size_type;
+
+	// std::cout << "Dense assignment\n";
+	checked_change_resource(src);
+
+	Matrix& matrix= static_cast<Matrix&>(*this);
+	for (size_type r= 0; r < num_rows(matrix); ++r)
+	    for (size_type c= 0; c < num_cols(matrix); ++c)
+		matrix[r][c]-= src[r][c];
+	return matrix;
+    }
+
+    // If sparse matrices are involved evaluate step-wise (or assignment from scalar)
+    template <typename Source>
+    Matrix& density_minus_assign(const Source& src, boost::mpl::false_)
+    {
+	// std::cout << "Sparse assignment\n";
+	return crtp_minus_assign<Source, Matrix>()(src, static_cast<Matrix&>(*this));
+    }
+
+    // For (compatible) dense matrices do a loop over all entries
+    template <typename Source>
+    Matrix& density_ele_rscale(const Source& src, boost::mpl::true_)
+    {
+	vampir_trace<3035> tracer;
+	// typedef typename Collection<Source>::size_type size_type;
+	typedef unsigned size_type;
+
+	// std::cout << "Dense assignment\n";
+	checked_change_resource(src);
+	// del checked_change_dim(num_rows(src), num_cols(src));
+
+	Matrix& matrix= static_cast<Matrix&>(*this);
+	for (size_type r= 0; r < num_rows(matrix); ++r)
+	    for (size_type c= 0; c < num_cols(matrix); ++c)
+		matrix[r][c]*= src[r][c];
+	return matrix;
+    }
+
+    // If sparse matrices are involved evaluate step-wise (or assignment from scalar)
+    template <typename Factor>
+    Matrix& density_ele_rscale(const Factor& alpha, boost::mpl::false_)
+    {
+	// std::cout << "Sparse assignment\n";
+	matrix_copy_ele_times(alpha, static_cast<Matrix&>(*this));
+	return static_cast<Matrix&>(*this);
+    }
+
+  public:
+
+    /// Check wether source and target have compatible resources, generalization of check_dim
+    /** For expressions like A= B + C, A can be set to the size of B and C if still is 0 by 0. **/
+    template <typename Src>
+    void check_resource(const Src& src) const
+    {	check_resource(src, typename mtl::traits::category<Matrix>::type());    }
+
+    // Default case just check_dim
+    template <typename Src>
+    void check_resource(const Src& src, tag::universe) const
+    {	check_dim(num_rows(src), num_cols(src));    }
+
+    /// Check wether source and target have compatible resources and wether target has already resources
+    /** For expressions like A+= B + C, A must be already larger then 0 by 0 and compatible to B and C. **/
+    //  Generalization with 2 arguments might be needed (check rows from first and columns from second)
+    template <typename Src>
+    void check_ready_resource(const Src& src) const
+    {
+	MTL_DEBUG_THROW_IF(num_rows(src) * num_cols(src) == 0, need_nonempty());
+	check_resource(src);
+    }
+
+    /// Check wether source and target have compatible resources and adapt empty target
+    /** For expressions like A= B + C, A can be set to the size of B and C if still is 0 by 0. **/
+    template <typename Src>
+    void checked_change_resource(const Src& src)
+    {	checked_change_resource(src, src);   }
+
+    /// Check whether source and target have compatible resources and adapt empty target
+    /** For expressions like A= B + C, A can be set to the size of B and C if still is 0 by 0. **/
+    template <typename Src1, typename Src2>
+    void checked_change_resource(const Src1& src1, const Src2& src2)
+    {   checked_change_resource_aux(src1, src2, typename mtl::traits::category<Matrix>::type());    }
+
+    template <typename Src1, typename Src2>
+    void checked_change_resource_aux(const Src1& src1, const Src2& src2, tag::universe)
+    {   checked_change_dim(num_rows(src1), num_cols(src2));  }
+
+
+    /// Check whether matrix sizes are compatible or if matrix is 0 by 0 change it to r by c.
+    /** Deprecated, superseded by checked_change_resource. **/
+    void checked_change_dim(SizeType r, SizeType c)
+    {
+	Matrix& matrix= static_cast<Matrix&>(*this);
+	matrix.check_dim(r, c);
+	matrix.change_dim(r, c);
+    }
+
+    /// Templated assignment implemented by functor to allow for partial specialization
+    // Despite there is only an untemplated assignment and despite the disable_if MSVC whines about ambiguity :-!
+    // Scalar assignment is also taking out because it has another return type
+    template <typename Source>
+    typename boost::disable_if_c<boost::is_same<Matrix, Source>::value
+                                   || boost::is_same<typename ashape::ashape<Source>::type, ashape::scal>::value,
+				 Matrix&>::type
+    operator=(const Source& src)
+    {
+	return density_assign(src, boost::mpl::bool_< boost::is_same<typename ashape::ashape<Matrix>::type,
+			                                             typename ashape::ashape<Source>::type>::value
+			                              && mtl::traits::eval_dense< mat_mat_asgn_expr<Matrix, Source> >::value >());
+    }
+
+    // Helper type for assigning scalars to handle both A= a; and A= a, b, c;
+    template <typename Source>
+    struct scalar_assign
+    {
+	scalar_assign(Source src, Matrix& matrix)
+	  : src(src), with_comma(false), r(0), c(0), matrix(matrix), ins(matrix, 1) {}
+
+	~scalar_assign()
+	{
+	    vampir_trace<3047> tracer;
+	    if (with_comma) {
+		MTL_DEBUG_THROW_IF(r != num_rows(matrix), incompatible_size("Not all matrix entries initialized!"));
+	    } else {
+		using std::min;
+		if (src == math::zero(src)) // it is already set to zero
+		    return;
+		// Otherwise set diagonal (if square)
+		for (SizeType i= 0, n= min(num_rows(matrix), num_cols(matrix)); i < n; i++)
+		    ins[i][i] << src;
+	    }
+	}
+
+	template <typename ValueSource>
+	scalar_assign& operator, (ValueSource val)
+	{
+	    if (!with_comma) {
+		with_comma= true;
+		assert(r == 0 && c == 0);
+		ins[r][c++] << src; // We haven't set v[0] yet
+		if (c == num_cols(matrix))
+		    c= 0, r++;
+	    }
+	    ins[r][c++] << val;
+	    if (c == num_cols(matrix))
+		c= 0, r++;
+	    return *this;
+	}
+
+	Source         src;
+	bool           with_comma;
+	SizeType       r, c;
+	Matrix&        matrix;
+	inserter<Matrix> ins;
+    };
+
+    template <typename Source>
+    typename boost::enable_if<boost::is_same<typename ashape::ashape<Source>::type, ashape::scal>,
+			      scalar_assign<Source> >::type
+    operator=(Source src)
+    {
+	Matrix& matrix= static_cast<Matrix&>(*this);
+	MTL_DEBUG_THROW_IF(num_rows(matrix) * num_cols(matrix) == 0,
+			   range_error("Trying to initialize a 0 by 0 matrix with a value"));
+	set_to_zero(matrix);
+	return scalar_assign<Source>(src, static_cast<Matrix&>(*this));
+    }
+
+    template <typename Source>
+    Matrix& operator+=(const Source& src)
+    {
+	return density_plus_assign(src, mtl::traits::eval_dense< mat_mat_asgn_expr<Matrix, Source> >());
+    }
+
+    template <typename Source>
+    Matrix& operator-=(const Source& src)
+    {
+	return density_minus_assign(src, mtl::traits::eval_dense< mat_mat_asgn_expr<Matrix, Source> >());
+    }
+
+    /// Scale matrix (in place) with scalar value or other matrix
+    template <typename Factor>
+    Matrix& operator*=(const Factor& alpha)
+    {
+	right_scale_inplace(static_cast<Matrix&>(*this), alpha);
+	return static_cast<Matrix&>(*this);
+    }
+
+    // Element-wise scaling from right (i.e. like *= as elementwise)
+    template <typename Factor>
+    Matrix& ele_rscale(const Factor& alpha)
+    {
+	return density_ele_rscale(alpha, mtl::traits::eval_dense< mat_mat_asgn_expr<Matrix, Factor> >());
+    }
+
+    /// Divide matrix (in place) by scalar value
+    // added by Hui Li
+    template <typename Factor>
+    Matrix& operator/=(const Factor& alpha)
+    {
+	divide_by_inplace(static_cast<Matrix&>(*this), alpha);
+	return static_cast<Matrix&>(*this);
+    }
+};
+
+
+
+template <typename Matrix, typename ValueType, typename SizeType>
+struct const_crtp_matrix_bracket
+{
+    template <typename T>
+    typename boost::disable_if_c<boost::is_same<T, mtl::irange>::value || boost::is_same<T, mtl::iset>::value,
+				 operations::bracket_proxy<Matrix, const Matrix&, ValueType> >::type
+    operator[] (const T& row) const
+    {
+	return operations::bracket_proxy<Matrix, const Matrix&, ValueType>(static_cast<const Matrix&>(*this), row);
+    }
+
+    // Compiler error (later) if no sub_matrix function (or row vector resp.) available
+    template <typename T>
+    typename boost::enable_if<boost::is_same<T, mtl::irange>, operations::range_bracket_proxy<Matrix, const Matrix&, const Matrix> >::type
+    operator[] (const T& row_range) const
+    {
+	return operations::range_bracket_proxy<Matrix, const Matrix&, const Matrix>(static_cast<const Matrix&>(*this), row_range);
+    }
+
+    operations::set_bracket_proxy<Matrix, const Matrix&, const Matrix>
+    operator[] (const iset& row_set) const
+    {
+	return operations::set_bracket_proxy<Matrix, const Matrix&, const Matrix>(static_cast<const Matrix&>(*this), row_set);
+    }
+};
+
+template <typename Matrix, typename ValueType, typename SizeType>
+struct crtp_matrix_bracket
+{
+    operations::bracket_proxy<Matrix, const Matrix&, const ValueType&>
+    operator[] (SizeType row) const
+    {
+        return operations::bracket_proxy<Matrix, const Matrix&, const ValueType&>(static_cast<const Matrix&>(*this), row);
+    }
+
+    template <typename T>
+    typename boost::disable_if_c<boost::is_same<T, mtl::irange>::value || boost::is_same<T, mtl::iset>::value,
+			       operations::bracket_proxy<Matrix, Matrix&, ValueType&> >::type
+    // operations::bracket_proxy<Matrix, Matrix&, ValueType&>
+    operator[] (const T& row)
+    {
+        return operations::bracket_proxy<Matrix, Matrix&, ValueType&>(static_cast<Matrix&>(*this), row);
+    }
+
+    // Compiler error (later) if no sub_matrix function available
+    operations::range_bracket_proxy<Matrix, const Matrix&, const Matrix>
+    operator[] (const irange& row_range) const
+    {
+	return operations::range_bracket_proxy<Matrix, const Matrix&, const Matrix>(static_cast<const Matrix&>(*this), row_range);
+    }
+
+    // Compiler error (later) if no sub_matrix function available
+    template <typename T>
+    typename boost::enable_if<boost::is_same<T, mtl::irange>, operations::range_bracket_proxy<Matrix, Matrix&, Matrix> >::type
+    // operations::range_bracket_proxy<Matrix, Matrix&, Matrix>
+    operator[] (const T& row_range)
+    {
+	return operations::range_bracket_proxy<Matrix, Matrix&, Matrix>(static_cast<Matrix&>(*this), row_range);
+    }
+
+    operations::set_bracket_proxy<Matrix, const Matrix&, const Matrix>
+    operator[] (const iset& row_set) const
+    {
+	return operations::set_bracket_proxy<Matrix, const Matrix&, const Matrix>(static_cast<const Matrix&>(*this), row_set);
+    }
+};
+
+template <typename Matrix, typename ValueType, typename SizeType>
+struct crtp_matrix_lvalue
+{
+    // Function must be overwritten by Matrix if m(row, col) does not return a reference
+    ValueType& lvalue(SizeType row, SizeType col)
+    {
+	return static_cast<Matrix&>(*this)(row, col);
+    }
+};
+
+template <typename Matrix, typename ValueType, typename SizeType>
+struct const_crtp_base_matrix
+  : public const_crtp_matrix_bracket<Matrix, ValueType, SizeType>
+{};
+
+template <typename Matrix, typename ValueType, typename SizeType>
+struct mutable_crtp_base_matrix
+  : public crtp_matrix_bracket<Matrix, ValueType, SizeType>,
+    public crtp_matrix_assign<Matrix, ValueType, SizeType>
+{};
+
+template <typename Matrix, typename ValueType, typename SizeType>
+struct crtp_base_matrix
+  : boost::mpl::if_<boost::is_const<Matrix>,
+		    const_crtp_base_matrix<Matrix, ValueType, SizeType>,
+		    mutable_crtp_base_matrix<Matrix, ValueType, SizeType>
+                   >::type
+{};
+
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_CRTP_BASE_MATRIX_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/dense2D.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_DENSE2D_INCLUDE
+#define MTL_DENSE2D_INCLUDE
+
+
+#include <algorithm>
+#include <boost/mpl/bool.hpp>
+#include <boost/type_traits/is_same.hpp>
+#include <boost/mpl/if.hpp>
+#include <boost/utility/enable_if.hpp>
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/matrix/crtp_base_matrix.hpp>
+#include <boost/numeric/mtl/matrix/base_sub_matrix.hpp>
+#include <boost/numeric/mtl/matrix/all_mat_expr.hpp>
+#include <boost/numeric/mtl/matrix/operators.hpp>
+#include <boost/numeric/mtl/detail/contiguous_memory_block.hpp>
+#include <boost/numeric/mtl/operation/set_to_zero.hpp>
+#include <boost/numeric/mtl/operation/compute_factors.hpp>
+#include <boost/numeric/mtl/operation/clone.hpp>
+#include <boost/numeric/mtl/operation/is_negative.hpp>
+#include <boost/numeric/mtl/utility/common_include.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/is_static.hpp>
+#include <boost/numeric/mtl/utility/irange.hpp>
+#include <boost/numeric/mtl/utility/dense_el_cursor.hpp>
+#include <boost/numeric/mtl/utility/static_assert.hpp>
+#include <boost/numeric/mtl/utility/strided_dense_el_cursor.hpp>
+#include <boost/numeric/mtl/utility/strided_dense_el_iterator.hpp>
+#include <boost/numeric/mtl/utility/transposed_orientation.hpp>
+#include <boost/numeric/linear_algebra/identity.hpp>
+
+#ifdef MTL_WITH_INITLIST
+# include <initializer_list>
+#endif
+
+// Forward declaration (for friend declaration)
+namespace mtl { namespace traits { namespace detail {
+    template <typename, typename, bool> struct dense2D_iterator_range_generator;
+}}}
+
+namespace mtl { namespace matrix {
+
+
+using std::size_t;
+
+// Forward declarations
+template <typename Value, typename Parameters> class dense2D;
+class dense2D_indexer;
+
+// Helper type
+struct dense2D_sub_ctor {};
+
+// Indexing for dense matrices
+class dense2D_indexer
+{
+    // helpers for public functions
+    size_t offset(size_t ldim, size_t r, size_t c, row_major) const
+    {
+	return r * ldim + c;
+    }
+    size_t offset(size_t ldim, size_t r, size_t c, col_major) const
+    {
+	return c * ldim + r;
+    }
+
+    size_t row(size_t offset, size_t ldim, row_major) const
+    {
+	return offset / ldim;
+    }
+    size_t row(size_t offset, size_t ldim, col_major) const
+    {
+	return offset % ldim;
+    }
+
+    size_t col(size_t offset, size_t ldim, row_major) const
+    {
+	return offset % ldim;
+    }
+    size_t col(size_t offset, size_t ldim, col_major) const
+    {
+	return offset / ldim;
+    }
+
+ public:
+    template <typename Value, class Parameters>
+    size_t operator() (const dense2D<Value, Parameters>& ma, size_t r, size_t c) const
+    {
+	typedef dense2D<Value, Parameters> matrix_type;
+	// convert into c indices
+	typename matrix_type::index_type my_index;
+	size_t my_r= index::change_from(my_index, r);
+	size_t my_c= index::change_from(my_index, c);
+	return offset(ma.ldim, my_r, my_c, typename matrix_type::orientation());
+    }
+
+    template <typename Value, class Parameters>
+    size_t row(const dense2D<Value, Parameters>& ma,
+	       typename dense2D<Value, Parameters>::key_type key) const
+    {
+	typedef dense2D<Value, Parameters> matrix_type;
+	// row with c-index for my orientation
+	size_t r= row(ma.offset(key), ma.ldim, typename matrix_type::orientation());
+	return index::change_to(typename matrix_type::index_type(), r);
+    }
+
+    template <typename Value, class Parameters>
+    size_t col(const dense2D<Value, Parameters>& ma,
+	       typename dense2D<Value, Parameters>::key_type key) const
+    {
+	typedef dense2D<Value, Parameters> matrix_type;
+	// column with c-index for my orientation
+	size_t c= col(ma.offset(key), ma.ldim, typename matrix_type::orientation());
+	return index::change_to(typename matrix_type::index_type(), c);
+    }
+    template <typename, typename> friend class dense2D;
+}; // dense2D_indexer
+
+
+namespace detail
+{
+
+    // Compute required memory
+    // Enabling mechanism to make sure that computation is valid
+    template <typename Parameters, bool Enable>
+    struct dense2D_array_size {
+	static std::size_t const value= 0;
+    };
+
+    template <typename Parameters>
+    struct dense2D_array_size<Parameters, true>
+    {
+	typedef typename Parameters::dimensions   dimensions;
+	MTL_STATIC_ASSERT((dimensions::is_static), "Size must be known at compile time.");
+	static std::size_t const value= dimensions::Num_Rows * dimensions::Num_Cols;
+    };
+
+    // return const-ref if matrix on stack and type itself if on heap
+    template <typename Matrix, bool on_stack>
+    struct ref_on_stack
+    {
+	typedef Matrix type;
+    };
+
+    template <typename Matrix>
+    struct ref_on_stack<Matrix, true>
+    {
+	typedef const Matrix& type;
+    };
+
+} // namespace detail
+
+
+/// Dense matrix type
+template <typename Value, typename Parameters = parameters<> >
+class dense2D
+  : public base_sub_matrix<Value, Parameters>,
+    public mtl::detail::contiguous_memory_block< Value, Parameters::on_stack,
+						 detail::dense2D_array_size<Parameters, Parameters::on_stack>::value >,
+    public crtp_base_matrix< dense2D<Value, Parameters>, Value, std::size_t >,
+    public mat_expr< dense2D<Value, Parameters> >
+{
+    typedef dense2D                                           self;
+    typedef base_sub_matrix<Value, Parameters>                super;
+    typedef mtl::detail::contiguous_memory_block<Value, Parameters::on_stack,
+						 detail::dense2D_array_size<Parameters, Parameters::on_stack>::value>     memory_base;
+    typedef mat_expr< dense2D<Value, Parameters> >            expr_base;
+    typedef crtp_base_matrix< self, Value, std::size_t >      crtp_base;
+    typedef crtp_matrix_assign< self, Value, std::size_t >    assign_base;
+  public:
+    typedef Parameters                        parameters;
+    typedef typename Parameters::orientation  orientation;
+    typedef typename Parameters::index        index_type;
+    typedef typename Parameters::dimensions   dim_type;
+    typedef Value                             value_type;
+    typedef const value_type&                 const_reference;
+    typedef value_type&                       reference;
+
+    typedef const value_type*                 const_pointer_type;
+    typedef const_pointer_type                key_type;
+    typedef std::size_t                       size_type;
+    typedef dense_el_cursor<Value>            el_cursor_type;
+    typedef dense2D_indexer                   indexer_type;
+
+    // Self-similar type unless dimension is fixed
+    // Not supported for the moment
+    typedef self                              sub_matrix_type;
+
+  protected:
+    // Obviously, the next 3 functions must be called after setting dimensions
+    void set_nnz() { this->my_nnz = this->num_rows() * this->num_cols(); }
+    void set_ldim(row_major) { ldim= this->num_cols(); }
+    void set_ldim(col_major) { ldim= this->num_rows(); }
+    void set_ldim() { set_ldim(orientation()); }
+
+    void init()
+    {
+      set_nnz(); set_ldim(); // set_to_zero(*this);
+    }
+
+  public:
+    /// Default constructor, if compile time matrix size allocate memory
+    dense2D() : super(), memory_base(dim_type().num_rows() * dim_type().num_cols())
+    { 	init();     }
+
+    /// Constructor that only sets dimensions, only for run-time dimensions
+    explicit dense2D(mtl::non_fixed::dimensions d)
+	: super(d), memory_base(d.num_rows() * d.num_cols())
+    { 	init();     }
+
+    /// Most common constructor from number of rows and columns
+    explicit dense2D(size_type num_rows, size_type num_cols)
+	: super(dim_type(num_rows, num_cols)),
+	  memory_base(num_rows * num_cols)
+    { 	init();     }
+
+    /// Constructor that sets dimensions and pointer to external data
+    explicit dense2D(mtl::non_fixed::dimensions d, value_type* a)
+      : super(d), memory_base(a, d.num_rows() * d.num_cols())
+    { 	init();     }
+
+    /// Constructor that sets dimensions and pointer to external data
+    explicit dense2D(size_type num_rows, size_type num_cols, value_type* a)
+      : super(mtl::non_fixed::dimensions(num_rows, num_cols)), memory_base(a, num_rows * num_cols)
+    { 	init();     }
+
+    /// Constructor for compile time matrix size
+    /** sets dimensions and pointer to external data **/
+    explicit dense2D(value_type* a)
+	: super(), memory_base(a, dim_type().num_rows() * dim_type().num_cols())
+    {
+	MTL_STATIC_ASSERT((dim_type::is_static), "Size must be known at compile time.");
+        init();
+    }
+
+    /// Default copy constructor
+    dense2D(const self& m)
+	: super(dim_type(m.num_rows(), m.num_cols())),
+	  memory_base(m)
+    {
+	// In case of sub-matrices we need m's ldim -> init doesn't work
+	this->my_nnz= m.my_nnz; ldim= m.ldim;
+    }
+
+    /// Clone constructor, copies every source including sub-matrices and other matrices with references
+    explicit dense2D(const self& m, clone_ctor)
+	: super(mtl::non_fixed::dimensions(m.num_rows(), m.num_cols())),
+	  memory_base(m, clone_ctor())
+    {
+	init();
+	*this= m;
+    }
+
+    /// General copy constructor, uses functionality from CRTP base
+    template <typename MatrixSrc>
+    explicit dense2D(const MatrixSrc& src)
+	: super(), memory_base(dim_type().num_rows() * dim_type().num_cols())
+    {
+	init();
+	*this= src;
+    }
+
+    /// Constructor for creating sub-matrices
+    template <typename MatrixSrc>
+    dense2D(MatrixSrc& matrix, dense2D_sub_ctor,
+	    size_type begin_r, size_type end_r, size_type begin_c, size_type end_c)
+      : super(mtl::non_fixed::dimensions(matrix.num_rows(), matrix.num_cols())),
+        memory_base(matrix.data, (end_r - begin_r) * (end_c - begin_c), true)
+    {
+	sub_matrix_constructor(matrix, begin_r, end_r, begin_c, end_c, boost::mpl::bool_<memory_base::on_stack>());
+    }
+
+#if defined(MTL_WITH_INITLIST) && defined(MTL_WITH_AUTO) && defined(MTL_WITH_RANGEDFOR)
+    /// Constructor for initializer list \p values
+    template <typename Value2>
+    dense2D(std::initializer_list<std::initializer_list<Value2> > values)
+      : super(mtl::non_fixed::dimensions(values.size(), values.size()? values.begin()->size() : 0)),
+    	memory_base(this->num_rows() * this->num_cols())
+    {
+    	init();
+	*this= values;
+    }
+#endif
+
+  private:
+    template <typename MatrixSrc>
+    void sub_matrix_constructor(MatrixSrc& matrix, size_type begin_r, size_type end_r,
+				size_type begin_c, size_type end_c, boost::mpl::false_)
+    {
+	matrix.check_ranges(begin_r, end_r, begin_c, end_c);
+
+	if(end_r <= begin_r || end_c <= begin_c)
+	    set_ranges(0, 0);
+	else {
+	    // Leading dimension doesn't change
+	    this->data += matrix.indexer(matrix, begin_r, begin_c);  // Takes care of indexing
+	    set_ranges(end_r - begin_r, end_c - begin_c);
+	}
+	this->my_nnz= matrix.nnz(); ldim= matrix.get_ldim();
+    }
+
+    template <typename MatrixSrc>
+    void sub_matrix_constructor(MatrixSrc&, size_type, size_type,
+				size_type, size_type, boost::mpl::true_)
+    {
+	MTL_THROW(logic_error("Matrices cannot be used as sub-matrices!"));
+    }
+
+  public:
+
+#ifdef MTL_WITH_MOVE
+    /// Move assignment for data on heap
+    self& operator=(self&& src)
+    {	return self_assign(src, boost::mpl::bool_<memory_base::on_stack>());    }
+
+    // {
+    // 	swap(*this, src);
+    // 	std::cout << "In dense2D::move_assignment\n";
+    // 	return *this;
+    // }
+
+     /// (Copy) Assignment
+    self& operator=(const self& src)
+    {	return self_assign(src, boost::mpl::true_());    }
+
+#else
+    /// (Copy) Assignment
+    self& operator=(typename detail::ref_on_stack<self, memory_base::on_stack>::type src)
+    {
+	return self_assign(src, boost::mpl::bool_<memory_base::on_stack>());
+    }
+#endif
+
+  private:
+    // Already copied for lvalues -> data can be stolen (need non-const ref)
+    self& self_assign(self& src, boost::mpl::false_)
+    {
+	// Self-copy would be an indication of an error
+	assert(this != &src);
+	// std::cout << "In move assignment: this* = \n" << *this << "src = \n" << src;
+
+	this->checked_change_dim(src.num_rows(), src.num_cols());
+	if (this->category == memory_base::view || src.category == memory_base::view)
+	    matrix_copy(src, *this);
+	else {
+	    if (this->num_rows() != src.num_rows() || this->num_cols() != src.num_cols()) {
+		super::change_dim(src.num_rows(), src.num_cols());
+		init();
+	    }
+	    memory_base::move_assignment(src);
+	}
+	// std::cout << "End of move assignment: this* = \n" << *this;
+	return *this;
+    }
+
+    // For matrices with data on stack (or lvalues in C++11)
+    self& self_assign(const self& src, boost::mpl::true_)
+    {
+	if (this != &src) {
+	    this->checked_change_dim(src.num_rows(), src.num_cols());
+	    matrix_copy(src, *this);
+	}
+	return *this;
+    }
+  public:
+
+    // import operators from CRTP base class
+#if 0 // def __PGI
+    using crtp_base::operator=;
+#else
+    using assign_base::operator=;
+#endif
+
+    /// Change dimension, can keep old data
+    void change_dim(size_type r, size_type c, bool keep_data = false)
+    {
+	change_dim(r, c, keep_data, mtl::traits::is_static<self>());
+    }
+
+  private:
+    void change_dim(size_type r, size_type c, bool keep_data, boost::mpl::false_)
+    {
+	if (r == this->num_rows() && c == this->num_cols())
+	    return;
+
+	self temp;
+	if (keep_data) {
+	    temp.super::change_dim(this->num_rows(), this->num_cols());
+	    temp.init();
+	    temp.memory_base::move_assignment(*this);
+	}
+	memory_base::realloc(r*c);
+	super::change_dim(r, c);
+	init();
+	if (keep_data) {
+	    if (r > temp.num_rows() || c > temp.num_cols()){
+		set_to_zero(*this);
+#if 0
+		irange rr(0, std::min(r,temp.num_rows())), cr(0, std::min(c,temp.num_cols()));
+		*this[rr][cr]= temp[rr][cr];
+#endif
+		sub_matrix(*this,0,std::min(r,temp.num_rows()),0,std::min(c,temp.num_cols()))
+		    = sub_matrix(temp,0,std::min(r,temp.num_rows()),0,std::min(c,temp.num_cols()));
+	    } else
+		*this = temp[irange(0, r)][irange(0, c)];
+	}
+    }
+
+    void change_dim(size_type MTL_DEBUG_ARG(r), size_type MTL_DEBUG_ARG(c), bool, boost::mpl::true_)
+    {	assert(r == this->num_rows() && c == this->num_cols());    }
+
+ public:
+    /// Check whether indices r and c are in range
+    bool check_indices(size_t r, size_t c) const
+    {	return r >= this->begin_row() && r < this->end_row() && c >= this->begin_col() && c < this->end_col();    }
+
+    /// Constant access to element
+    const_reference operator() (size_t r, size_t c) const
+    {
+	MTL_DEBUG_THROW_IF(is_negative(r) || r >= this->num_rows() || is_negative(c) || c >= this->num_cols(), index_out_of_range());
+        return this->data[indexer(*this, r, c)];
+    }
+
+    /// Mutable access to element
+    value_type& operator() (size_t r, size_t c)
+    {
+	MTL_DEBUG_THROW_IF(is_negative(r) || r >= this->num_rows() || is_negative(c) || c >= this->num_cols(), index_out_of_range());
+	return this->data[indexer(*this, r, c)];
+    }
+
+    // offset regarding c-style indices
+    size_t c_offset(size_t r, size_t c) const
+    {	return indexer.offset(ldim, r, c, orientation());    }
+
+    /// Get lower dimension [advanced]
+    size_type get_ldim() const
+    {	return ldim;    }
+
+    /// Swap two matrices
+    friend void swap(self& matrix1, self& matrix2)
+    {
+	swap(static_cast<memory_base&>(matrix1), static_cast<memory_base&>(matrix2));
+	swap(static_cast<super&>(matrix1), static_cast<super&>(matrix2));
+	std::swap(matrix1.ldim, matrix2.ldim);
+    }
+
+    void crop() {} ///< Delete structural zeros; only dummy here
+
+    /// Address of first data entry (mutable); to be used with care. [advanced]
+    value_type* address_data() { return this->data; }
+    /// Address of first data entry (constant); to be used with care. [advanced]
+    const value_type* address_data() const { return this->data; }
+
+    /// Whether data is stored in strides
+    bool has_strided_data() const { return this->category != this->own; }
+
+  protected:
+
+    // Set ranges from begin_r to end_r and begin_c to end_c
+    void set_ranges(size_type begin_r, size_type end_r, size_type begin_c, size_type end_c)
+    {
+	super::set_ranges(begin_r, end_r, begin_c, end_c);
+	set_nnz();
+    }
+
+    // Set ranges to a num_row x num_col matrix, keeps indexing
+    void set_ranges(size_type num_rows, size_type num_cols)
+    {
+	set_ranges(this->begin_row(), this->begin_row() + num_rows,
+		   this->begin_col(), this->begin_col() + num_cols);
+    }
+
+  public:
+
+    indexer_type  indexer;
+
+    friend class dense2D_indexer;
+
+#if !defined(_MSC_VER) || _MSC_VER != 1400 // Bug in MSVC 2005
+    template <typename> friend struct sub_matrix_t;
+    template <typename, typename> friend struct mtl::traits::range_generator;
+    template <typename, typename, bool> friend struct mtl::traits::detail::dense2D_iterator_range_generator;
+
+  protected:
+#endif
+
+    // Leading dimension is minor dimension in original matrix
+    // Opposed to other dims doesn't change in sub-matrices
+    size_type     ldim;
+
+}; // dense2D
+
+
+// ================
+// Free functions
+// ================
+
+
+/// Number of rows
+template <typename Value, typename Parameters>
+typename dense2D<Value, Parameters>::size_type
+inline num_rows(const dense2D<Value, Parameters>& matrix)
+{
+    return matrix.num_rows();
+}
+
+/// Number of columns
+template <typename Value, typename Parameters>
+typename dense2D<Value, Parameters>::size_type
+inline num_cols(const dense2D<Value, Parameters>& matrix)
+{
+    return matrix.num_cols();
+}
+
+/// Size of the matrix, i.e. the number of row times columns
+template <typename Value, typename Parameters>
+typename dense2D<Value, Parameters>::size_type
+inline size(const dense2D<Value, Parameters>& matrix)
+{
+    return matrix.num_cols() * matrix.num_rows();
+}
+
+
+}
+
+using matrix::dense2D;
+
+} // namespace mtl::matrix
+
+
+namespace mtl { namespace traits {
+
+
+    // VC 8.0 finds ambiguity with mtl::tag::dense2D (I wonder why, especially here)
+    using mtl::matrix::dense2D;
+
+    // ================
+    // Range generators
+    // For cursors
+    // ================
+
+    template <typename Value, typename Parameters>
+    struct range_generator<glas::tag::all, dense2D<Value, Parameters> >
+      : detail::dense_element_range_generator<dense2D<Value, Parameters>,
+					      dense_el_cursor<Value>, complexity_classes::linear_cached>
+    {};
+
+    template <typename Value, typename Parameters>
+    struct range_generator<glas::tag::nz, dense2D<Value, Parameters> >
+      : detail::dense_element_range_generator<dense2D<Value, Parameters>,
+					      dense_el_cursor<Value>, complexity_classes::linear_cached>
+    {};
+
+    namespace detail
+    {
+	// complexity of dense row cursor depends on storage scheme
+	// if orientation is row_major then complexity is cached_linear, otherwise linear
+	template <typename Orientation> struct dense2D_rc {};
+	template<> struct dense2D_rc<row_major>
+	{
+	    typedef complexity_classes::linear_cached type;
+	};
+	template<> struct dense2D_rc<col_major>
+	{
+	    typedef complexity_classes::linear type;
+	};
+
+	// Complexity of column cursor is of course opposite
+	template <typename Orientation> struct dense2D_cc
+	    : dense2D_rc<typename mtl::traits::transposed_orientation<Orientation>::type>
+	{};
+    }
+
+    template <typename Value, typename Parameters>
+    struct range_generator<glas::tag::row, dense2D<Value, Parameters> >
+	: detail::all_rows_range_generator<dense2D<Value, Parameters>,
+					   typename detail::dense2D_rc<typename Parameters::orientation>::type>
+    {};
+
+    // For a cursor pointing to some row give the range of elements in this row
+    template <typename Value, typename Parameters>
+    struct range_generator<glas::tag::nz,
+			   detail::sub_matrix_cursor<dense2D<Value, Parameters>, glas::tag::row, 2> >
+    {
+	typedef dense2D<Value, Parameters>                                            matrix;
+	typedef typename matrix::size_type                                            size_type;
+	typedef detail::sub_matrix_cursor<matrix, glas::tag::row, 2>               cursor;
+
+	// linear for col_major and linear_cached for row_major
+	typedef typename detail::dense2D_rc<typename Parameters::orientation>::type   complexity;
+	static int const                                                              level = 1;
+
+	typedef typename boost::mpl::if_<
+	    boost::is_same<typename Parameters::orientation, row_major>
+	  , dense_el_cursor<Value>
+	  , strided_dense_el_cursor<Value>
+	>::type type;
+
+      private:
+
+	type dispatch(cursor const& c, size_type col, row_major) const
+	{
+	    return type(c.ref, c.key, col);
+	}
+	type dispatch(cursor const& c, size_type col, col_major) const
+	{
+	    return type(c.ref, c.key, col, c.ref.ldim);
+	}
+
+      public:
+
+	type begin(cursor const& c) const
+	{
+	    return dispatch(c, c.ref.begin_col(), typename matrix::orientation());
+	}
+	type end(cursor const& c) const
+	{
+	    return dispatch(c, c.ref.end_col(), typename matrix::orientation());
+	}
+	type lower_bound(cursor const& c, size_type position) const
+	{
+	    return dispatch(c, std::min(c.ref.end_col(), position), typename matrix::orientation());
+	}
+    };
+
+    template <typename Value, typename Parameters>
+    struct range_generator<glas::tag::all,
+			   detail::sub_matrix_cursor<dense2D<Value, Parameters>, glas::tag::row, 2> >
+        : range_generator<glas::tag::nz,
+			  detail::sub_matrix_cursor<dense2D<Value, Parameters>, glas::tag::row, 2> >
+    {};
+
+
+    template <typename Value, typename Parameters>
+    struct range_generator<glas::tag::col, dense2D<Value, Parameters> >
+	: detail::all_cols_range_generator<dense2D<Value, Parameters>,
+					   typename detail::dense2D_cc<typename Parameters::orientation>::type>
+    {};
+
+    // For a cursor pointing to some row give the range of elements in this row
+    template <typename Value, typename Parameters>
+    struct range_generator<glas::tag::nz,
+			   detail::sub_matrix_cursor<dense2D<Value, Parameters>, glas::tag::col, 2> >
+    {
+	typedef dense2D<Value, Parameters>                                            matrix;
+	typedef typename matrix::size_type                                            size_type;
+	typedef detail::sub_matrix_cursor<matrix, glas::tag::col, 2>               cursor;
+	typedef typename detail::dense2D_cc<typename Parameters::orientation>::type   complexity;
+	static int const                                                              level = 1;
+
+	typedef typename boost::mpl::if_<
+	    boost::is_same<typename Parameters::orientation, col_major>
+	  , dense_el_cursor<Value>
+	  , strided_dense_el_cursor<Value>
+	>::type type;
+
+      private:
+	type dispatch(cursor const& c, size_type row, col_major) const
+	{
+	    return type(c.ref, row, c.key);
+	}
+	type dispatch(cursor const& c, size_type row, row_major) const
+	{
+	    return type(c.ref, row, c.key, c.ref.ldim);
+	}
+
+      public:
+	type begin(cursor const& c) const
+	{
+	    return dispatch(c, c.ref.begin_row(), typename matrix::orientation());
+	}
+	type end(cursor const& c) const
+	{
+	    return dispatch(c, c.ref.end_row(), typename matrix::orientation());
+	}
+	type lower_bound(cursor const& c, size_type position) const
+	{
+	    return dispatch(c, std::min(c.ref.end_row(), position), typename matrix::orientation());
+	}
+    };
+
+    template <typename Value, typename Parameters>
+    struct range_generator<glas::tag::all,
+			   detail::sub_matrix_cursor<dense2D<Value, Parameters>, glas::tag::col, 2> >
+      : public range_generator<glas::tag::nz,
+			       detail::sub_matrix_cursor<dense2D<Value, Parameters>, glas::tag::col, 2> >
+    {};
+
+// =============
+// For iterators
+// =============
+
+
+    namespace detail {
+
+        // Traversal along major dimension first and then along minor
+        template <typename OuterTag, typename Orientation>
+        struct major_traversal
+        {
+	    static const bool value= false;
+        };
+
+        template <> struct major_traversal<glas::tag::row, row_major>
+        {
+	    static const bool value= true;
+        };
+
+        template <> struct major_traversal<glas::tag::col, col_major>
+        {
+	    static const bool value= true;
+        };
+
+
+        template <typename OuterTag, typename Matrix, bool is_const>
+        struct dense2D_iterator_range_generator
+        {
+	    typedef Matrix                                                                matrix_type;
+	    typedef typename matrix_type::size_type                                       size_type;
+	    typedef typename matrix_type::value_type                                      value_type;
+	    typedef typename matrix_type::parameters                                      parameters;
+	    typedef detail::sub_matrix_cursor<matrix_type, OuterTag, 2>                   cursor;
+
+	    // if traverse first along major dimension then memory access is contiguous (otherwise strided)
+	    typedef typename boost::mpl::if_<
+		major_traversal<OuterTag, typename parameters::orientation>
+	      , complexity_classes::linear_cached
+	      , complexity_classes::linear
+	    >::type                                                                       complexity;
+	    static int const                                                              level = 1;
+
+	    // if traverse first along major dimension use pointer otherwise strided iterator
+	    typedef typename boost::mpl::if_<
+		major_traversal<OuterTag, typename parameters::orientation>
+	      , typename boost::mpl::if_c<
+    	            is_const
+		  , const value_type*
+    	          , value_type*
+		>::type
+	      , typename boost::mpl::if_c<
+    	            is_const
+		  , strided_dense_el_const_iterator<value_type>
+    	          , strided_dense_el_iterator<value_type>
+    	        >::type
+	    >::type type;
+
+        private:
+	    // if traverse first along major dim. then return address as pointer
+	    type dispatch(cursor const& c, size_type row, size_type col, complexity_classes::linear_cached) const
+	    {
+		matrix_type& ma= const_cast<matrix_type&>(c.ref);
+		return ma.elements() + ma.indexer(ma, row, col); // &ref[row][col];
+	    }
+
+	    // otherwise strided
+	    type dispatch(cursor const& c, size_type row, size_type col, complexity_classes::linear) const
+	    {
+		// cast const away (is dirty and should be improved later (cursors must distinct constness))
+		matrix_type& ref= const_cast<matrix_type&>(c.ref);
+		return type(ref, row, col, ref.ldim);
+	    }
+
+	    type begin_dispatch(cursor const& c, glas::tag::row) const
+	    {
+		return dispatch(c, c.key, c.ref.begin_col(), complexity());
+	    }
+
+	    type end_dispatch(cursor const& c, glas::tag::row) const
+	    {
+		return dispatch(c, c.key, c.ref.end_col(), complexity());
+	    }
+
+
+	    type begin_dispatch(cursor const& c, glas::tag::col) const
+	    {
+		return dispatch(c, c.ref.begin_row(), c.key, complexity());
+	    }
+
+	    type end_dispatch(cursor const& c, glas::tag::col) const
+	    {
+		return dispatch(c, c.ref.end_row(), c.key, complexity());
+	    }
+
+        public:
+
+	    type begin(cursor const& c) const
+	    {
+		return begin_dispatch(c, OuterTag());
+	    }
+
+	    type end(cursor const& c) const
+	    {
+		return end_dispatch(c, OuterTag());
+	    }
+        };
+
+    } // namespace detail
+
+
+    template <typename Value, typename Parameters, typename OuterTag>
+    struct range_generator<tag::iter::nz,
+			   detail::sub_matrix_cursor<dense2D<Value, Parameters>, OuterTag, 2> >
+      : public detail::dense2D_iterator_range_generator<OuterTag, dense2D<Value, Parameters>, false>
+    {};
+
+    template <typename Value, typename Parameters, typename OuterTag>
+    struct range_generator<tag::iter::all,
+			   detail::sub_matrix_cursor<dense2D<Value, Parameters>, OuterTag, 2> >
+      : public detail::dense2D_iterator_range_generator<OuterTag, dense2D<Value, Parameters>, false>
+    {};
+
+    template <typename Value, typename Parameters, typename OuterTag>
+    struct range_generator<tag::const_iter::nz,
+			   detail::sub_matrix_cursor<dense2D<Value, Parameters>, OuterTag, 2> >
+      : public detail::dense2D_iterator_range_generator<OuterTag, dense2D<Value, Parameters>, true>
+    {};
+
+    template <typename Value, typename Parameters, typename OuterTag>
+    struct range_generator<tag::const_iter::all,
+			   detail::sub_matrix_cursor<dense2D<Value, Parameters>, OuterTag, 2> >
+      : public detail::dense2D_iterator_range_generator<OuterTag, dense2D<Value, Parameters>, true>
+    {};
+
+
+}} // namespace mtl::traits
+
+namespace mtl { namespace matrix {
+
+    // ==========
+    // Sub matrix
+    // ==========
+
+    template <typename Value, typename Parameters>
+    struct sub_matrix_t<dense2D<Value, Parameters> >
+    {
+        typedef dense2D<Value, Parameters>      matrix_type;
+	// copy orientation, ignore index, set dimension to non-fixed and on_stack to false
+	typedef parameters<typename Parameters::orientation> para_type;
+
+        typedef dense2D<Value, para_type>       sub_matrix_type;
+        typedef sub_matrix_type const           const_sub_matrix_type;
+        typedef typename matrix_type::size_type size_type;
+
+        sub_matrix_type operator()(matrix_type& matrix, size_type begin_r, size_type end_r, size_type begin_c, size_type end_c)
+        {
+	    return sub_matrix_type(matrix, dense2D_sub_ctor(), begin_r, end_r, begin_c, end_c);
+        }
+
+        const_sub_matrix_type
+        operator()(matrix_type const& matrix, size_type begin_r, size_type end_r, size_type begin_c, size_type end_c)
+        {
+	    // To minimize code duplication, we use the non-const version
+	    sub_matrix_type tmp((*this)(const_cast<matrix_type&>(matrix), begin_r, end_r, begin_c, end_c));
+	    return tmp;
+        }
+    };
+
+}} // mtl::matrix
+
+namespace mtl {
+
+    // Enable cloning of dense matrices
+    template <typename Value, typename Parameters>
+    struct is_clonable< mtl::matrix::dense2D<Value, Parameters> > : boost::mpl::true_ {};
+
+} // namespace mtl
+
+
+
+namespace math {
+
+    // Multiplicative identities of matrices
+    template <typename Value, typename Parameters>
+    struct identity_t< mult<mtl::matrix::dense2D<Value, Parameters> >, mtl::matrix::dense2D<Value, Parameters> >
+        : public std::binary_function< mult<mtl::matrix::dense2D<Value, Parameters> >,
+				       mtl::matrix::dense2D<Value, Parameters>,
+				       mtl::matrix::dense2D<Value, Parameters> >
+    {
+        typedef mtl::matrix::dense2D<Value, Parameters>  matrix_type;
+
+        matrix_type operator() (const mult<matrix_type>&, const matrix_type& ref) const
+        {
+	    matrix_type tmp(ref);
+	    tmp= one(typename matrix_type::value_type());
+	    return tmp;
+        }
+    };
+
+} // namespace math
+
+
+#endif // MTL_DENSE2D_INCLUDE
+
+
+/*
+Limitations:
+- with compile-time constant dimension, submatrices are not supported (would violate self-similarity)
+- Element cursor doesn't work for sub-matrices (not contiguous)
+*/

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/diagonal_setup.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_DIAGONAL_SETUP_INCLUDE
+#define MTL_DIAGONAL_SETUP_INCLUDE
+
+#include <boost/numeric/mtl/matrix/inserter.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/operation/set_to_zero.hpp>
+#include <boost/numeric/linear_algebra/identity.hpp>
+
+namespace mtl { namespace matrix {
+
+/// Setup a matrix to a multiple of the unity matrix
+/** Intended for sparse matrices but works also with dense matrices.
+    If the value is 0 the matrix is only zeroed out, whereby
+    a sparse matrix will be empty after this operation,
+    i.e. the zeros on the diagonal are not explicitly stored.
+    The diagonal in its generalized form is the set of entries with equal row and column
+    index (since r6843, older revision considered it erroneous to store
+    a non-zero scalar to a non-square matrix).
+ **/
+template <typename Matrix, typename Value>
+inline void diagonal_setup(Matrix& matrix, const Value& value)
+{
+    using std::min;
+    if (num_rows(matrix) == 0 || num_cols(matrix) == 0)
+	return;
+
+    set_to_zero(matrix);
+    inserter<Matrix>      ins(matrix, 1);
+    for (typename Collection<Matrix>::size_type i= 0, n= min(num_rows(matrix), num_cols(matrix)); i < n; ++i)
+	ins[i][i] << value;
+}
+
+}} // namespace mtl::matrix
+
+#endif // MTL_DIAGONAL_SETUP_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/dimension.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_DIMENSIONS_INCLUDE
+#define MTL_DIMENSIONS_INCLUDE
+
+#include <iostream>
+#include <cassert>
+#include <boost/mpl/if.hpp>
+#include <boost/utility/enable_if.hpp>
+
+namespace mtl {
+
+// dimension is a type for declaring matrix dimensions
+// num_rows() and num_cols() return the number or rows and columns
+// is_static says whether it is declared at compile time or not
+
+// Compile time version
+namespace fixed
+{
+
+    /// Compile-time dimensions
+    template <std::size_t Rows, std::size_t Cols>
+    struct dimensions
+    {
+	typedef std::size_t  size_type;
+
+        static size_type const Num_Rows= Rows;
+        static size_type const Num_Cols= Cols;
+
+	// To have the same interface as fixed
+#ifndef NDEBUG
+	/// Constructor does not need arguments but if given they are compared against the template arguments in debug mode
+	explicit dimensions(size_type r= Rows, size_type c= Cols)
+	{
+	    assert(r == Rows); assert(c == Cols);
+	}
+#else
+	explicit dimensions(size_type, size_type) {}
+	explicit dimensions(size_type) {}
+	explicit dimensions() {}
+#endif
+
+	size_type num_rows() const { return Rows; } ///< Number of rows
+	size_type num_cols() const { return Cols; } ///< Number of columns
+
+	/// To check whether dimensions are static
+	static bool const is_static= true;
+
+	/// Transposed dimension (type)
+	typedef dimensions<Cols, Rows> transposed_type;
+	transposed_type transpose() const
+	{
+	    return transposed_type();
+	}
+    };
+
+    /// Output of dimensions
+    template <std::size_t R, std::size_t C>
+    inline std::ostream& operator<< (std::ostream& stream, dimensions<R, C>)
+    {
+	return stream << R << 'x' << C;
+    }
+
+} // namespace fixed
+
+namespace non_fixed
+{
+    /// Run-time dimensions
+    struct dimensions
+    {
+	typedef std::size_t  size_type;
+
+	/// Constructor
+	dimensions(size_type r= 0, size_type c= 0) : r(r), c(c) {}
+
+	/// Assignment
+	dimensions& operator=(const dimensions& x)
+	{
+	    r= x.r; c= x.c; return *this;
+	}
+	size_type num_rows() const { return r; } ///< Number of rows
+	size_type num_cols() const { return c; } ///< Number of columns
+
+	/// Transposed dimension
+	typedef dimensions transposed_type;
+	transposed_type transpose()
+	{
+	    return transposed_type(c, r);
+	}
+
+	/// To check whether dimensions are static
+	static bool const is_static= false;
+    protected:
+	size_type r, c;
+    };
+
+    /// Output of dimensions
+    inline std::ostream& operator<< (std::ostream& stream, dimensions d)
+    {
+	return stream << d.num_rows() << 'x' << d.num_cols();
+    }
+
+} // namespace non_fixed
+
+} // namespace mtl
+
+#endif // MTL_DIMENSIONS_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/element.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+//
+// Algorithm inspired by Nick Vannieuwenhoven, written by Cornelius Steinhardt
+
+
+#ifndef MTL_ELEMENT_INCLUDE
+#define MTL_ELEMENT_INCLUDE
+
+#include <vector>
+#include <set>
+#include <algorithm>
+#include <iostream>
+#include <functional>
+
+#include <boost/unordered_set.hpp>
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/matrix/dense2D.hpp>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+#include <boost/numeric/mtl/utility/make_copy_or_reference.hpp>
+#include <boost/numeric/itl/pc/comparators.hpp>
+
+
+namespace mtl { namespace matrix {
+
+///  A class representing an element with ValType equals the type of the numeric values
+template <typename ValType>
+class element
+{
+
+public:
+    /// The type of this element.
+    typedef element<ValType> element_type;
+
+    /// The value type of the matrix and rhs elements.
+    typedef ValType value_type;
+
+    /// The type of a set of neighbors.
+    typedef std::vector<element_type*, std::allocator<element_type*> > neighbor_collection_type; // trouble with Clang 4.2 w/o allocator
+
+    /// An iterator over the neighbors of this element.
+    typedef typename neighbor_collection_type::iterator neighbor_iterator;
+
+    /// The type of an unordered set of neighbors.
+    typedef typename boost::unordered_set<
+	element_type*,
+	compare::address_hasher<element_type>,
+	compare::address_compare_equal<element_type>
+       > neighbor_set_type;
+
+    /// The type of the iterator over an unorderd set of neighbors.
+    typedef typename neighbor_set_type::iterator neighbor_set_iterator_type;
+
+    /// The type of matrix.
+    typedef mtl::matrix::dense2D<value_type> matrix_type;
+
+    /// The type of the index vector.
+    typedef mtl::vector::dense_vector<int> index_type;
+
+
+/*******************************************************************************
+ * Constructors
+ ******************************************************************************/
+
+
+    /**
+     * Constructs the element using the memory specified by the two references.
+     *
+     * p_indices: 	a reference to the memory where the indices may be stored.
+     * p_values:	a reference to the memory where the values may be stored.
+     */
+  public:
+    element(int p_sequence_number, const index_type& p_indices, const matrix_type& p_values)
+      : m_indices(p_indices),
+	m_values(p_values),
+	m_sequence_number(p_sequence_number)
+    {}
+
+    element()
+      :	m_sequence_number(-1) {}
+
+    element(const element_type& other)
+      :	m_sequence_number(-1)
+    {	*this = other;	}
+
+    /// Deep copy the given element.
+    void operator=(const element_type& other)
+    {
+	m_sequence_number = other.m_sequence_number;
+	m_neighbors = other.m_neighbors;
+
+	m_indices = other.m_indices;
+	m_values =  other.m_values;
+    }
+
+    /// Returns the unique identifier of this element.
+    inline int get_id() const {	return m_sequence_number; }
+    inline int& get_id()      {return m_sequence_number; }
+
+    /// Returns the number of variables.
+    inline int nb_vars() const { return size(m_indices); }
+
+    /// Returns the number of values.
+    inline int nb_values() const { return nb_vars()*nb_vars(); }
+
+    /// Reference to the matrix of values.
+    inline matrix_type& get_values() { return m_values; }
+
+    /// Constant reference to the matrix of values.
+    inline const matrix_type& get_values() const { return m_values; }
+
+    /// Reference to the indices.
+    inline const index_type& get_indices() const { return m_indices; }
+
+    /// Mutable reference to the indices.
+    inline index_type& get_indices() {	return m_indices; }
+
+    /// The actual number of non-zero values.
+    int nnz() const
+    {
+	const value_type    zero= math::zero(value_type());
+	int nbr_nz = 0;
+	for (int r = 0; r < nb_vars(); ++r)
+	    for (int c = 0; c < nb_vars(); ++c)
+		nbr_nz +=  get_values()(r,c) != zero;
+	return nbr_nz;
+    }
+
+    /// Reference to the set of neighbors.
+    neighbor_collection_type& get_neighbors() { return m_neighbors;	}
+
+    /// Reference to the set of neighbors.
+    const neighbor_collection_type& get_neighbors() const { return m_neighbors; }
+
+    /// Number of neighbors this element is connected to.
+    int get_nb_neighbors() const {	return int(m_neighbors.size()); }
+
+    /// Add new neighbors, max 6 at the time
+    void add_neighbors(element* n1, element* n2= 0, element* n3= 0,
+			element* n4= 0, element* n5= 0, element* n6= 0)
+    {
+	m_neighbors.push_back(n1);
+	if (n2) {
+	    m_neighbors.push_back(n2);
+	    if (n3) m_neighbors.push_back(n3);
+	    if (n4) m_neighbors.push_back(n4);
+	    if (n5) m_neighbors.push_back(n5);
+	    if (n6) m_neighbors.push_back(n6);
+	}
+    }
+
+
+/*******************************************************************************
+ * Useful Inspector Methods
+ ******************************************************************************/
+
+    /// The set of nodes that is incident to the element.
+    boost::unordered_set<int> get_incident_nodes() const
+    {
+	boost::unordered_set<int> nodes(2 * get_nb_neighbors());
+	for(typename neighbor_collection_type::const_iterator neigh_it = m_neighbors.begin();
+	    neigh_it != m_neighbors.end(); ++neigh_it) {
+	    element_type& neigh = **neigh_it;
+	    nodes.insert(neigh.get_indices().begin(), neigh.get_indices().end());
+	}
+	// Remove the nodes of the element.
+	for( int i = 0; i < nb_vars(); ++i )
+	    nodes.erase( get_indices()(i) );
+	return nodes;
+    }
+
+
+    /// Get the set of level-k neighbors, for a given k.
+    neighbor_set_type get_level_neighbors(const int level = 1)
+    {
+	neighbor_set_type result( get_nb_neighbors() * level );
+
+	if (level > 0) {
+	    result.insert( m_neighbors.begin(), m_neighbors.end() );
+	    if (level > 1) {
+		for(int i = 0; i < get_nb_neighbors(); ++i) {
+		    neighbor_set_type neighs(m_neighbors[i]->get_level_neighbors(level-1));
+		    result.insert( neighs.begin(), neighs.end() );
+		}
+		result.erase( this );
+	    }
+	}
+	return result;
+    }
+
+/*******************************************************************************
+ * Manipulation
+ ******************************************************************************/
+  public:
+    /// Permutes the rows and the columns of the element coefficient matrix along
+    /// with the indices such that the latter are sorted in ascending order.
+    void sort_indices() {
+	if(size(m_indices) == 0) {
+	    assert(size(m_values) == 0);
+	    return;
+	}
+
+	bool sorted = true;
+	for(int i = 0; i < nb_vars()-1; ++i) {
+	    sorted &= (get_indices()(i) < get_indices()(i+1));
+	}
+	if(sorted) {
+	    return;
+	}
+
+	index_type  orig_index( get_indices() );
+	matrix_type orig_matrix( get_values() );
+
+	std::sort(
+		  &(get_indices()(0)),
+		  &(get_indices()(0))+nb_vars()
+		  );
+
+	index_type orig_offset( nb_vars() );
+	orig_offset = -1;
+	for(int i = 0; i < nb_vars(); ++i) {
+	    int seek_idx = get_indices()(i);
+	    int j = 0;
+	    for(; (j < nb_vars()) && (orig_index(j) != seek_idx); ++j){};
+	    orig_offset(i) = j;
+	}
+
+	matrix_type& values = get_values();
+	for(int r = 0; r < nb_vars(); ++r) {
+	    for(int c = 0; c < nb_vars(); ++c) {
+		values(r,c) = orig_matrix( orig_offset(r), orig_offset(c) );
+	    }
+	}
+
+#ifndef NDEBUG
+	sorted = true;
+	for(int i = 0; i < nb_vars()-1; ++i) {
+	    sorted &= (get_indices()(i) < get_indices()(i+1));
+	}
+	assert(sorted);
+#endif
+    }
+
+
+  public:
+    /// Removes the given set of nodes from the element
+    template< class Vector >
+    void remove_nodes(const Vector& nodes, element_type& el) {
+	if(size(m_indices) == 0) {
+	    assert(size(m_values) == 0);
+	    return;
+	}
+	if(nb_vars() == 0) {
+	    return;
+	}
+
+#ifndef NDEBUG
+	bool sorted = true;
+	for(unsigned int i = 1; i < mtl::size(nodes); ++i) {
+	    sorted &= ( nodes[i-1] < nodes[i] );
+	}
+	assert(sorted);
+#endif
+
+	const int nb_nodes = mtl::size(nodes);
+
+	// Count number of remaining variables.
+	int new_nb_nodes = nb_vars();
+	{
+	    int i = 0, j = 0;
+	    while( i < nb_vars() && j < nb_nodes ) {
+		const int diff = get_indices()(i) - nodes[j];
+		if( diff < 0 ) {
+		    ++i;
+		} else if( diff > 0 ) {
+		    ++j;
+		} else {
+		    --new_nb_nodes;
+		    ++i;
+		    ++j;
+		}
+	    }
+	}
+	assert(new_nb_nodes >= 0);
+
+	// Construct new index array.
+	index_type index;
+	index_type local_index(new_nb_nodes);
+	if(new_nb_nodes > 0) {
+	    index.change_dim(new_nb_nodes);
+//	    index = new index_type(new_nb_nodes);
+	    int i = 0, j = 0, pos = 0;
+	    while( i < nb_vars() && j < nb_nodes ) {
+		const int diff = get_indices()(i) - nodes[j];
+		if( diff < 0 ) {
+		    assert( pos < new_nb_nodes );
+	//	    (*index)(pos) = get_indices()(i);
+		    index[pos] = get_indices()(i);
+		    local_index(pos) = i;
+		    ++pos;
+		    ++i;
+		} else if( diff > 0 ) {
+		    ++j;
+		} else {
+		    ++i;
+		    ++j;
+		}
+	    }
+	    while( i < nb_vars() ) {
+		assert( pos < new_nb_nodes );
+//		(*index)(pos) = get_indices()(i);
+		index[pos] = get_indices()(i);
+		local_index(pos) = i;
+		++pos;
+		++i;
+	    }
+	} else {
+//	    index = new index_type(0);
+	}
+
+	matrix_type values;
+	if(new_nb_nodes > 0) {
+	    values.change_dim( new_nb_nodes, new_nb_nodes );
+	    matrix_type tmp(get_values()), tmp2(new_nb_nodes, new_nb_nodes);
+	    for(unsigned int i=0;i<size(local_index);i++){
+		for(unsigned int j=0;j<size(local_index);j++){
+		    tmp2[i][j]=tmp[local_index(i)][local_index(j)];
+		}
+	    }
+	    values = tmp2;
+	} else {
+//	std::cout<< "ELSE\n";
+//	    values = new matrix_type(0,0);
+	}
+	// Update the neighborhood.
+	std::set<int, std::less<int>, std::allocator<int> > remove_neighs; // trouble with Clang 4.2 w/o allocator
+	for(
+	    neighbor_iterator neigh_it = m_neighbors.begin();
+	    neigh_it != m_neighbors.end();
+	    ++neigh_it
+	    ) {
+	    element_type& neigh = **neigh_it;
+
+	    // Search a matching index.
+	    bool connected = false;
+	    {
+		int i = 0, j = 0;
+		while(
+		      (i < new_nb_nodes) &&
+		      (j < neigh.nb_vars()) &&
+		      !connected
+		      ) {
+//		    const int diff = (*index)(i) - neigh.get_indices()(j);
+		    const int diff = index[i] - neigh.get_indices()(j);
+		    if(diff < 0) {
+			++i;
+		    } else if(diff > 0) {
+			++j;
+		    } else {
+			connected = true;
+		    }
+		}
+	    }
+
+	    // If not found, then remove ourself from the neighbors and vice versa.
+	    if(!connected) {
+		neighbor_iterator pos = std::find(neigh.get_neighbors().begin(), neigh.get_neighbors().end(), this );
+		if( (pos != neigh.get_neighbors().end()) && (&neigh != &el) ) {
+		    neigh.get_neighbors().erase(pos);
+		}
+		remove_neighs.insert( neigh.get_id() );
+	    }
+	}
+
+	// Remove the neighbors we're no longer connected to.
+	for(std::set<int>::iterator it = remove_neighs.begin(); it != remove_neighs.end(); ++it) {
+	    const int seek_seq_nbr = *it;
+	    for (std::size_t j = 0; j < m_neighbors.size(); ++j)
+		if (m_neighbors[j] != 0 && m_neighbors[j]->get_id() == seek_seq_nbr) {
+		    m_neighbors.erase( m_neighbors.begin()+j );
+		    break;
+		}
+	}
+
+	if(new_nb_nodes == 0) {
+	    m_neighbors.clear();
+	}
+
+	m_indices.change_dim(0);
+	m_values.change_dim(0, 0);
+	m_indices = index;
+	m_values = values;
+    }
+
+    /// Absorbs the values of the given matrix with the given index.
+    template< class Matrix, class Vector >
+    void absorb(Matrix& other_values, Vector& other_indices)
+    {
+	const value_type    zero= math::zero(value_type());
+#ifndef NDEBUG
+	bool sorted = true;
+	for(unsigned int i = 1; i < size(other_indices); ++i) {
+	    sorted &= ( other_indices(i-1) < other_indices(i) );
+	}
+	assert(sorted);
+#endif
+
+	const int other_idx_size = size( other_indices );
+
+	// Determine set of common indices.
+	const int max_common_idx =
+	    (nb_vars() < other_idx_size) ? nb_vars() : other_idx_size;
+	mtl::vector::dense_vector<int> my_idx( max_common_idx );
+	mtl::vector::dense_vector<int> ot_idx( max_common_idx );
+	int offset = 0;
+	for(int i = 0, j = 0; i < nb_vars() && j < other_idx_size; ) {
+	    int diff = (get_indices()(i) - other_indices(j));
+	    if(diff == 0) {
+		my_idx(offset) = i;
+		ot_idx(offset) = j;
+		++offset;
+		++i;
+		++j;
+	    } else if(diff < 0) {
+		++i;
+	    } else {
+		++j;
+	    }
+	}
+
+	// Absorb the values.
+	for(int i = 0; i < offset; ++i) {
+	    for(int j = i; j < offset; ++j) {
+		get_values()( my_idx(i), my_idx(j) ) += other_values( ot_idx(i), ot_idx(j) );
+		other_values( ot_idx(i), ot_idx(j) ) = zero;
+		get_values()( my_idx(j), my_idx(i) ) += other_values( ot_idx(j), ot_idx(i) );
+		other_values( ot_idx(j), ot_idx(i) ) = zero;
+	    }
+	}
+    }
+
+    /// Removes the numerical values from the element.
+    void clear() {
+	m_neighbors.clear();
+	m_neighbors.resize(1);
+	matrix_type empty;
+	swap(m_values, empty);
+	m_indices.change_dim(0);
+    }
+
+    template< class ValueType > friend class element_structure;
+
+  private:
+    /// The set of neighbors of the element.
+    neighbor_collection_type m_neighbors;
+
+    /// The set of indices of this element.
+    index_type m_indices;
+
+    /// The [Size x Size] element matrix.
+    matrix_type m_values;
+
+    /// A unique sequence number for the element, indicating it's order relative to other elements.
+    int m_sequence_number;
+
+    int *dummy;
+};
+
+
+/// Print an element to an output stream.
+template<typename OStream, class ValueType>
+OStream& operator<<(OStream& out, element<ValueType>& el)
+{
+    out << "ID: " << el.get_id() << "\n";
+    if(el.nb_vars() > 0) {
+	out << "Indices: (" << el.get_indices()(0);
+	for(int i = 1; i < el.nb_vars(); ++i) {
+	    out << ", " << el.get_indices()(i);
+	}
+	out << ")\n";
+    } else {
+	out << "Indices: ()\n";
+    }
+    out << "Neighbors: (";
+    if (el.nb_vars() > 0)
+	for(int i = 0; i < el.get_nb_neighbors(); ++i)
+	    out << el.get_neighbors()[i]->get_id() << (i+1 < el.get_nb_neighbors()? ", " : ")\n");
+    out << "Values: \n" << el.get_values();
+    return out;
+}
+
+  }} // mtl::matrix
+
+#endif // MTL_ELEMENT_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/element_array.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_ELEMENT_ARRAY_INCLUDE
+#define MTL_ELEMENT_ARRAY_INCLUDE
+
+namespace mtl { namespace matrix {
+
+template <typename Array, typename Rows, typename Cols>
+struct element_array_t
+{
+    explicit element_array_t(const Array& array, const Rows& rows, const Cols& cols)
+	: array(array), rows(rows), cols(cols)
+    {}
+
+    const Array& array;
+    const Rows&  rows;
+    const Cols&  cols;
+};
+
+
+template <typename Array, typename Rows, typename Cols>
+element_array_t<Array, Rows, Cols>
+inline element_array(const Array& array, const Rows& rows, const Cols& cols)
+{
+    return element_array_t<Array, Rows, Cols>(array, rows, cols);
+}
+
+template <typename Array, typename Rows>
+element_array_t<Array, Rows, Rows>
+inline element_array(const Array& array, const Rows& rows)
+{
+    return element_array_t<Array, Rows, Rows>(array, rows, rows);
+}
+
+} // namespace matrix
+
+using matrix::element_array;
+
+} // namespace mtl
+
+#endif // MTL_ELEMENT_ARRAY_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/element_matrix.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_ELEMENT_MATRIX_INCLUDE
+#define MTL_ELEMENT_MATRIX_INCLUDE
+
+namespace mtl { namespace matrix {
+
+
+template <typename Matrix, typename Rows, typename Cols>
+struct element_matrix_t
+{
+    explicit element_matrix_t(const Matrix& matrix, const Rows& rows, const Cols& cols)
+	: matrix(matrix), rows(rows), cols(cols)
+    {}
+
+    const Matrix&  matrix;
+    const Rows&    rows;
+    const Cols&    cols;
+};
+
+
+template <typename Matrix, typename Rows, typename Cols>
+element_matrix_t<Matrix, Rows, Cols>
+inline element_matrix(const Matrix& matrix, const Rows& rows, const Cols& cols)
+{
+    return element_matrix_t<Matrix, Rows, Cols>(matrix, rows, cols);
+}
+
+template <typename Matrix, typename Rows>
+element_matrix_t<Matrix, Rows, Rows>
+inline element_matrix(const Matrix& matrix, const Rows& rows)
+{
+    return element_matrix_t<Matrix, Rows, Rows>(matrix, rows, rows);
+}
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_ELEMENT_MATRIX_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/element_structure.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+//
+// Algorithm inspired by Nick Vannieuwenhoven, written by Cornelius Steinhardt
+
+
+
+#ifndef MTL_ELEMENT_STRUCTURE_INCLUDE
+#define MTL_ELEMENT_STRUCTURE_INCLUDE
+
+#include <iostream>
+#include <ostream>
+
+#include <boost/numeric/mtl/matrix/element.hpp>
+#include <boost/numeric/mtl/utility/ashape.hpp>
+
+namespace mtl { namespace matrix {
+
+#if 0   ///TODO need for write elmement_structur
+namespace print {
+
+	template< class Type >
+	struct print_type {
+		template<class Stream>
+		static void print(Stream&);
+	};
+
+	template< >
+	struct print_type<double> {
+		template<class Stream>
+		static void print(Stream& str) {
+			str << "double\n";
+		}
+	};
+
+	template< >
+	struct print_type<float> {
+		template<class Stream>
+		static void print(Stream& str) {
+			str << "double\n";
+		}
+	};
+
+	template< >
+	struct print_type<std::complex<double> > {
+		template<class Stream>
+		static void print(Stream& str) {
+			str << "complex\n";
+		}
+	};
+
+	template< class Type >
+	struct print_value {
+		template<class Stream>
+		static void print(Stream&);
+	};
+
+	template< >
+	struct print_value<double> {
+		template<class Stream>
+		static void print(Stream& str, double& val) {
+			str << std::scientific << std::setprecision(15);
+			str << val;
+		}
+	};
+
+	template< >
+	struct print_value<float> {
+		template<class Stream>
+		static void print(Stream& str, float& val) {
+			str << std::scientific << std::setprecision(15);
+			str << double(val);
+		}
+	};
+
+	template< >
+	struct print_value<std::complex<double> > {
+		template<class Stream>
+		static void print(Stream& str, std::complex<double>& val) {
+			str << std::scientific << std::setprecision(15);
+			str << val.real() << "\t" << val.imag() << "\t";
+		}
+	};
+}
+  #endif
+
+/**
+ * A generic abstract base class for meshes. It describes the concept of a mesh.
+ */
+template< class ValueType >
+class element_structure
+{
+
+public:
+    /// Type of the numerical values of the element coefficient matrices.
+    typedef ValueType value_type;
+
+    /// Type of the element.
+    typedef element<value_type> element_type;
+
+    /// Type of index arrays.
+    typedef typename element_type::index_type index_type;
+    // Type of the indices themselves
+    typedef typename Collection<index_type>::value_type ii_type;
+
+    /// Type of the iterator over the elements of the mesh.
+    typedef element_type* element_iterator;
+
+    /// Type of this class.
+    typedef element_structure<ValueType> this_type;
+    typedef this_type                    self;
+
+    /// Standard constructor.
+    element_structure(int total_elements= 0, int total_vars= 0, element_type* elements= 0)
+      : m_total_elements(total_elements), m_total_vars(total_vars),
+	m_elements(elements), index_heap(0), value_heap(0)
+    { }
+
+    /// consume elements into element_structure
+    void consume(int total_elements, int total_vars, element_type* elements)
+    {
+	m_total_elements= total_elements;
+	m_total_vars= total_vars;
+	delete[] m_elements;
+	m_elements= elements;
+	delete[] index_heap; index_heap= 0;
+	delete[] value_heap; value_heap= 0;
+    }
+
+
+    /// Copy the given mesh.
+    element_structure(this_type const& other)
+      : m_total_elements(other.m_total_elements),
+	m_total_vars(other.m_total_vars),
+	m_elements(m_total_elements == 0 ? 0 : new element_type[m_total_elements]),
+	index_heap(0), value_heap(0)
+    {
+	typedef typename element_type::neighbor_collection_type neigh_coll_type;
+
+	int j = 0;
+	bool ordered = true;
+	for(element_iterator it = other.element_begin(); it != other.element_end(); ++it) {
+	    // Deep copy the elements.
+	    m_elements[j] = *it;
+	    ordered &= (it->get_id() == j);
+	    ++j;
+	}
+	assert( ordered );
+	// Reconstruct the network of neighbors.
+	for(element_iterator it = this->element_begin(); it != this->element_end(); ++it) {
+	    neigh_coll_type new_neighs;
+	    neigh_coll_type& old_neighs = it->get_neighbors();
+	    for(int i = 0; i < it->get_nb_neighbors(); ++i) {
+		element_type& neigh = *(old_neighs[i]);
+		int pos = neigh.get_id();
+		new_neighs.push_back( this->m_elements+pos );
+	    }
+	    old_neighs.assign(new_neighs.begin(), new_neighs.end());
+	}
+    }
+
+    /// Destructor
+    ~element_structure() { delete[] m_elements; delete[] index_heap; delete[] value_heap; }
+
+    /// make compakt memory block from elements
+    void make_compact()
+    {
+	assert(index_heap == 0); assert(value_heap == 0); // might be relaxed later
+
+	int total_indices= 0, total_values= 0;
+	for (int i= 0; i < m_total_elements; i++) {
+	    total_indices+= m_elements[i].nb_vars();
+	    total_values+= m_elements[i].nb_values();
+	}
+	index_heap= new ii_type[total_indices];
+	value_heap= new value_type[total_values];
+
+	int index_pos= 0, value_pos= 0;
+	for (int i= 0; i < m_total_elements; i++) {
+	    element_type& element= m_elements[i];
+	    int s= element.nb_vars();
+	    index_type index_tmp(s, index_heap + index_pos);
+	    index_tmp= element.get_indices();
+	    swap(index_tmp, element.get_indices());
+	    index_pos+= s;
+
+	    typename element_type::matrix_type value_tmp(s, s, value_heap + value_pos);
+	    value_tmp= element.get_values();
+	    swap(value_tmp, element.get_values());
+	    value_pos+= s * s;
+	}
+	assert(total_indices == index_pos);
+	assert(total_values  == value_pos);
+    }
+
+    /*******************************************************************************
+     * Inspector Members
+     ******************************************************************************/
+
+    /// Total number of elements in the grid.
+    int get_total_elements() const { return m_total_elements; }
+
+    /// Total number of variables.
+    int get_total_vars() const { return m_total_vars;   }
+
+    /// Total number of non-zero values.
+    int get_total_nnz() const
+    {
+	int nnz = 0;
+	for(element_iterator it = element_begin(); it != element_end(); ++it) {
+	    nnz += it->nnz();
+	}
+	return nnz;
+    }
+
+    /// Iterator to the first element.
+    element_iterator element_begin() const { return m_elements + 0;  }
+
+    /// An iterator to the element past the last element.
+    element_iterator element_end() const { return m_elements + this->get_total_elements();   }
+
+#if 1
+    /// Writes the elements to the specified file.  TODO at the moment very slow
+    void write_to_file(const std::string& filename)
+    {
+	//using namespace print;
+
+	std::ofstream file(filename.c_str());
+
+	// Write header information.
+	file << get_total_elements() << "\n";
+	file << this->get_total_vars() << "\n";
+	//print_type<value_type>::print(file);
+
+	// Write element matrices.
+	for(element_iterator it = element_begin(); it != element_end(); ++it) {
+	    // Write indices.
+	    for(int i = 0; i < it->nb_vars()-1; ++i) {
+		file << it->get_indices()(i) << " ";
+	    }
+	    file << it->get_indices()(it->nb_vars()-1) << "\n";
+
+	    // Write values.
+	    for(int r = 0; r < it->nb_vars(); ++r) {
+		for(int c = 0; c < it->nb_vars()-1; ++c) {
+		    //print_value<value_type>::print(file, it->get_values()(r,c));
+		    file << it->get_values()(r,c); file << " ";
+		}
+		//print_value<value_type>::print(file, it->get_values()(r,it->nb_vars()-1));
+		file << it->get_values()(r,it->nb_vars()-1);
+		file << "\n";
+	    }
+	    file << "\n";
+	}
+    }
+#endif
+
+    int           m_total_elements; ///< The total number of elements.
+    int           m_total_vars; ///< The total number of variables.
+    element_type* m_elements; ///< The elements of the grid, stored consecutively.
+
+    ii_type*      index_heap;
+    value_type*   value_heap;
+};
+
+template <typename ValueType>
+inline std::size_t num_rows(const element_structure<ValueType>& A)
+{   return A.get_total_vars(); }
+
+template <typename ValueType>
+inline std::size_t num_cols(const element_structure<ValueType>& A)
+{   return A.get_total_vars(); }
+
+template <typename ValueType>
+inline std::size_t size(const element_structure<ValueType>& A)
+{   return A.get_total_vars() * A.get_total_vars(); }
+
+
+template <typename ValueType>
+inline void swap(element_structure<ValueType>& x, element_structure<ValueType>& y)
+{
+    swap(x.m_total_elements, y.m_total_elements);
+    swap(x.m_total_vars, y.m_total_vars);
+    swap(x.m_elements, y.m_elements);
+}
+
+
+}} // mtl::matrix
+
+
+#endif // MTL_ELEMENT_STRUCTURE_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/ell_matrix.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG, www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also tools/license/license.mtl.txt in the distribution.
+
+#ifndef MTL_MATRIX_ELL_MATRIX_INCLUDE
+#define MTL_MATRIX_ELL_MATRIX_INCLUDE
+
+#include <vector>
+#include <cassert>
+
+#include <boost/static_assert.hpp>
+#include <boost/numeric/mtl/matrix/parameter.hpp>
+#include <boost/numeric/mtl/matrix/compressed2D.hpp>
+#include <boost/numeric/mtl/utility/wrapped_object.hpp>
+#include <boost/numeric/mtl/utility/is_row_major.hpp>
+#include <boost/numeric/mtl/operation/std_output_operator.hpp>
+
+
+namespace mtl { namespace matrix {
+
+/// Matrix in Ell-Pack format; still in early stage, to be used with care (if at all)
+template <typename Value, typename Parameters = matrix::parameters<> >
+class ell_matrix
+  : public base_matrix<Value, Parameters>,
+    public const_crtp_base_matrix< ell_matrix<Value, Parameters>, Value, typename Parameters::size_type >,
+    public crtp_matrix_assign< ell_matrix<Value, Parameters>, Value, typename Parameters::size_type >,
+    public mat_expr< ell_matrix<Value, Parameters> >
+{
+    BOOST_STATIC_ASSERT((mtl::traits::is_row_major<Parameters>::value));
+
+    typedef base_matrix<Value, Parameters>             super;
+    typedef ell_matrix                                 self;
+    typedef mat_expr< ell_matrix<Value, Parameters> >  expr_base;
+
+    void set_stride()
+    {   my_stride= (this->dim1() + alignment - 1) / alignment * alignment; }
+
+  public:
+    typedef Parameters                                 parameters;
+    typedef typename Parameters::orientation           orientation;
+    typedef typename Parameters::dimensions            dimensions;
+    typedef Value                                      value_type;
+    typedef value_type                                 const_reference;
+
+    typedef typename Parameters::size_type             size_type;
+    typedef crtp_matrix_assign<self, Value, size_type> assign_base;
+
+    static const unsigned alignment=                   32; // TBD: make more flexible later
+
+    /// Default constructor
+    explicit ell_matrix ()
+      : super(non_fixed::dimensions(0, 0)), my_slots(0), inserting(false)
+    {  set_stride(); }
+
+    /// Construct matrix of size \p num_rows times \p num_cols
+    explicit ell_matrix (size_type num_rows, size_type num_cols)
+      : super(non_fixed::dimensions(num_rows, num_cols)), my_slots(0), inserting(false)
+    {  set_stride(); }
+
+    /// Print internal representation
+    template <typename OStream>
+    void print_internal(OStream& os) const
+    {
+#     ifdef MTL_HAS_STD_OUTPUT_OPERATOR
+	os << "indices = " << indices << '\n';
+	os << "values  = " << data << '\n';
+#     endif
+    }
+
+    /// Entry in row \p r and column \p c
+    value_type operator()(size_type r, size_type c) const
+    {
+	for (size_type k= r, i= 0; i < my_slots; ++i, k+= my_stride)
+	    if (indices[k] == c)
+		return data[k];
+	return value_type(0);
+    }
+
+    const std::vector<size_type>&  ref_minor() const { return indices; } ///< Refer index vector [advanced]
+          std::vector<size_type>&  ref_minor()       { return indices; } ///< Refer index vector [advanced]
+    const std::vector<value_type>& ref_data()  const { return data; } ///< Refer data vector [advanced]
+          std::vector<value_type>& ref_data()        { return data; } ///< Refer data vector [advanced]
+
+    size_type stride() const { return my_stride; } /// Stride [advanced]
+    size_type slots() const { return my_slots; } /// Slots, i.e. maximum number of entries per row/column
+
+    void make_empty()
+    {	my_slots= 0; indices.resize(0); data.resize(0);    }
+
+    void change_dim(size_type r, size_type c)
+    {
+	if (this->num_rows() != r || this->num_cols() != c) {
+	    super::change_dim(r, c);
+	    set_stride();
+	    make_empty();
+	}
+    }
+
+  protected:
+    void allocate_slots(size_type s)
+    {
+	my_slots= s;
+	size_type size= my_stride * s;
+	indices.resize(size); data.resize(size);
+    }
+
+    template <typename V, typename P, typename Updater> friend struct ell_matrix_inserter;
+
+    std::vector<value_type> data;
+    std::vector<size_type>  indices;
+    size_type               my_stride, my_slots;
+    bool                    inserting;
+};
+
+
+template <typename Value, typename Parameters, typename Updater = mtl::operations::update_store<Value> >
+struct ell_matrix_inserter
+  : wrapped_object<compressed2D<Value, Parameters> >,
+    compressed2D_inserter<Value, Parameters, Updater>
+{
+    typedef typename Parameters::size_type    size_type;
+    typedef Value                             value_type;
+    typedef ell_matrix<Value, Parameters>     matrix_type;
+    typedef compressed2D<Value, Parameters>   compressed_type;
+    typedef wrapped_object<compressed_type>   wrapped_type;
+    typedef compressed2D_inserter<Value, Parameters, Updater>   base_inserter;
+
+    explicit ell_matrix_inserter(matrix_type& A, size_type slot_size = 5)
+      : wrapped_type(num_rows(A), num_cols(A)),
+	base_inserter(wrapped_type::wrapped_object_member, slot_size),
+	A(A)
+    {
+	A.inserting= true;
+    }
+
+    ~ell_matrix_inserter()
+    {
+	this->finish();
+	const compressed_type& B= this->wrapped_object_member;
+	// std::cout << "Finished insertion!\nA (compressed2D) is:\n" << B;
+
+	size_type max_slots= 0;
+	for (size_type i= 0; i < B.dim1(); ++i) {
+	    size_type s= this->starts[i+1] - this->starts[i];
+	    if (s > max_slots)
+		max_slots= s;
+	}
+	A.allocate_slots(max_slots);
+
+	for (size_type i= 0; i < B.dim1(); ++i) {
+	    size_type patch_entries= max_slots - (this->starts[i+1] - this->starts[i]), k= i,
+		      patch_index= 0;
+	    for (size_type j= this->starts[i]; j < this->starts[i+1]; ++j, k+= A.my_stride) {
+		patch_index= A.indices[k]= B.ref_minor()[j];
+		A.data[k]= B.data[j];
+	    }
+	    for (size_type j= 0; j < patch_entries; ++j, k+= A.my_stride) {
+		A.indices[k]= patch_index;
+		A.data[k]= value_type(0);
+	    }
+	}
+	A.my_nnz= B.nnz();
+	A.inserting= false;
+    }
+
+    matrix_type& A;
+};
+
+// ================
+// Free functions
+// ================
+
+template <typename Value, typename Parameters>
+typename ell_matrix<Value, Parameters>::size_type
+inline num_rows(const ell_matrix<Value, Parameters>& matrix)
+{
+    return matrix.num_rows();
+}
+
+template <typename Value, typename Parameters>
+typename ell_matrix<Value, Parameters>::size_type
+inline num_cols(const ell_matrix<Value, Parameters>& matrix)
+{
+    return matrix.num_cols();
+}
+
+template <typename Value, typename Parameters>
+// typename ell_matrix<Value, Parameters>::size_type risks overflow
+std::size_t
+inline size(const ell_matrix<Value, Parameters>& matrix)
+{
+    return std::size_t(matrix.num_cols()) * std::size_t(matrix.num_rows());
+}
+
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_ELL_MATRIX_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/hermitian_view.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_HERMITIAN_VIEW_INCLUDE
+#define MTL_HERMITIAN_VIEW_INCLUDE
+
+#include <iostream>
+#include <boost/shared_ptr.hpp>
+#include <boost/numeric/mtl/matrix/map_view.hpp>
+#include <boost/numeric/mtl/matrix/transposed_view.hpp>
+#include <boost/numeric/mtl/operation/conj.hpp>
+#include <boost/numeric/mtl/operation/matrix_bracket.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+
+
+namespace mtl { namespace matrix {
+
+template <class Matrix>
+struct hermitian_view
+  : private transposed_view<Matrix>,
+    public map_view<mtl::sfunctor::conj<typename Matrix::value_type>,
+		    transposed_view<Matrix> >
+{
+    typedef transposed_view<Matrix>                                trans_base;
+    typedef mtl::sfunctor::conj<typename Matrix::value_type>       functor_type;
+    typedef map_view<functor_type, transposed_view<Matrix> >       base;
+    typedef hermitian_view                                         self;
+    typedef const Matrix&                                          const_ref_type;
+    typedef typename Collection<Matrix>::size_type                 size_type;
+    typedef typename Collection<Matrix>::value_type                value_type;
+
+    typedef typename OrientedCollection<trans_base>::orientation   orientation; // Should not be needed because defined in Collection (bug in g++???)
+
+    hermitian_view(const Matrix& matrix)
+      : trans_base(const_cast<Matrix&>(matrix)),
+	base(functor_type(), static_cast<trans_base&>(*this))
+    {}
+
+#if 0
+    hermitian_view(boost::shared_ptr<Matrix> p)
+	: trans_base(p), base(functor_type(), static_cast<trans_base&>(*this))
+    {}
+#endif
+
+    typename base::value_type operator()(size_type r, size_type c) const { return base::operator()(r, c); }
+
+    operations::bracket_proxy<self, const self&, value_type>
+    operator[] (size_type r) const
+    {
+	return operations::bracket_proxy<self, const self&, value_type>(*this, r);
+    }
+
+    friend size_type inline num_rows(const self& A) { return num_rows((const base&)(A)); }
+    friend size_type inline num_cols(const self& A) { return num_cols((const base&)(A)); }
+
+    const_ref_type const_ref() const
+    {
+	// make two statements because nvcc cannot handle ref.ref
+	const transposed_view<Matrix>& r1= base::ref;
+	return r1.ref;
+    }
+
+    size_type nnz() const { return base::nnz(); }
+
+    friend inline std::ostream& operator<<(std::ostream& os, const self& A) { return os << (const base&)(A); }
+};
+
+// If not defined ambigous between map_view and transposed_view
+template <class Matrix>
+inline std::size_t size(const hermitian_view<Matrix>& A)
+{
+    return num_rows(A) * num_rows(A);
+}
+
+// TBD submatrix of Hermitian (not trivial)
+
+
+}} // namespace mtl::matrix
+
+
+
+// Traits for Hermitian views
+namespace mtl { namespace traits {
+
+template <typename Matrix>
+struct row< mtl::matrix::hermitian_view<Matrix> >
+  : public row< mtl::matrix::map_view<sfunctor::conj<typename Matrix::value_type>,
+				      mtl::matrix::transposed_view<Matrix> > >
+{};
+
+template <typename Matrix>
+struct col< mtl::matrix::hermitian_view<Matrix> >
+    : public col< mtl::matrix::map_view<sfunctor::conj<typename Matrix::value_type>,
+				   mtl::matrix::transposed_view<Matrix> > >
+{};
+
+template <typename Matrix>
+struct const_value< mtl::matrix::hermitian_view<Matrix> >
+    : public const_value< mtl::matrix::map_view<sfunctor::conj<typename Matrix::value_type>,
+					   mtl::matrix::transposed_view<Matrix> > >
+{};
+
+template <typename Tag, typename Matrix>
+struct range_generator< Tag, mtl::matrix::hermitian_view<Matrix> >
+    : public range_generator< Tag, mtl::matrix::map_view<sfunctor::conj<typename Matrix::value_type>,
+						    mtl::matrix::transposed_view<Matrix> > >
+{};
+
+template <typename Matrix>
+struct range_generator< tag::major, mtl::matrix::hermitian_view<Matrix> >
+    : public range_generator< tag::major, mtl::matrix::map_view<sfunctor::conj<typename Matrix::value_type>,
+							   mtl::matrix::transposed_view<Matrix> > >
+{};
+
+
+}} // mtl::traits
+
+#endif // MTL_HERMITIAN_VIEW_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/hessian_setup.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_HESSIAN_MATRIX_UTILITIES_INCLUDE
+#define MTL_HESSIAN_MATRIX_UTILITIES_INCLUDE
+
+#include <cmath>
+#include <iostream>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/matrix/inserter.hpp>
+#include <boost/numeric/mtl/operation/entry_similar.hpp>
+
+namespace mtl { namespace matrix {
+
+/// Fills a matrix A with A[i][j] = factor * (i + j)
+/** Intended for dense matrices.
+    Works on sparse matrices with inserter but is very expensive. **/
+template <typename Matrix, typename Value>
+void hessian_setup(Matrix& A, Value factor)
+{
+    typedef typename Collection<Matrix>::value_type    value_type;
+    typedef typename Collection<Matrix>::size_type     size_type;
+
+    inserter<Matrix> ins(A, num_cols(A));
+
+    for (size_type r= 0; r < num_rows(A); r++)
+	for (size_type c= 0; c < num_cols(A); c++)
+	    ins[r][c] << factor * (value_type(r) + value_type(c));
+}
+
+namespace impl {
+
+    /*
+    - Check matrix product C = A * B with:
+      - A is MxN, B is NxL, C is MxL
+      - with matrices a_ij = i+j, b_ij = 2(i+j);
+      - c_ij = 1/3 N (1 - 3i - 3j + 6ij - 3N + 3iN + 3jN + 2N^2).
+
+    */
+    // Not really generic
+    template <typename Value>
+    double inline hessian_product_i_j (Value i, Value j, Value N)
+    {
+        return 1.0/3.0 * N * (1.0 - 3*i - 3*j + 6*i*j - 3*N + 3*i*N + 3*j*N + 2*N*N);
+    }
+
+    template <typename Value>
+    inline bool similar_values(Value x, Value y)
+    {
+        using std::abs; using std::max;
+        return abs(x - y) / max(abs(x), abs(y)) < 0.000001;
+    }
+
+    template <typename Matrix>
+    void inline check_entry(Matrix const& C, unsigned long r, unsigned long c,
+			    unsigned long reduced_dim, double factor)
+    {
+	if (!entry_similar(C, r, c, factor * hessian_product_i_j(r, c, reduced_dim), 0.00001)) {
+	    std::cerr << "Result in C[" << r << "][" << c << "] should be "
+		      << factor * hessian_product_i_j(r, c, reduced_dim)
+		      << " but is " << C[r][c] << "\n";
+	    MTL_THROW(unexpected_result());
+	}
+    }
+
+} // impl
+
+
+/// Check if matrix C is A * B with A and B set by hessian_setup
+/** C has dimensions M x L and reduced_dim is N, see hessian_setup. **/
+template <typename Matrix>
+void check_hessian_matrix_product(Matrix const& C, typename Matrix::size_type reduced_dim, double factor= 1.0)
+{
+    if (num_rows(C) * num_cols(C) == 0) return; // otherwise out of range
+
+    impl::check_entry(C, 0, 0, reduced_dim, factor);
+    impl::check_entry(C, 0, num_cols(C)-1, reduced_dim, factor);
+    impl::check_entry(C, num_rows(C)-1, 0, reduced_dim, factor);
+    impl::check_entry(C, num_rows(C)-1, num_cols(C)-1, reduced_dim, factor);
+    impl::check_entry(C, num_rows(C)/2, num_cols(C)/2, reduced_dim, factor);
+}
+
+} // namespace matrix;
+
+using matrix::hessian_setup;
+using matrix::check_hessian_matrix_product;
+
+} // namespace mtl
+
+#endif // MTL_HESSIAN_MATRIX_UTILITIES_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/identity.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_MATRIX_IDENTITY_INCLUDE
+#define MTL_MATRIX_IDENTITY_INCLUDE
+
+// #include <boost/numeric/linear_algebra/identity.hpp>
+// #include <boost/numeric/mtl/mtl_fwd.hpp>
+// #include <boost/numeric/mtl/matrix/parameter.hpp>
+// #include <boost/numeric/mtl/matrix/diagonal_setup.hpp>
+
+#include <boost/numeric/mtl/matrix/identity2D.hpp>
+
+namespace mtl { namespace matrix {
+
+inline identity2D identity(std::size_t nrows, std::size_t ncols)
+{
+    return identity2D(nrows, ncols);
+}
+
+
+inline identity2D identity(std::size_t nrows)
+{
+    return identity2D(nrows, nrows);
+}
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_IDENTITY_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/identity2D.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG, www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also tools/license/license.mtl.txt in the distribution.
+
+#ifndef MTL_MATRIX_IDENTITY2D_INCLUDE
+#define MTL_MATRIX_IDENTITY2D_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/ashape.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/irange.hpp>
+#include <boost/numeric/mtl/vector/mat_cvec_multiplier.hpp>
+
+namespace mtl { namespace matrix {
+
+/// Matrix-free linear operator for identity
+struct identity2D
+{
+    /// Constructor for \p m by \p m matrix
+    identity2D(std::size_t m) : m(m), n(m) {}
+
+    /// Constructor for \p m by \p n matrix
+    identity2D(std::size_t m, std::size_t n) : m(m), n(n) {}
+
+    /// Member function that realizes the multiplication
+    template <typename VectorIn, typename VectorOut, typename Assign>
+    void mult(const VectorIn& v, VectorOut& w, Assign) const
+    {
+	MTL_DEBUG_THROW_IF(std::size_t(size(v)) != n, incompatible_size());
+	MTL_DEBUG_THROW_IF(size(w) != 0 && std::size_t(size(w)) != m, incompatible_size());
+
+	if (size(w) == 0)
+	    w.change_dim(m);
+
+	if (m == n)
+	    Assign::first_update(w, v);
+	else if (m < n)
+	    Assign::first_update(w, v[irange(m)]);
+	else {
+	    VectorOut w1(w[irange(n)]), w2(w[irange(n, imax)]);
+	    Assign::first_update(w1, v);
+	    Assign::init(w2);
+	}
+    }
+
+    /// Multiplication is procastinated until we know where the product goes
+    template <typename VectorIn>
+    vector::mat_cvec_multiplier<identity2D, VectorIn> operator*(const VectorIn& v) const
+    {	return vector::mat_cvec_multiplier<identity2D, VectorIn>(*this, v);    }
+
+    std::size_t m, n;
+};
+
+inline std::size_t size(const identity2D& A) { return A.m * A.n; } ///< Matrix size
+inline std::size_t num_rows(const identity2D& A) { return A.m; } ///< Number of rows
+inline std::size_t num_cols(const identity2D& A) { return A.n; } ///< Number of columns
+
+}} // namespace mtl::matrix
+
+namespace mtl {
+
+    template <>
+    struct Collection<matrix::identity2D>
+    {
+	typedef double         value_type;
+	typedef std::size_t    size_type;
+    };
+
+    namespace ashape {
+	template <> struct ashape_aux<matrix::identity2D>
+	{	typedef nonscal type;    };
+    }
+}
+
+#endif // MTL_MATRIX_IDENTITY2D_INCLUDE