AMDiS/lib/mtl4/boost/numeric/mtl/matrix/implicit_dense.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_IMPLICIT_DENSE_INCLUDE
+#define MTL_MATRIX_IMPLICIT_DENSE_INCLUDE
+
+#include <vector>
+#include <boost/numeric/linear_algebra/inverse.hpp>
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/matrix/crtp_base_matrix.hpp>
+#include <boost/numeric/mtl/matrix/mat_expr.hpp>
+#include <boost/numeric/mtl/concept/std_concept.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+
+#ifdef MTL_HAS_MPI
+#   include <boost/numeric/mtl/matrix/distributed.hpp>
+#   include <boost/numeric/mtl/matrix/inserter.hpp>
+#   include <boost/mpi/collectives.hpp>
+#endif
+
+namespace mtl { namespace matrix {
+
+template <typename Functor>
+class implicit_dense
+  : public const_crtp_base_matrix< implicit_dense<Functor>,
+				   typename Functor::result_type, typename Functor::size_type >,
+    public mat_expr< implicit_dense<Functor> >
+{
+    typedef implicit_dense                             self;
+  public:
+    typedef mtl::tag::row_major                        orientation; // for completeness
+    typedef typename Functor::result_type              value_type;
+    typedef typename Functor::result_type              const_reference;
+    typedef typename Functor::size_type                size_type;
+
+    // Should not be needed, to be removed after transposed_view is cleaned up
+    typedef index::c_index                             index_type;
+    typedef mtl::traits::detail::matrix_element_key<self> key_type;
+    typedef mtl::non_fixed::dimensions                 dim_type;
+
+    explicit implicit_dense (const Functor& functor) : my_functor(functor) {}
+
+    value_type operator() (size_type r, size_type c) const { return my_functor(r, c); }
+
+    size_type nnz() const { return dim1() * dim2(); }
+    size_type dim1() const { return num_rows(my_functor); }
+    size_type dim2() const { return num_cols(my_functor); }
+
+    friend size_type inline num_rows(const self& A) { return num_rows(A.my_functor); }
+    friend size_type inline num_cols(const self& A) { return num_cols(A.my_functor); }
+
+    Functor& functor() { return my_functor; }
+    Functor const& functor() const { return my_functor; }
+
+  private:
+    Functor           my_functor;
+};
+
+template <typename Functor>
+typename Functor::size_type inline size(const implicit_dense<Functor>& A)
+{
+    return num_rows(A) * num_cols(A);
+}
+
+// ==========
+// Sub matrix
+// ==========
+
+// To do later
+
+}} // namespace mtl::matrix
+
+
+namespace mtl { namespace traits {
+
+    // ================
+    // Range generators
+    // For cursors
+    // ================
+
+    template <typename Functor>
+    struct range_generator<glas::tag::row, mtl::matrix::implicit_dense<Functor> >
+      : detail::all_rows_range_generator<mtl::matrix::implicit_dense<Functor>, complexity_classes::linear>
+    {};
+
+    template <typename Functor>
+    struct range_generator<glas::tag::major, mtl::matrix::implicit_dense<Functor> >
+      : range_generator<glas::tag::row, mtl::matrix::implicit_dense<Functor> >
+    {};
+
+    template <typename Functor>
+    struct range_generator<glas::tag::nz,
+			   detail::sub_matrix_cursor<mtl::matrix::implicit_dense<Functor>, glas::tag::row, 2> >
+      : detail::all_cols_in_row_range_generator<detail::sub_matrix_cursor<mtl::matrix::implicit_dense<Functor>, glas::tag::row, 2> >
+    {};
+
+    template <typename Functor>
+    struct range_generator<glas::tag::col, mtl::matrix::implicit_dense<Functor> >
+      : detail::all_cols_range_generator<mtl::matrix::implicit_dense<Functor>, complexity_classes::linear>
+    {};
+
+    template <typename Functor>
+    struct range_generator<glas::tag::nz,
+			   detail::sub_matrix_cursor<mtl::matrix::implicit_dense<Functor>, glas::tag::col, 2> >
+      : detail::all_rows_in_col_range_generator<detail::sub_matrix_cursor<mtl::matrix::implicit_dense<Functor>, glas::tag::col, 2> >
+    {};
+
+
+    template <typename Tag, typename Value>
+    struct range_generator<Tag, mtl::matrix::ones_matrix<Value> >
+      : public range_generator<Tag, mtl::matrix::implicit_dense<mtl::matrix::ones_functor<Value> > >
+    {};
+
+    template <typename Value>
+    struct range_generator<glas::tag::major, mtl::matrix::ones_matrix<Value> >
+      : public range_generator<glas::tag::major, mtl::matrix::implicit_dense<mtl::matrix::ones_functor<Value> > >
+    {};
+
+    template <typename Tag, typename Value>
+    struct range_generator<Tag, mtl::matrix::hilbert_matrix<Value> >
+      : public range_generator<Tag, mtl::matrix::implicit_dense<mtl::matrix::hilbert_functor<Value> > >
+    {};
+
+    template <typename Value>
+    struct range_generator<glas::tag::major, mtl::matrix::hilbert_matrix<Value> >
+      : public range_generator<glas::tag::major, mtl::matrix::implicit_dense<mtl::matrix::hilbert_functor<Value> > >
+    {};
+
+    template <typename Tag, typename Vector1, typename Vector2>
+    struct range_generator<Tag, mtl::matrix::outer_product_matrix<Vector1, Vector2> >
+      : public range_generator<Tag, mtl::matrix::implicit_dense<mtl::matrix::outer_product_functor<Vector1, Vector2> > >
+    {};
+
+    template <typename Vector1, typename Vector2>
+    struct range_generator<glas::tag::major, mtl::matrix::outer_product_matrix<Vector1, Vector2> >
+      : public range_generator<glas::tag::major, mtl::matrix::implicit_dense<mtl::matrix::outer_product_functor<Vector1, Vector2> > >
+    {};
+
+}} // mtl::traits
+
+// =============
+// Some functors
+// =============
+
+
+namespace mtl { namespace matrix {
+
+template <typename Value= int>
+class ones_functor
+{
+    typedef ones_functor    self;
+  public:
+    typedef std::size_t    size_type;
+    typedef Value          result_type;
+
+    ones_functor(size_type nr, size_type nc) : nr(nr), nc(nc) {}
+
+    friend size_type inline num_rows(const self& A) { return A.nr; }
+    friend size_type inline num_cols(const self& A) { return A.nc; }
+
+    result_type operator()(size_type, size_type) const { return Value(1); }
+
+  private:
+    size_type nr, nc;
+};
+
+
+template <typename Value= double>
+class hilbert_functor
+{
+    typedef hilbert_functor    self;
+  public:
+    typedef std::size_t    size_type;
+    typedef Value          result_type;
+
+    hilbert_functor(size_type nr, size_type nc) : nr(nr), nc(nc) {}
+
+    friend size_type inline num_rows(const self& A) { return A.nr; }
+    friend size_type inline num_cols(const self& A) { return A.nc; }
+
+    result_type operator()(size_type r, size_type c) const
+    {
+	using math::reciprocal;
+	return reciprocal(Value(r + c + 1));
+    }
+  private:
+    size_type nr, nc;
+};
+
+
+template <typename Vector1, typename Vector2>
+class outer_product_functor
+{
+    typedef outer_product_functor    self;
+  public:
+    typedef std::size_t              size_type;
+    typedef typename Multiplicable<typename Collection<Vector1>::value_type,
+				   typename Collection<Vector2>::value_type>::result_type   result_type;
+
+    outer_product_functor(const Vector1& v1, const Vector2& v2) : my_v1(v1), my_v2(v2) {}
+    outer_product_functor(size_type r, size_type c) : my_v1(r), my_v2(c) {}
+
+    friend size_type inline num_rows(const self& A) { return size(A.my_v1); }
+    friend size_type inline num_cols(const self& A) { return size(A.my_v2); }
+
+    result_type operator()(size_type r, size_type c) const { return my_v1[r] * my_v2[c]; }
+
+    Vector1&       v1()       { return my_v1; }
+    Vector1 const& v1() const { return my_v1; }
+    Vector2&       v2()       { return my_v2; }
+    Vector2 const& v2() const { return my_v2; }
+
+  private:
+    Vector1  my_v1; // keeps copy
+    Vector2  my_v2;
+};
+
+// ======================
+// Some implicit matrices
+// ======================
+
+
+template <typename Value= int>
+class ones_matrix
+  : public implicit_dense<ones_functor<Value> >
+{
+  public:
+    typedef ones_functor<Value>                functor_type;
+    typedef typename functor_type::size_type  size_type;
+    typedef implicit_dense<functor_type>      base;
+
+    ones_matrix(size_type r, size_type c) : base(functor_type(r, c)) {}
+};
+
+
+template <typename Value= double>
+class hilbert_matrix
+  : public implicit_dense<hilbert_functor<Value> >
+{
+  public:
+    typedef hilbert_functor<Value>            functor_type;
+    typedef typename functor_type::size_type  size_type;
+    typedef implicit_dense<functor_type>      base;
+
+    hilbert_matrix(size_type r, size_type c) : base(functor_type(r, c)) {}
+};
+
+
+template <typename Vector1, typename Vector2>
+class outer_product_matrix
+  : public implicit_dense<outer_product_functor<Vector1, Vector2> >
+{
+    typedef outer_product_matrix self;
+  public:
+    typedef outer_product_functor<Vector1, Vector2>  functor_type;
+    typedef typename functor_type::size_type  size_type;
+    typedef implicit_dense<functor_type>      base;
+
+    outer_product_matrix(const Vector1& v1, const Vector2& v2) : base(functor_type(v1, v2)) {}
+    outer_product_matrix(size_type r, size_type c) : base(functor_type(r, c)) {}
+
+    Vector1&       v1()       { return this->functor().v1(); }
+    Vector1 const& v1() const { return this->functor().v1(); }
+    Vector2&       v2()       { return this->functor().v2(); }
+    Vector2 const& v2() const { return this->functor().v2(); }
+};
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_IMPLICIT_DENSE_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/indirect.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_INDIRECT_INCLUDE
+#define MTL_MATRIX_INDIRECT_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/detail/range_generator.hpp>
+#include <boost/numeric/mtl/utility/complexity.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/iset.hpp>
+#include <boost/numeric/mtl/operation/is_negative.hpp>
+#include <boost/numeric/mtl/matrix/mat_expr.hpp>
+
+namespace mtl { namespace matrix {
+
+/// Class for indirect access to matrix
+/** So far only constant access, mutable access for appropriate matrices might be added later. **/
+template <typename Matrix>
+struct indirect
+  : mat_expr<indirect<Matrix> >,
+    const_crtp_base_matrix<indirect<Matrix>, typename Collection<Matrix>::value_type, typename Collection<Matrix>::size_type>
+{
+    typedef indirect                                self;
+    typedef Matrix                                  other;
+    typedef typename Collection<Matrix>::value_type value_type;
+    typedef typename Collection<Matrix>::size_type  size_type;
+    // if implementation uses const_reference -> change collection accordingly
+
+
+    /// Construct from constant matrix reference and isets for rows and columns
+    indirect(const Matrix& ref, const iset& rows, const iset& cols) : ref(ref), rows(rows), cols(cols) {}
+
+    friend size_type inline num_rows(const self& A) { return A.rows.size(); } ///< Number of rows
+    friend size_type inline num_cols(const self& A) { return A.cols.size(); } ///< Number of colums
+    size_type nnz() const { return num_rows(*this) * num_cols(*this); } ///< Number of non-zeros
+
+    size_type dim1() const { return rows.size(); } ///< Dimension 1 is equal to number of rows
+    size_type dim2() const { return cols.size(); } ///< Dimension 2 is equal to number of columns
+
+    value_type operator() (size_type r, size_type c) const
+    {   return (ref(rows[r], cols[c]));    } ///< Read A[r][c]
+
+  private:
+    const Matrix& ref;
+    iset          rows, cols;
+};
+
+template <typename Matrix>
+inline std::size_t size(const indirect<Matrix>& A)
+{     return num_rows(A) * num_rows(A); }
+
+}} // namespace mtl::matrix
+
+
+// -- Range generators in utility/range_generator.hpp
+// -- Property maps in utility/property_maps.hpp
+
+
+#endif // MTL_MATRIX_INDIRECT_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/inserter.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_INSERTER_INCLUDE
+#define MTL_MATRIX_INSERTER_INCLUDE
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/operation/update.hpp>
+#include <boost/numeric/mtl/detail/trivial_inserter.hpp>
+
+
+
+namespace mtl { namespace matrix {
+
+/// Matrix inserter
+/** The matrix inserter has two template arguments: the type of the target matrix and
+    an update functor.
+    The update functor determines how an existing entry is updated: overwritten, added,
+    subtracted...
+    The default is to overwrite existing entries.
+**/
+template <typename Matrix,
+	  typename Updater = mtl::operations::update_store<typename Matrix::value_type> >
+struct inserter
+  : public mtl::detail::trivial_inserter<Matrix, Updater>
+{
+    typedef mtl::detail::trivial_inserter<Matrix, Updater>     base;
+    typedef typename Matrix::size_type                         size_type;
+
+    explicit inserter(Matrix& matrix, size_type slot_size = 0) : base(matrix, slot_size) {}
+};
+
+
+template <typename Value, typename Parameters, typename Updater>
+struct inserter<compressed2D<Value, Parameters>, Updater>
+  : compressed2D_inserter<Value, Parameters, Updater>
+{
+    typedef compressed2D<Value, Parameters>                    matrix_type;
+    typedef typename matrix_type::size_type                    size_type;
+    typedef compressed2D_inserter<Value, Parameters, Updater > base;
+
+    explicit inserter(matrix_type& matrix, size_type slot_size = 5) : base(matrix, slot_size) {}
+};
+
+template <typename Value, typename Parameters, typename Updater>
+struct inserter<coordinate2D<Value, Parameters>, Updater>
+  : coordinate2D_inserter<coordinate2D<Value, Parameters>, Updater>
+{
+    typedef coordinate2D<Value, Parameters>                                  matrix_type;
+    typedef typename Parameters::size_type                                   size_type;
+    typedef coordinate2D_inserter<coordinate2D<Value, Parameters>, Updater>  base;
+
+    explicit inserter(matrix_type& matrix, size_type slot_size= 1) : base(matrix, slot_size) {}
+};
+
+template <typename Value, typename Parameters, typename Updater>
+struct inserter<sparse_banded<Value, Parameters>, Updater>
+  : sparse_banded_inserter<Value, Parameters, Updater>
+{
+    typedef sparse_banded<Value, Parameters>                                 matrix_type;
+    typedef typename Parameters::size_type                                   size_type;
+    typedef sparse_banded_inserter<Value, Parameters, Updater>  base;
+
+    explicit inserter(matrix_type& matrix, size_type slot_size= 1) : base(matrix, slot_size) {}
+};
+
+template <typename Value, typename Parameters, typename Updater>
+struct inserter<ell_matrix<Value, Parameters>, Updater>
+  : ell_matrix_inserter<Value, Parameters, Updater>
+{
+    typedef ell_matrix<Value, Parameters>                    matrix_type;
+    typedef typename matrix_type::size_type                    size_type;
+    typedef ell_matrix_inserter<Value, Parameters, Updater > base;
+
+    explicit inserter(matrix_type& matrix, size_type slot_size = 5) : base(matrix, slot_size) {}
+};
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_INSERTER_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/laplacian_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_LAPLACIAN_SETUP_INCLUDE
+#define MTL_LAPLACIAN_SETUP_INCLUDE
+
+#include <boost/numeric/mtl/matrix/inserter.hpp>
+#include <boost/numeric/mtl/operation/set_to_zero.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace mtl { namespace matrix {
+
+/// Setup a matrix according to a Laplacian equation on a 2D-grid using a five-point-stencil
+/** Intended for sparse matrices but works also with dense matrices. Changes the size of
+    the matrix \f$m\cdot n\times m\cdot n\f$. **/
+template <typename Matrix>
+inline void laplacian_setup(Matrix& A, unsigned m, unsigned n)
+{
+    vampir_trace<3063> tracer;
+    A.change_dim(m*n, m*n);
+    set_to_zero(A);
+    inserter<Matrix>      ins(A, 5);
+
+    for (unsigned i= 0; i < m; i++)
+	for (unsigned j= 0; j < n; j++) {
+	    typename Collection<Matrix>::value_type four(4.0), minus_one(-1.0);
+	    unsigned row= i * n + j;
+	    ins(row, row) << four;
+	    if (j < n-1) ins(row, row+1) << minus_one;
+	    if (i < m-1) ins(row, row+n) << minus_one;
+	    if (j > 0) ins(row, row-1) << minus_one;
+	    if (i > 0) ins(row, row-n) << minus_one;
+	}
+}
+
+}} // namespace mtl::matrix
+
+#endif // MTL_LAPLACIAN_SETUP_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/lower.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_LOWER_INCLUDE
+#define MTL_MATRIX_LOWER_INCLUDE
+
+namespace mtl { namespace matrix {
+
+namespace traits {
+
+    template <typename Matrix>
+    struct lower
+    {
+	typedef typename traits::bands<Matrix>::type type;
+    };
+}
+
+/// Lower triangular matrix
+template <typename Matrix>
+typename traits::lower<Matrix>::type
+inline lower(const Matrix& A)
+{
+    return bands(A, std::numeric_limits<long>::min(), 1);
+}
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_LOWER_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/make_fast_multi_vector_expr.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_MAKE_FAST_MULTI_VECTOR_EXPR_INCLUDE
+#define MTL_MATRIX_MAKE_FAST_MULTI_VECTOR_EXPR_INCLUDE
+
+#include <boost/numeric/mtl/utility/fast_multi_vector_expr.hpp>
+#include <boost/numeric/mtl/matrix/mat_mat_plus_expr.hpp>
+#include <boost/numeric/mtl/matrix/mat_mat_minus_expr.hpp>
+#include <boost/numeric/mtl/matrix/map_view.hpp>
+#include <boost/numeric/mtl/vector/vec_vec_plus_expr.hpp>
+#include <boost/numeric/mtl/vector/vec_vec_minus_expr.hpp>
+#include <boost/numeric/mtl/vector/map_view.hpp>
+
+namespace mtl { namespace matrix {
+
+
+template <typename E1, typename E2>
+typename mtl::traits::fast_multi_vector_expr< mv_mv_plus_expr<E1, E2> >::type
+inline make_fast_multi_vector_expr(const mv_mv_plus_expr<E1, E2>& expr)
+{
+    typedef typename mtl::traits::fast_multi_vector_expr< mv_mv_plus_expr<E1, E2> >::type type;
+    return type(make_fast_multi_vector_expr(expr.first), make_fast_multi_vector_expr(expr.second));
+}
+
+template <typename E1, typename E2>
+typename mtl::traits::fast_multi_vector_expr< mv_mv_minus_expr<E1, E2> >::type
+inline make_fast_multi_vector_expr(const mv_mv_minus_expr<E1, E2>& expr)
+{
+    typedef typename mtl::traits::fast_multi_vector_expr< mv_mv_minus_expr<E1, E2> >::type type;
+    return type(make_fast_multi_vector_expr(expr.first), make_fast_multi_vector_expr(expr.second));
+}
+
+template <typename Functor, typename Matrix>
+typename mtl::traits::fast_multi_vector_expr< map_view<Functor, Matrix> >::type
+inline make_fast_multi_vector_expr(const map_view<Functor, Matrix>& expr)
+{
+    typedef typename mtl::traits::fast_multi_vector_expr< map_view<Functor, Matrix> >::type type;
+    return type(expr.functor, make_fast_multi_vector_expr(expr.ref));
+}
+
+template <typename Value1, typename Matrix>
+typename mtl::traits::fast_multi_vector_expr< scaled_view<Value1, Matrix> >::type
+inline make_fast_multi_vector_expr(const scaled_view<Value1, Matrix>& expr)
+{
+    typedef typename mtl::traits::fast_multi_vector_expr< scaled_view<Value1, Matrix> >::type type;
+    return type(expr.functor.value(), make_fast_multi_vector_expr(expr.ref));
+}
+
+template <typename Value1, typename Matrix>
+typename mtl::traits::fast_multi_vector_expr< rscaled_view<Value1, Matrix> >::type
+inline make_fast_multi_vector_expr(const rscaled_view<Value1, Matrix>& expr)
+{
+    typedef typename mtl::traits::fast_multi_vector_expr< rscaled_view<Value1, Matrix> >::type type;
+    return type(make_fast_multi_vector_expr(expr.ref), expr.functor.value());
+}
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_MAKE_FAST_MULTI_VECTOR_EXPR_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/map_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_MAP_VIEW_INCLUDE
+#define MTL_MAP_VIEW_INCLUDE
+
+#include <utility>
+#include <boost/utility/enable_if.hpp>
+#include <boost/mpl/if.hpp>
+#include <boost/type_traits/is_integral.hpp>
+#include <boost/shared_ptr.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/is_multi_vector_expr.hpp>
+#include <boost/numeric/mtl/matrix/crtp_base_matrix.hpp>
+#include <boost/numeric/mtl/operation/sub_matrix.hpp>
+#include <boost/numeric/mtl/operation/sfunctor.hpp>
+#include <boost/numeric/mtl/operation/tfunctor.hpp>
+#include <boost/numeric/mtl/operation/conj.hpp>
+#include <boost/numeric/mtl/operation/imag.hpp>
+#include <boost/numeric/mtl/operation/real.hpp>
+#include <boost/numeric/mtl/matrix/mat_expr.hpp>
+#include <boost/numeric/mtl/vector/map_view.hpp>
+
+
+namespace mtl { namespace matrix { namespace detail {
+    // Forward declaration for friend declaration
+    template <typename, typename> struct map_value;
+
+    template <typename Functor, typename Matrix> struct map_vector {}; // not defined in general
+    template <typename Functor, typename Vector>
+    struct map_vector<Functor, mtl::matrix::multi_vector<Vector> >
+    {
+	typedef mtl::vector::map_view<Functor, Vector> type;
+    };
+
+}}}
+
+namespace mtl { namespace matrix {
+
+template <typename Functor, typename Matrix>
+struct map_view
+  : public const_crtp_base_matrix< map_view<Functor, Matrix>,
+				   typename Functor::result_type, typename Matrix::size_type >,
+    public mat_expr< map_view<Functor, Matrix> >
+{
+    typedef map_view                                   self;
+    typedef mat_expr< self >                           expr_base;
+    typedef Matrix                                     other;
+    typedef const Matrix&                              const_ref_type;
+    typedef typename Matrix::orientation               orientation;
+
+    typedef typename Functor::result_type              value_type;
+    typedef typename Functor::result_type              const_reference;
+
+    typedef typename Matrix::key_type                  key_type;
+    typedef typename Matrix::size_type                 size_type;
+    typedef typename Matrix::index_type                index_type;
+    typedef typename Matrix::dim_type                  dim_type;
+
+    struct dummy { typedef void type; };
+    typedef typename boost::mpl::eval_if<mtl::traits::is_multi_vector_expr<Matrix>, detail::map_vector<Functor, Matrix>, dummy>::type vector_type;
+
+    map_view (const Functor& functor, const other& ref) : functor(functor), ref(ref) {}
+
+    map_view (const Functor& functor, boost::shared_ptr<Matrix> p)
+      : functor(functor), my_copy(p), ref(*p) {}
+
+#ifdef MTL_WITH_CPP11_MOVE
+  map_view (self&& that) : functor(that.functor), my_copy(std::move(that.my_copy)), ref(that.ref) {}
+  map_view (const self& that) : functor(that.functor), ref(that.ref) { assert(that.my_copy.use_count() == 0); }
+#endif
+
+
+    value_type operator() (size_type r, size_type c) const
+    {
+        return functor(ref(r, c));
+    }
+    // for multi_vector, needs enable_if since only defined for multi_vector
+    template <typename S>
+    typename boost::lazy_enable_if<boost::is_integral<S>, detail::map_vector<Functor, Matrix> >::type
+    vector(S c) const
+    {
+	return typename detail::map_vector<Functor, Matrix>::type(functor, ref.vector(c));
+    }
+
+    size_type dim1() const { return ref.dim1(); }
+    size_type dim2() const { return ref.dim2(); }
+    dim_type dimensions() const { return ref.dimensions(); }
+
+    size_type begin_row() const { return ref.begin_row(); }
+    size_type end_row() const { return ref.end_row(); }
+    size_type begin_col() const { return ref.begin_col(); }
+    size_type end_col() const {	return ref.end_col(); }
+
+    size_type nnz() const { return ref.nnz(); }
+
+    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);     }
+    template <typename, typename> friend struct detail::map_value;
+
+  protected:
+    boost::shared_ptr<Matrix>           my_copy;
+  public:
+    Functor           functor;
+    const other&      ref;
+};
+
+template <typename Functor, typename Matrix>
+inline std::size_t size(const map_view<Functor, Matrix>& A)
+{     return num_rows(A) * num_rows(A); }
+
+// ==========
+// Sub matrix
+// ==========
+
+template <typename Functor, typename Matrix>
+struct sub_matrix_t< mtl::matrix::map_view<Functor, Matrix> >
+{
+    typedef mtl::matrix::map_view<Functor, Matrix>                                     view_type;
+
+    // Mapping of sub-matrix type
+    typedef typename sub_matrix_t<Matrix>::const_sub_matrix_type                       ref_sub_type;
+    typedef mtl::matrix::map_view<Functor, ref_sub_type>                               const_sub_matrix_type;
+    typedef typename view_type::size_type                                              size_type;
+
+    const_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 map_view
+	// Functor is copied from view
+	pointer_type p(new ref_sub_type(sub_matrix(view.ref, begin_r, end_r, begin_c, end_c)));
+	return const_sub_matrix_type(view.functor, p);
+    }
+};
+
+
+}} // namespace mtl::matrix
+
+
+namespace mtl { namespace traits {
+
+    namespace detail {
+
+
+	template <typename Functor, typename Matrix>
+	struct map_value
+	{
+	    typedef typename Matrix::key_type                      key_type;
+	    typedef typename mtl::matrix::map_view<Functor, Matrix>::value_type value_type;
+
+	    map_value(mtl::matrix::map_view<Functor, Matrix> const& map_matrix)
+	      : map_matrix(map_matrix), its_value(map_matrix.ref)
+	    {}
+
+	    value_type operator() (key_type const& key) const
+	    {
+		return map_matrix.functor(its_value(key));
+	    }
+
+	  protected:
+	    mtl::matrix::map_view<Functor, Matrix> const&        map_matrix;
+		typename mtl::traits::const_value<Matrix>::type its_value;
+        };
+
+
+
+	template <typename Functor, typename Matrix>
+	struct mapped_row
+	{
+	    typedef typename Matrix::key_type   key_type;
+	    typedef typename Matrix::size_type  size_type;
+
+		explicit mapped_row(const mtl::matrix::map_view<Functor, Matrix>& view) : its_row(view.ref) {}
+		explicit mapped_row(const mtl::matrix::banded_view<Matrix>& view) : its_row(view.ref) {}
+
+	    size_type operator() (key_type const& key) const
+	    {
+		return its_row(key);
+	    }
+
+	  protected:
+	    typename row<Matrix>::type  its_row;
+        };
+
+
+        template <typename Functor, typename Matrix>
+        struct mapped_col
+        {
+	    typedef typename Matrix::key_type   key_type;
+	    typedef typename Matrix::size_type  size_type;
+
+	    mapped_col(const mtl::matrix::map_view<Functor, Matrix>& view) : its_col(view.ref) {}
+	    mapped_col(const mtl::matrix::banded_view<Matrix>& view) : its_col(view.ref) {}
+
+	    size_type operator() (key_type const& key) const
+	    {
+		return its_col(key);
+	    }
+
+          protected:
+	    typename col<Matrix>::type  its_col;
+        };
+
+    } // namespace detail
+
+    template <typename Functor, typename Matrix>
+    struct row<mtl::matrix::map_view<Functor, Matrix> >
+    {
+	typedef detail::mapped_row<Functor, Matrix>   type;
+    };
+
+    template <typename Functor, typename Matrix>
+    struct col<mtl::matrix::map_view<Functor, Matrix> >
+    {
+	typedef detail::mapped_col<Functor, Matrix>   type;
+    };
+
+    template <typename Functor, typename Matrix>
+    struct const_value<mtl::matrix::map_view<Functor, Matrix> >
+    {
+	typedef detail::map_value<Functor, Matrix>  type;
+    };
+
+
+    // ================
+    // Range generators
+    // ================
+
+    // Use range_generator of original matrix
+    template <typename Tag, typename Functor, typename Matrix>
+    struct range_generator<Tag, mtl::matrix::map_view<Functor, Matrix> >
+	: public detail::referred_range_generator<mtl::matrix::map_view<Functor, Matrix>,
+						  range_generator<Tag, Matrix> >
+    {};
+
+    // To disambigue
+    template <typename Functor, typename Matrix>
+    struct range_generator<tag::major, mtl::matrix::map_view<Functor, Matrix> >
+	: public detail::referred_range_generator<mtl::matrix::map_view<Functor, Matrix>,
+						  range_generator<tag::major, Matrix> >
+    {};
+
+
+}} // mtl::traits
+
+
+namespace mtl { namespace matrix {
+
+template <typename Scaling, typename Matrix>
+struct scaled_view
+  : public map_view<tfunctor::scale<Scaling, typename Matrix::value_type>, Matrix>
+{
+    typedef tfunctor::scale<Scaling, typename Matrix::value_type>  functor_type;
+    typedef map_view<functor_type, Matrix>                         base;
+    typedef scaled_view                                            self;
+
+    scaled_view(const Scaling& scaling, const Matrix& matrix)
+      : base(functor_type(scaling), matrix)
+    {}
+
+    scaled_view(const Scaling& scaling, boost::shared_ptr<Matrix> p)
+      : base(functor_type(scaling), p)
+    {}
+
+#ifdef MTL_WITH_CPP11_MOVE
+    scaled_view (self&& that) : base(that) {}
+    scaled_view (const self& that) : base(that) {}
+#endif
+};
+
+// rscaled_view -- added by Hui Li
+template <typename Matrix, typename RScaling>
+struct rscaled_view
+  : public map_view<tfunctor::rscale<typename Matrix::value_type,RScaling>, Matrix>
+{
+    typedef tfunctor::rscale<typename Matrix::value_type, RScaling>  functor_type;
+    typedef map_view<functor_type, Matrix>                          base;
+    typedef rscaled_view                                            self;
+
+    rscaled_view(const Matrix& matrix, const RScaling& rscaling)
+      : base(functor_type(rscaling),matrix)
+    {}
+
+    rscaled_view(boost::shared_ptr<Matrix> p, const RScaling& rscaling)
+      : base(functor_type(rscaling), p)
+    {}
+
+#ifdef MTL_WITH_CPP11_MOVE
+    rscaled_view (self&& that) : base(that) {}
+    rscaled_view (const self& that) : base(that) {}
+#endif
+};
+
+// divide_by_view -- added by Hui Li
+template <typename Matrix, typename Divisor>
+struct divide_by_view
+  : public map_view<tfunctor::divide_by<typename Matrix::value_type,Divisor>, Matrix>
+{
+    typedef tfunctor::divide_by<typename Matrix::value_type, Divisor>  functor_type;
+    typedef map_view<functor_type, Matrix>                             base;
+    typedef divide_by_view                                             self;
+
+    divide_by_view(const Matrix& matrix,const Divisor& div)
+      : base(functor_type(div), matrix)
+    {}
+
+    divide_by_view(boost::shared_ptr<Matrix> p, const Divisor& div)
+      : base(functor_type(div), p)
+    {}
+
+#ifdef MTL_WITH_CPP11_MOVE
+    divide_by_view (self&& that) : base(that) {}
+    divide_by_view (const self& that) : base(that) {}
+#endif
+};
+
+template <typename Matrix>
+struct conj_view
+  : public map_view<mtl::sfunctor::conj<typename Matrix::value_type>, Matrix>
+{
+    typedef mtl::sfunctor::conj<typename Matrix::value_type>            functor_type;
+    typedef map_view<functor_type, Matrix>                              base;
+    typedef conj_view                                                   self;
+
+    conj_view(const Matrix& matrix) : base(functor_type(), matrix) {}
+    conj_view(boost::shared_ptr<Matrix> p) : base(functor_type(), p) {}
+
+#ifdef MTL_WITH_CPP11_MOVE
+    conj_view (self&& that) : base(that) {}
+    conj_view (const self& that) : base(that) {}
+#endif
+};
+
+template <typename Matrix>
+struct imag_view
+  : public map_view<mtl::sfunctor::imag<typename Matrix::value_type>, Matrix>
+{
+    typedef mtl::sfunctor::imag<typename Matrix::value_type>            functor_type;
+    typedef map_view<functor_type, Matrix>                              base;
+    typedef imag_view                                                   self;
+
+    imag_view(const Matrix& matrix) : base(functor_type(), matrix) {}
+    imag_view(boost::shared_ptr<Matrix> p) : base(functor_type(), p) {}
+
+#ifdef MTL_WITH_CPP11_MOVE
+    imag_view (self&& that) : base(that) {}
+    imag_view (const self& that) : base(that) {}
+#endif
+};
+
+template <typename Matrix>
+struct negate_view
+  : public map_view<mtl::sfunctor::negate<typename Matrix::value_type>, Matrix>
+{
+    typedef mtl::sfunctor::negate<typename Matrix::value_type>            functor_type;
+    typedef map_view<functor_type, Matrix>                              base;
+    typedef negate_view                                                   self;
+
+    negate_view(const Matrix& matrix) : base(functor_type(), matrix) {}
+    negate_view(boost::shared_ptr<Matrix> p) : base(functor_type(), p) {}
+
+#ifdef MTL_WITH_CPP11_MOVE
+    negate_view (self&& that) : base(that) {}
+    negate_view (const self& that) : base(that) {}
+#endif
+};
+
+template <typename Matrix>
+struct real_view
+  : public map_view<mtl::sfunctor::real<typename Matrix::value_type>, Matrix>
+{
+    typedef mtl::sfunctor::real<typename Matrix::value_type>            functor_type;
+    typedef map_view<functor_type, Matrix>                              base;
+    typedef real_view                                                   self;
+
+    real_view(const Matrix& matrix) : base(functor_type(), matrix) {}
+    real_view(boost::shared_ptr<Matrix> p) : base(functor_type(), p) {}
+
+#ifdef MTL_WITH_CPP11_MOVE
+    real_view (self&& that) : base(that) {}
+    real_view (const self& that) : base(that) {}
+#endif
+};
+
+template <typename Scaling, typename Matrix>
+struct sub_matrix_t< mtl::matrix::scaled_view<Scaling, Matrix> >
+  : public sub_matrix_t< mtl::matrix::map_view<tfunctor::scale<Scaling, typename Matrix::value_type>,
+					       Matrix> >
+{};
+
+template <typename Matrix>
+struct sub_matrix_t< mtl::matrix::conj_view<Matrix> >
+  : public sub_matrix_t< mtl::matrix::map_view<sfunctor::conj<typename Matrix::value_type>, Matrix> >
+{};
+
+template <typename Matrix, typename RScaling>
+struct sub_matrix_t< mtl::matrix::rscaled_view<Matrix, RScaling> >
+  : public sub_matrix_t< mtl::matrix::map_view<tfunctor::rscale<typename Matrix::value_type, RScaling>,
+					       Matrix> >
+{};
+
+template <typename Matrix, typename Divisor>
+struct sub_matrix_t< mtl::matrix::divide_by_view<Matrix, Divisor> >
+  : public sub_matrix_t< mtl::matrix::map_view<tfunctor::divide_by<typename Matrix::value_type, Divisor>,
+					       Matrix> >
+{};
+
+
+}} // namespace mtl::matrix
+
+namespace mtl { namespace sfunctor {
+
+    template <typename Matrix>
+    struct conj_aux<Matrix, tag::matrix>
+    {
+	typedef matrix::conj_view<Matrix> result_type;
+
+	static inline result_type apply(const Matrix& matrix)
+	{
+	    return result_type(matrix);
+	}
+
+	result_type operator() (const Matrix& matrix) const
+	{
+	    return apply(matrix);
+	}
+    };
+
+}} // namespace mtl::sfunctor
+
+// Traits for specific views
+namespace mtl { namespace traits {
+
+template <typename Scaling, typename Matrix>
+struct row< mtl::matrix::scaled_view<Scaling, Matrix> >
+  : public row< mtl::matrix::map_view<tfunctor::scale<Scaling, typename Matrix::value_type>,
+				      Matrix> >
+{};
+
+template <typename Matrix>
+struct row< mtl::matrix::conj_view<Matrix> >
+  : public row< mtl::matrix::map_view<sfunctor::conj<typename Matrix::value_type>, Matrix> >
+{};
+
+template <typename Matrix, typename RScaling>
+struct row< mtl::matrix::rscaled_view<Matrix, RScaling> >
+  : public row< mtl::matrix::map_view<tfunctor::rscale<typename Matrix::value_type, RScaling>,
+				      Matrix> >
+{};
+
+template <typename Matrix, typename Divisor>
+struct row< mtl::matrix::divide_by_view<Matrix, Divisor> >
+  : public row< mtl::matrix::map_view<tfunctor::divide_by<typename Matrix::value_type, Divisor>,
+				      Matrix> >
+{};
+
+
+template <typename Scaling, typename Matrix>
+struct col< mtl::matrix::scaled_view<Scaling, Matrix> >
+  : public col< mtl::matrix::map_view<tfunctor::scale<Scaling, typename Matrix::value_type>,
+				      Matrix> >
+{};
+
+template <typename Matrix>
+struct col< mtl::matrix::conj_view<Matrix> >
+  : public col< mtl::matrix::map_view<sfunctor::conj<typename Matrix::value_type>, Matrix> >
+{};
+
+template <typename Matrix, typename RScaling>
+struct col< mtl::matrix::rscaled_view<Matrix, RScaling> >
+  : public col< mtl::matrix::map_view<tfunctor::rscale<typename Matrix::value_type, RScaling>,
+				      Matrix> >
+{};
+
+template <typename Matrix, typename Divisor>
+struct col< mtl::matrix::divide_by_view<Matrix, Divisor> >
+  : public col< mtl::matrix::map_view<tfunctor::divide_by<typename Matrix::value_type, Divisor>,
+				      Matrix> >
+{};
+
+
+
+
+
+template <typename Scaling, typename Matrix>
+struct const_value< mtl::matrix::scaled_view<Scaling, Matrix> >
+  : public const_value< mtl::matrix::map_view<tfunctor::scale<Scaling, typename Matrix::value_type>,
+					      Matrix> >
+{};
+
+template <typename Matrix>
+struct const_value< mtl::matrix::conj_view<Matrix> >
+  : public const_value< mtl::matrix::map_view<sfunctor::conj<typename Matrix::value_type>, Matrix> >
+{};
+
+template <typename Matrix, typename RScaling>
+struct const_value< mtl::matrix::rscaled_view<Matrix, RScaling> >
+  : public const_value< mtl::matrix::map_view<tfunctor::rscale<typename Matrix::value_type, RScaling>,
+					      Matrix> >
+{};
+
+template <typename Matrix, typename Divisor>
+struct const_value< mtl::matrix::divide_by_view<Matrix, Divisor> >
+  : public const_value< mtl::matrix::map_view<tfunctor::divide_by<typename Matrix::value_type, Divisor>,
+					      Matrix> >
+{};
+
+
+
+
+
+template <typename Tag, typename Scaling, typename Matrix>
+struct range_generator< Tag, mtl::matrix::scaled_view<Scaling, Matrix> >
+    : public range_generator< Tag, mtl::matrix::map_view<tfunctor::scale<Scaling, typename Matrix::value_type>,
+						    Matrix> >
+{};
+
+template <typename Tag, typename Matrix>
+struct range_generator< Tag, mtl::matrix::conj_view<Matrix> >
+    : public range_generator< Tag, mtl::matrix::map_view<sfunctor::conj<typename Matrix::value_type>, Matrix> >
+{};
+
+template <typename Tag, typename Matrix, typename RScaling>
+struct range_generator< Tag, mtl::matrix::rscaled_view<Matrix, RScaling> >
+    : public range_generator< Tag, mtl::matrix::map_view<tfunctor::rscale<typename Matrix::value_type, RScaling>,
+						    Matrix> >
+{};
+
+template <typename Tag, typename Matrix, typename Divisor>
+struct range_generator< Tag, mtl::matrix::divide_by_view<Matrix, Divisor> >
+    : public range_generator< Tag, mtl::matrix::map_view<tfunctor::divide_by<typename Matrix::value_type, Divisor>,
+						    Matrix> >
+{};
+
+
+
+template <typename Scaling, typename Matrix>
+struct range_generator< tag::major, mtl::matrix::scaled_view<Scaling, Matrix> >
+    : public range_generator< tag::major, mtl::matrix::map_view<tfunctor::scale<Scaling, typename Matrix::value_type>,
+						    Matrix> >
+{};
+
+template <typename Matrix>
+struct range_generator< tag::major, mtl::matrix::conj_view<Matrix> >
+    : public range_generator< tag::major, mtl::matrix::map_view<sfunctor::conj<typename Matrix::value_type>, Matrix> >
+{};
+
+template <typename Matrix, typename RScaling>
+struct range_generator< tag::major, mtl::matrix::rscaled_view<Matrix, RScaling> >
+    : public range_generator< tag::major, mtl::matrix::map_view<tfunctor::rscale<typename Matrix::value_type, RScaling>,
+						    Matrix> >
+{};
+
+template <typename Matrix, typename Divisor>
+struct range_generator< tag::major, mtl::matrix::divide_by_view<Matrix, Divisor> >
+    : public range_generator< tag::major, mtl::matrix::map_view<tfunctor::divide_by<typename Matrix::value_type, Divisor>,
+						    Matrix> >
+{};
+
+
+
+}} // mtl::traits
+
+
+#endif // MTL_MAP_VIEW_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/mapped_inserter.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_MAPPED_INSERTER_INCLUDE
+#define MTL_MATRIX_MAPPED_INSERTER_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/operation/update.hpp>
+
+namespace mtl { namespace matrix {
+
+/// Inserter with shifted row and column indices
+/** The main work is performed by the underlying base inserter whose type is given as template
+    argument. **/
+template <typename BaseInserter, typename Mapper >
+class mapped_inserter
+{
+  public:
+    typedef mapped_inserter                                self;
+    typedef typename BaseInserter::matrix_type              matrix_type;
+    typedef typename Collection<matrix_type>::size_type     size_type;
+    typedef operations::update_proxy<BaseInserter, size_type>   proxy_type;
+
+    /// Constructor with matrix \p A, the mapping, and the slot size
+    mapped_inserter(matrix_type& A, Mapper& map, size_type slot_size= 0)
+      : ins(A, slot_size), map(map) {}
+
+  private:
+    struct bracket_proxy
+    {
+	bracket_proxy(BaseInserter& ref, Mapper& map, size_type row)
+          : ref(ref),map(map), row(row) {}
+
+	proxy_type operator[](size_type col)
+	{   return proxy_type(ref, row, map.col(col)); }
+
+	BaseInserter&   ref;
+        Mapper&         map;
+	size_type row;
+    };
+
+  public:
+    /// To be used in ins[r][c] << value;
+    bracket_proxy operator[] (size_type row)
+    {	return bracket_proxy(ins, map, map.row(row));    }
+
+    /// To be used in ins(r, c) << value;
+    proxy_type operator() (size_type row, size_type col)
+    {	return proxy_type(ins, map.row(row),map.col(col));    }
+
+    // update, modify and operator<< are used from BaseInserter
+
+private:
+    BaseInserter ins;
+    Mapper&     map;
+};
+
+template< typename BaseInserter, typename Mapper, typename Elt, typename Parameters >
+mapped_inserter< BaseInserter, Mapper >& operator<<(mapped_inserter< BaseInserter, Mapper >& minserter,
+    const compressed2D< Elt, Parameters >& rhs)
+{
+
+  using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
+  namespace traits= mtl::traits;
+  typedef compressed2D< Elt, Parameters >   Matrix;
+
+  typename traits::row<Matrix>::type                                 row(rhs);
+  typename traits::col<Matrix>::type                                 col(rhs);
+  typename traits::const_value<Matrix>::type                         value(rhs);
+
+  typedef typename traits::range_generator<major, Matrix>::type      cursor_type;
+  typedef typename traits::range_generator<nz, cursor_type>::type    icursor_type;
+
+  for (cursor_type cursor = begin<major>(rhs), cend = end<major>(rhs); cursor != cend; ++cursor)
+    for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor)
+      minserter[row(*icursor)][col(*icursor)]<< value(*icursor);
+
+  return minserter;
+}
+
+} // namespace matrix
+
+} // namespace mtl
+
+#endif // MTL_MATRIX_MAPPED_INSERTER_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/mat_expr.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_MAT_EXPR_INCLUDE
+#define MTL_MAT_EXPR_INCLUDE
+
+namespace mtl { namespace matrix {
+
+/// Base class for CRTP with matrices
+template <typename Matrix>
+struct mat_expr
+{
+    typedef Matrix   ref_type;
+};
+
+
+/// Base class for CRTP with dense matrices
+template <typename Matrix>
+struct dmat_expr
+    : public mat_expr<Matrix>
+{
+    typedef mat_expr<Matrix> base;
+};
+
+
+/// Base class for CRTP with sparse matrices
+template <typename Matrix>
+struct smat_expr
+    : public mat_expr<Matrix>
+{
+    typedef mat_expr<Matrix> base;
+};
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MAT_EXPR_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/mat_mat_asgn_expr.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_MAT_MAT_ASGN_EXPR_INCLUDE
+#define MTL_MATRIX_MAT_MAT_ASGN_EXPR_INCLUDE
+
+namespace mtl { namespace matrix {
+
+
+// Currently only dummy to be used in traits::eval_dense
+template <typename E1, typename E2>
+struct mat_mat_asgn_expr {};
+
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_MAT_MAT_ASGN_EXPR_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/mat_mat_ele_times_expr.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_MAT_MAT_ELE_TIMES_EXPR_INCLUDE
+#define MTL_MAT_MAT_ELE_TIMES_EXPR_INCLUDE
+
+#include <boost/shared_ptr.hpp>
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/matrix/mat_mat_op_expr.hpp>
+#include <boost/numeric/mtl/operation/sfunctor.hpp>
+#include <boost/numeric/linear_algebra/identity.hpp>
+
+
+namespace mtl { namespace matrix {
+
+template <typename E1, typename E2>
+struct mat_mat_ele_times_expr
+    : public mat_mat_op_expr< E1, E2, mtl::sfunctor::times<typename E1::value_type, typename E2::value_type> >,
+      public mat_expr< mat_mat_ele_times_expr<E1, E2> >
+{
+    typedef mat_mat_op_expr< E1, E2, mtl::sfunctor::times<typename E1::value_type, typename E2::value_type> > op_base;
+    typedef mat_expr< mat_mat_ele_times_expr<E1, E2> >                                                        crtp_base;
+    typedef mat_mat_ele_times_expr               self;
+    typedef E1                                   first_argument_type ;
+    typedef E2                                   second_argument_type ;
+    typedef typename E1::orientation             orientation;
+    typedef mtl::non_fixed::dimensions           dim_type;
+    typedef typename E1::key_type                key_type;
+
+
+    mat_mat_ele_times_expr( E1 const& v1, E2 const& v2 )
+	: op_base( v1, v2 ), crtp_base(*this), first(v1), second(v2)
+    {}
+
+    first_argument_type const&  first ;
+    second_argument_type const& second ;
+};
+
+template <typename E1, typename E2>
+std::size_t inline num_rows(const mat_mat_ele_times_expr<E1, E2>& expr)
+{ return num_rows(expr.first); }
+
+template <typename E1, typename E2>
+std::size_t inline num_cols(const mat_mat_ele_times_expr<E1, E2>& expr)
+{ return num_cols(expr.second); }
+
+template <typename E1, typename E2>
+std::size_t inline size(const mat_mat_ele_times_expr<E1, E2>& expr)
+{ return num_rows(expr) * num_cols(expr); }
+
+}} // Namespace mtl::matrix
+
+#endif // MTL_MAT_MAT_ELE_TIMES_EXPR_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/mat_mat_minus_expr.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_MAT_MAT_MINUS_EXPR_INCLUDE
+#define MTL_MAT_MAT_MINUS_EXPR_INCLUDE
+
+
+#include <boost/numeric/mtl/matrix/mat_mat_op_expr.hpp>
+#include <boost/numeric/mtl/vector/vec_vec_pmop_expr.hpp>
+#include <boost/numeric/mtl/operation/sfunctor.hpp>
+
+namespace mtl {namespace matrix {
+
+template <typename E1, typename E2>
+struct mat_mat_minus_expr
+  : public mat_mat_op_expr< E1, E2, mtl::sfunctor::minus<typename E1::value_type, typename E2::value_type> >,
+    public mat_expr< mat_mat_minus_expr<E1, E2> >
+{
+    typedef mat_mat_op_expr< E1, E2, mtl::sfunctor::minus<typename E1::value_type, typename E2::value_type> > op_base;
+    typedef mat_expr< mat_mat_minus_expr<E1, E2> >                                                       crtp_base;
+    typedef E1                                   first_argument_type ;
+    typedef E2                                   second_argument_type ;
+
+    mat_mat_minus_expr( E1 const& v1, E2 const& v2 )
+      : op_base( v1, v2 ), crtp_base(*this), first(v1), second(v2)
+    {}
+
+    first_argument_type const&  first ;
+    second_argument_type const& second ;
+};
+
+template <typename E1, typename E2>
+struct mv_mv_minus_expr
+  : mat_mat_minus_expr<E1, E2>
+{
+    typedef mat_mat_minus_expr< E1, E2 > base;
+    typedef typename E1::vector_type V1;
+    typedef typename E2::vector_type V2;
+    typedef mtl::vector::vec_vec_pmop_expr< V1, V2, mtl::sfunctor::minus<typename V1::value_type, typename V2::value_type> > vector_type;
+
+    mv_mv_minus_expr( E1 const& v1, E2 const& v2 )
+      : base( v1, v2 )
+    {}
+
+    vector_type vector(std::size_t c) const { return vector_type(this->first.vector(c), this->second.vector(c)); }
+};
+
+}} // Namespace mtl::matrix
+
+
+#endif // MTL_MAT_MAT_MINUS_EXPR_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/mat_mat_op_expr.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_MAT_MAT_OP_EXPR_INCLUDE
+#define MTL_MAT_MAT_OP_EXPR_INCLUDE
+
+#include <boost/numeric/mtl/matrix/mat_expr.hpp>
+#include <boost/numeric/mtl/matrix/crtp_base_matrix.hpp>
+
+namespace mtl { namespace matrix {
+
+
+template <typename E1, typename E2, typename SFunctor>
+struct mat_mat_op_expr
+    : public const_crtp_matrix_bracket< mat_mat_op_expr<E1, E2, SFunctor>,
+					typename SFunctor::result_type,
+					typename E1::size_type >
+   // : public mat_expr< mat_mat_op_expr<E1, E2, SFunctor> >
+{
+    // typedef mat_expr< mat_mat_op_expr<E1, E2, SFunctor> > expr_base;
+    typedef mat_mat_op_expr                               self;
+
+    typedef typename SFunctor::result_type      value_type;
+
+    // temporary solution
+    typedef typename E1::size_type               size_type;
+
+    typedef value_type                           const_dereference_type ;
+
+    typedef E1                                   first_argument_type ;
+    typedef E2                                   second_argument_type ;
+
+    mat_mat_op_expr( first_argument_type const& v1, second_argument_type const& v2 )
+	: // expr_base( *this ),
+	  first( v1 ), second( v2 )
+    {
+#if 0
+	first.delay_assign();
+	second.delay_assign();
+#endif
+    }
+
+    void delay_assign() const {}
+
+    void check_shape() const {} // consistency of shapes depend on operation
+
+    const_dereference_type operator() (size_type row, size_type col) const
+    {
+        return SFunctor::apply( first(row, col), second(row, col) ) ;
+	// return SFunctor::apply( first[row][col], second[row][col] ) ;
+    }
+
+#if 0
+    template <typename> friend size_type size(const self&);
+    template <typename> friend size_type num_rows(const self&);
+    template <typename> friend size_type num_cols(const self&);
+#endif
+
+    first_argument_type const&  first ;
+    second_argument_type const& second ;
+};
+
+
+template <typename E1, typename E2, typename SFunctor>
+typename mat_mat_op_expr<E1, E2, SFunctor>::size_type
+inline size(mat_mat_op_expr<E1, E2, SFunctor> const& expr)
+{
+    expr.check_shape();
+    return size(expr.first) ;
+}
+
+template <typename E1, typename E2, typename SFunctor>
+typename mat_mat_op_expr<E1, E2, SFunctor>::size_type
+inline num_rows(mat_mat_op_expr<E1, E2, SFunctor> const& expr)
+{
+    expr.check_shape();
+    return num_rows(expr.first) ;
+}
+
+template <typename E1, typename E2, typename SFunctor>
+typename mat_mat_op_expr<E1, E2, SFunctor>::size_type
+inline num_cols(mat_mat_op_expr<E1, E2, SFunctor> const& expr)
+{
+    expr.check_shape();
+    return num_cols(expr.first) ;
+}
+
+
+}} // namespace mtl
+
+#endif // MTL_MAT_MAT_OP_EXPR_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/mat_mat_plus_expr.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_MAT_MAT_PLUS_EXPR_INCLUDE
+#define MTL_MAT_MAT_PLUS_EXPR_INCLUDE
+
+#include <boost/numeric/mtl/matrix/mat_mat_op_expr.hpp>
+#include <boost/numeric/mtl/vector/vec_vec_pmop_expr.hpp>
+#include <boost/numeric/mtl/operation/sfunctor.hpp>
+
+namespace mtl { namespace matrix {
+
+template <typename E1, typename E2>
+struct mat_mat_plus_expr
+  : public mat_mat_op_expr< E1, E2, mtl::sfunctor::plus<typename E1::value_type, typename E2::value_type> >,
+    public mat_expr< mat_mat_plus_expr<E1, E2> >
+{
+    typedef mat_mat_op_expr< E1, E2, mtl::sfunctor::plus<typename E1::value_type, typename E2::value_type> > op_base;
+    typedef mat_expr< mat_mat_plus_expr<E1, E2> >                                                       crtp_base;
+    typedef E1                                   first_argument_type ;
+    typedef E2                                   second_argument_type ;
+
+    mat_mat_plus_expr( E1 const& v1, E2 const& v2 )
+      : op_base( v1, v2 ), crtp_base(*this), first(v1), second(v2)
+    {}
+
+    first_argument_type const&  first ;
+    second_argument_type const& second ;
+};
+
+// Same as mat_mat_plus_expr for pair of dense matrix expressions
+// Future versions will probably provide more efficient implementations for it
+template <typename E1, typename E2>
+struct dmat_dmat_plus_expr
+  : public mat_mat_op_expr< E1, E2, mtl::sfunctor::plus<typename E1::value_type, typename E2::value_type> >
+{
+    typedef mat_mat_op_expr< E1, E2, mtl::sfunctor::plus<typename E1::value_type, typename E2::value_type> > base;
+    dmat_dmat_plus_expr( E1 const& v1, E2 const& v2 )
+      : base( v1, v2 )
+    {}
+};
+
+template <typename E1, typename E2>
+struct mv_mv_plus_expr
+  : mat_mat_plus_expr<E1, E2>
+{
+    typedef mat_mat_plus_expr< E1, E2 > base;
+    typedef typename E1::vector_type V1;
+    typedef typename E2::vector_type V2;
+    typedef mtl::vector::vec_vec_pmop_expr< V1, V2, mtl::sfunctor::plus<typename V1::value_type, typename V2::value_type> > vector_type;
+
+    mv_mv_plus_expr( E1 const& v1, E2 const& v2 )
+      : base( v1, v2 )
+    {}
+
+    vector_type vector(std::size_t c) const { return vector_type(this->first.vector(c), this->second.vector(c)); }
+};
+
+
+
+
+}} // Namespace mtl::matrix
+
+#endif // MTL_MAT_MAT_PLUS_EXPR_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/mat_mat_times_expr.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_MAT_MAT_TIMES_EXPR_INCLUDE
+#define MTL_MAT_MAT_TIMES_EXPR_INCLUDE
+
+#include <boost/mpl/if.hpp>
+#include <boost/mpl/and.hpp>
+#include <boost/shared_ptr.hpp>
+#include <boost/type_traits/is_base_of.hpp>
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/matrix/mat_mat_op_expr.hpp>
+#include <boost/numeric/mtl/operation/sfunctor.hpp>
+#include <boost/numeric/mtl/operation/compute_factors.hpp>
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/mtl/matrix/parameter.hpp>
+#include <boost/numeric/mtl/matrix/dense2D.hpp>
+#include <boost/numeric/mtl/matrix/compressed2D.hpp>
+#include <boost/numeric/mtl/concept/std_concept.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+
+
+namespace mtl { namespace matrix {
+
+template <typename E1, typename E2>
+struct mat_mat_times_expr
+  : public mat_mat_op_expr< E1, E2, mtl::sfunctor::times<typename Collection<E1>::value_type, typename Collection<E2>::value_type> >,
+    public mat_expr< mat_mat_times_expr<E1, E2> >
+{
+    typedef mat_mat_op_expr< E1, E2, mtl::sfunctor::times<typename Collection<E1>::value_type, typename Collection<E2>::value_type> > op_base;
+    typedef mat_expr< mat_mat_times_expr<E1, E2> >                                                       crtp_base;
+    typedef mat_mat_times_expr                   self;
+    typedef E1                                   first_argument_type ;
+    typedef E2                                   second_argument_type ;
+    typedef typename E1::orientation             orientation;
+    typedef mtl::non_fixed::dimensions           dim_type;
+    typedef typename E1::key_type                key_type;
+
+
+    typedef typename Collection<E1>::value_type  first_value_type;
+    typedef typename Collection<E2>::value_type  second_value_type;
+    typedef typename Multiplicable<first_value_type, second_value_type>::result_type result_value_type;
+
+#if 0 // Just an idea
+    typedef typename boost::mpl::if_<
+                boost::mpl::and_<
+	            boost::is_base_of<tag::sparse, typename traits::category<E1>::type>
+	          , boost::is_base_of<tag::sparse, typename traits::category<E2>::type>
+		>
+	      , compressed2D<result_value_type>
+	      , dense2D<result_value_type, parameters<> >
+            >::type                                       evaluated_result_type;
+
+    // Convert into matrix
+
+    operator evaluated_result_type() const
+    {
+	return evaluated_result_type(first * second);
+    }
+#endif
+
+    mat_mat_times_expr( E1 const& v1, E2 const& v2 )
+      : op_base( v1, v2 ), crtp_base(*this), first(v1), second(v2)
+    {}
+
+    // To prevent that cout << A * B prints the element-wise product, suggestion by Hui Li
+    // It is rather inefficient, esp. for multiple products (complexity increases with the number of arguments :-!)
+    //    or sparse matrices.
+    // Better compute your product first and print it then when compute time is an issue,
+    // this is ONLY for convenience.
+    result_value_type
+    operator()(std::size_t r, std::size_t c) const
+    {
+	using math::zero;
+	MTL_THROW_IF(num_cols(first) != num_rows(second), incompatible_size());
+
+	result_value_type ref, sum(zero(ref));
+	for (std::size_t i= 0; i < num_cols(first); i++)
+	    sum+= first(r, i) * second(i, c);
+	return sum;
+    }
+
+
+    result_value_type
+    operator()(std::size_t r, std::size_t c)
+    {
+	return (*const_cast<const self*>(this))(r, c);
+    }
+
+    first_argument_type const&  first ;
+    second_argument_type const& second ;
+};
+
+template <typename E1, typename E2>
+std::size_t inline num_rows(const mat_mat_times_expr<E1, E2>& expr)
+{ return num_rows(expr.first); }
+
+template <typename E1, typename E2>
+std::size_t inline num_cols(const mat_mat_times_expr<E1, E2>& expr)
+{ return num_cols(expr.second); }
+
+template <typename E1, typename E2>
+std::size_t inline size(const mat_mat_times_expr<E1, E2>& expr)
+{ return num_rows(expr) * num_cols(expr); }
+
+}} // Namespace mtl::matrix
+
+#endif // MTL_MAT_MAT_TIMES_EXPR_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/mat_negate_expr.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_MAT_NEGATE_EXPR_INCLUDE
+#define MTL_MAT_NEGATE_EXPR_INCLUDE
+
+#include <boost/numeric/mtl/matrix/map_view.hpp>
+
+namespace mtl { namespace matrix {
+
+template <typename E1>
+inline negate_view< E1 >
+operator- (const mat_expr<E1>& e1)
+{
+    return negate_view< E1 >(static_cast<const E1&>(e1));
+}
+
+} } // Namespace mtl::matrix
+
+#endif
+

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/morton_dense.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_MORTON_DENSE_INCLUDE
+#define MTL_MORTON_DENSE_INCLUDE
+
+#include <boost/type_traits/is_same.hpp>
+#include <boost/mpl/if.hpp>
+
+#include <boost/numeric/mtl/utility/common_include.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/matrix/parameter.hpp>
+#include <boost/numeric/mtl/matrix/crtp_base_matrix.hpp>
+#include <boost/numeric/mtl/matrix/base_sub_matrix.hpp>
+#include <boost/numeric/mtl/detail/contiguous_memory_block.hpp>
+#include <boost/numeric/mtl/detail/dilated_int.hpp>
+#include <boost/numeric/mtl/utility/iterator_adaptor.hpp>
+#include <boost/numeric/mtl/operation/set_to_zero.hpp>
+#include <boost/numeric/mtl/matrix/mat_expr.hpp>
+#include <boost/numeric/mtl/operation/print_matrix.hpp>
+#include <boost/numeric/mtl/operation/compute_factors.hpp>
+#include <boost/numeric/mtl/operation/clone.hpp>
+
+#ifdef MTL_WITH_INITLIST
+# include <initializer_list>
+#endif
+
+
+namespace mtl { namespace matrix {
+
+// Helper type
+struct morton_dense_sub_ctor {};
+
+template <std::size_t BitMask>
+struct morton_dense_key
+{
+    typedef std::size_t                               size_type;
+    typedef dilated_int<std::size_t, BitMask, true>   dilated_row_t;
+    typedef dilated_int<std::size_t, ~BitMask, true>  dilated_col_t;
+    typedef morton_dense_key                          self;
+
+    morton_dense_key(size_type my_row, size_type my_col)
+	: my_row(my_row), my_col(my_col), dilated_row(my_row), dilated_col(my_col)
+    {}
+
+    bool operator== (self const& x) const
+    {
+	return my_row == x.my_row && my_col == x.my_col;
+    }
+
+    bool operator!= (self const& x)
+    {
+	return !(*this == x);
+    }
+
+    size_type row() const
+    {
+	return my_row;
+    }
+
+    size_type col() const
+    {
+	return my_col;
+    }
+
+    self& advance_row(int row_inc)
+    {
+	dilated_row.advance(row_inc);
+	// potential addition of signed and unsigned
+	my_row+= row_inc;
+	return *this;
+    }
+
+    self& advance_col(int col_inc)
+    {
+	dilated_col.advance(col_inc);
+	// potential addition of signed and unsigned
+	my_col+= col_inc;
+	return *this;
+    }
+
+    self& advance(int row_inc, int col_inc)
+    {
+	advance_row(row_inc);
+	advance_col(col_inc);
+	return *this;
+    }
+
+public:
+    size_type                    my_row, my_col;
+    dilated_row_t                dilated_row;
+    dilated_col_t                dilated_col;
+};
+
+template <std::size_t BitMask>
+struct morton_dense_el_cursor
+  : public morton_dense_key<BitMask>
+{
+    typedef std::size_t                               size_type;
+    typedef dilated_int<std::size_t, ~BitMask, true>  dilated_col_t;
+    typedef morton_dense_el_cursor                    self;
+    typedef morton_dense_key<BitMask>                 base;
+    typedef base                                      key_type;
+
+    morton_dense_el_cursor(size_type my_row, size_type my_col, size_type num_cols)
+	: base(my_row, my_col), num_cols(num_cols)
+    {}
+
+    self& operator++ ()
+    {
+	++this->my_col; ++this->dilated_col;
+	if (this->my_col == num_cols) {
+	    this->my_col= 0; this->dilated_col= dilated_col_t(0);
+	    ++this->my_row; ++this->dilated_row;
+	}
+	return *this;
+    }
+
+    base& operator* ()
+    {
+	return *this;
+    }
+
+    const base& operator* () const
+    {
+	return *this;
+    }
+
+protected:
+    size_t                       num_cols;
+};
+
+template <std::size_t BitMask>
+struct morton_dense_row_cursor
+    : public morton_dense_key<BitMask>
+{
+    typedef std::size_t                               size_type;
+    typedef morton_dense_row_cursor                   self;
+    typedef morton_dense_key<BitMask>                 base;
+    typedef base                                      key_type;
+
+    morton_dense_row_cursor(size_type my_row, size_type my_col)
+	: base(my_row, my_col)
+    {}
+
+    self& operator++ ()
+    {
+	++this->my_row; ++this->dilated_row;
+	return *this;
+    }
+
+    self& operator+=(int inc)
+    {
+	this->advance_row(inc);
+	return *this;
+    }
+
+    self& operator-- ()
+    {
+	--this->my_row; --this->dilated_row;
+	return *this;
+    }
+
+    self& operator-=(int dec)
+    {
+	this->advance_row(-dec);
+	return *this;
+    }
+
+    self operator+ (int inc) const
+    {
+	self tmp(*this);
+	tmp.advance_row(inc);
+	return tmp;
+    }
+
+    base& operator* ()
+    {
+	return *this;
+    }
+
+    const base& operator* () const
+    {
+	return *this;
+    }
+};
+
+template <std::size_t BitMask>
+struct morton_dense_col_cursor
+    : public morton_dense_key<BitMask>
+{
+    typedef std::size_t                               size_type;
+    typedef morton_dense_col_cursor                   self;
+    typedef morton_dense_key<BitMask>                 base;
+    typedef base                                      key_type;
+
+    morton_dense_col_cursor(size_type my_row, size_type my_col)
+	: base(my_row, my_col)
+    {}
+
+    self& operator++ ()
+    {
+	++this->my_col; ++this->dilated_col;
+	return *this;
+    }
+
+    self& operator+=(int inc)
+    {
+	this->advance_col(inc);
+	return *this;
+    }
+
+    self& operator-- ()
+    {
+	--this->my_col; --this->dilated_col;
+	return *this;
+    }
+
+    self& operator-=(int dec)
+    {
+	this->advance_col(-dec);
+	return *this;
+    }
+
+    self operator+ (int inc) const
+    {
+	self tmp(*this);
+	tmp.advance_col(inc);
+	return tmp;
+    }
+
+    base& operator* ()
+    {
+	return *this;
+    }
+
+    const base& operator* () const
+    {
+	return *this;
+    }
+};
+
+
+template <typename Matrix>
+struct morton_dense_row_const_iterator
+    : utilities::const_iterator_adaptor<typename mtl::traits::const_value<Matrix>::type, morton_dense_row_cursor<Matrix::mask>,
+					typename Matrix::value_type>
+{
+    static const std::size_t                          mask= Matrix::mask;
+    typedef morton_dense_row_cursor<mask>               cursor_type;
+    typedef typename mtl::traits::const_value<Matrix>::type  map_type;
+    typedef typename Matrix::value_type                 value_type;
+    typedef typename Matrix::size_type                  size_type;
+    typedef utilities::const_iterator_adaptor<map_type, cursor_type, value_type> base;
+
+    morton_dense_row_const_iterator(const Matrix& matrix, size_type row, size_type col)
+	: base(map_type(matrix), cursor_type(row, col))
+    {}
+};
+
+
+template <typename Matrix>
+struct morton_dense_row_iterator
+    : utilities::iterator_adaptor<typename mtl::traits::value<Matrix>::type, morton_dense_row_cursor<Matrix::mask>,
+				  typename Matrix::value_type>
+{
+    static const std::size_t                          mask= Matrix::mask;
+    typedef morton_dense_row_cursor<mask>               cursor_type;
+    typedef typename mtl::traits::value<Matrix>::type   map_type;
+    typedef typename Matrix::value_type                 value_type;
+    typedef typename Matrix::size_type                  size_type;
+    typedef utilities::iterator_adaptor<map_type, cursor_type, value_type> base;
+
+    morton_dense_row_iterator(Matrix& matrix, size_type row, size_type col)
+	:  base(map_type(matrix), cursor_type(row, col))
+    {}
+};
+
+
+template <typename Matrix>
+struct morton_dense_col_const_iterator
+    : utilities::const_iterator_adaptor<typename mtl::traits::const_value<Matrix>::type, morton_dense_col_cursor<Matrix::mask>,
+					typename Matrix::value_type>
+{
+    static const std::size_t                          mask= Matrix::mask;
+    typedef morton_dense_col_cursor<mask>               cursor_type;
+    typedef typename mtl::traits::const_value<Matrix>::type  map_type;
+    typedef typename Matrix::value_type                 value_type;
+    typedef typename Matrix::size_type                  size_type;
+    typedef utilities::const_iterator_adaptor<map_type, cursor_type, value_type> base;
+
+    morton_dense_col_const_iterator(const Matrix& matrix, size_type row, size_type col)
+	: base(map_type(matrix), cursor_type(row, col))
+    {}
+};
+
+
+template <typename Matrix>
+struct morton_dense_col_iterator
+    : utilities::iterator_adaptor<typename mtl::traits::value<Matrix>::type, morton_dense_col_cursor<Matrix::mask>,
+				  typename Matrix::value_type>
+{
+    static const std::size_t                          mask= Matrix::mask;
+    typedef morton_dense_col_cursor<mask>               cursor_type;
+    typedef typename mtl::traits::value<Matrix>::type   map_type;
+    typedef typename Matrix::value_type                 value_type;
+    typedef typename Matrix::size_type                  size_type;
+    typedef utilities::iterator_adaptor<map_type, cursor_type, value_type> base;
+
+    morton_dense_col_iterator(Matrix& matrix, size_type row, size_type col)
+      : base(map_type(matrix), cursor_type(row, col))  {}
+};
+
+
+
+/// Dense Morton-order matrix
+template <typename Elt, std::size_t BitMask, typename Parameters = mtl::matrix::parameters<> >
+class morton_dense
+  : public base_sub_matrix<Elt, Parameters>,
+    public mtl::detail::contiguous_memory_block<Elt, false>,
+    public crtp_base_matrix< morton_dense<Elt, BitMask, Parameters>, Elt, std::size_t >,
+    public mat_expr< morton_dense<Elt, BitMask, Parameters> >
+{
+    typedef morton_dense                                               self;
+    typedef base_sub_matrix<Elt, Parameters>                           super;
+    typedef mtl::detail::contiguous_memory_block<Elt, false>           memory_base;
+    typedef crtp_matrix_assign< self, Elt, std::size_t >               assign_base;
+    typedef mat_expr< morton_dense<Elt, BitMask, Parameters> >         expr_base;
+
+  public:
+
+    typedef Parameters                        parameters;
+    typedef typename Parameters::orientation  orientation;
+    typedef typename Parameters::index        index_type;
+    typedef typename Parameters::dimensions   dim_type;
+    typedef Elt                               value_type;
+    typedef const value_type&                 const_reference;
+    typedef value_type&                       reference;
+    typedef typename parameters::size_type    size_type;
+    const static std::size_t                  mask= BitMask;
+
+    //  implement cursor for morton matrix, somewhere
+    //  also, morton indexer?
+
+    typedef morton_dense_key<BitMask>          key_type;
+    typedef morton_dense_el_cursor<BitMask>    el_cursor_type;
+
+    typedef dilated_int<std::size_t, BitMask, true>   dilated_row_t;
+    typedef dilated_int<std::size_t, ~BitMask, true>  dilated_col_t;
+
+  protected:
+
+    // ranges of rows and columns
+    dilated_row_t            my_begin_row, my_end_row;
+    dilated_col_t            my_begin_col, my_end_col;
+
+    // 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);
+	my_begin_row= begin_r; my_end_row= end_r;
+	my_begin_col= begin_c; my_end_col= 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);
+    }
+
+    void init(size_type num_rows, size_type num_cols)
+    {
+	set_ranges(num_rows, num_cols);
+	// set_to_zero(*this);
+    }
+
+  public:
+    /// Default constructor
+    /** If compile time matrix size allocate memory. **/
+    morton_dense() : memory_base(memory_need(dim_type().num_rows(), dim_type().num_cols()))
+    {
+	init(dim_type().num_rows(), dim_type().num_cols());
+    }
+
+    /// Construction from run-time dimension type
+    explicit morton_dense(mtl::non_fixed::dimensions d)
+	: memory_base(memory_need(d.num_rows(), d.num_cols()))
+    {
+	init(d.num_rows(), d.num_cols());
+    }
+
+    /// Construction of matrix of dimension \p num_rows by \p num_cols
+    morton_dense(size_type num_rows, size_type num_cols)
+	: memory_base(memory_need(num_rows, num_cols))
+    {
+	init(num_rows, num_cols);
+    }
+
+    /// Construction of matrix with dimension \p d using pointer \p a to external data
+    explicit morton_dense(mtl::non_fixed::dimensions d, value_type* a)
+      : memory_base(a, memory_need(d.num_rows(), d.num_cols()))
+    {
+	set_ranges(d.num_rows(), d.num_cols());
+    }
+
+    /// Construction of \p num_rows by \p num_cols matrix with pointer \p a to external data
+    explicit morton_dense(size_type num_rows, size_type num_cols, value_type* a)
+      : memory_base(a, memory_need(num_rows, num_cols))
+    {
+	set_ranges(num_rows, num_cols);
+    }
+
+    /// Construction of matrix with static dimension using pointer \p a to external data
+    explicit morton_dense(value_type* a)
+	: memory_base(a, memory_need(dim_type().num_rows(), dim_type().num_cols()))
+    {
+	BOOST_ASSERT((dim_type::is_static));
+	set_ranges(dim_type().num_rows(), dim_type().num_cols());
+    }
+
+    /// Copy constructor
+    morton_dense(const self& m)
+      : super(m), memory_base(m)
+    {
+	set_ranges(m.num_rows(), m.num_cols());
+    }
+
+    /// Clone constructor
+    explicit morton_dense(const self& m, clone_ctor)
+      : memory_base(m, clone_ctor())
+    {
+	init(m.num_rows(), m.num_cols());
+	*this= m;
+    }
+
+    /// Templated copy constructor
+    template <typename MatrixSrc>
+    explicit morton_dense(const MatrixSrc& src)
+      : memory_base(memory_need(dim_type().num_rows(), dim_type().num_cols()))
+    {
+	init(dim_type().num_rows(), dim_type().num_cols());
+	*this= src;
+    }
+
+#if defined(MTL_WITH_INITLIST) && defined(MTL_WITH_AUTO) && defined(MTL_WITH_RANGEDFOR)
+    /// Constructor for initializer list \p values
+    template <typename Value2>
+    morton_dense(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->num_rows(), this->num_cols());
+	*this= values;
+    }
+#endif
+
+
+    /// Construct a sub-matrix as a view
+    explicit morton_dense(self& matrix, morton_dense_sub_ctor,
+			  size_type begin_r, size_type end_r, size_type begin_c, size_type end_c)
+      : memory_base(matrix.data, memory_need(end_r - begin_r, end_c - begin_c), true) // View constructor
+    {
+	matrix.check_ranges(begin_r, end_r, begin_c, end_c);
+
+	if (begin_r >= end_r || begin_c >= end_c) {
+	    set_ranges(0, 0);
+	    return;
+	}
+
+	// Check whether sub-matrix is contigous memory block
+	// by comparing the address of the last and the first element in the entire and the sub-matrix
+	MTL_DEBUG_THROW_IF(&matrix[end_r-1][end_c-1] - &matrix[begin_r][begin_c]
+			   != &matrix[end_r-begin_r-1][end_c-begin_c-1] - &matrix[0][0],
+			   range_error("This sub-matrix cannot be used because it is split in memory"));
+	// Check with David if this is a sufficient condition (it is a necessary at least)
+
+	dilated_row_t  dilated_row(begin_r);
+	dilated_col_t  dilated_col(begin_c);
+
+	// Set new start address within masked matrix
+	this->data += dilated_row.dilated_value() + dilated_col.dilated_value();
+	set_ranges(end_r - begin_r, end_c - begin_c);
+    }
+
+    /// Change dimension to \p num_rows by \p num_cols
+    void change_dim(size_type num_rows, size_type num_cols)
+    {
+	set_ranges(num_rows, num_cols);
+	this->realloc(memory_need(num_rows, num_cols));
+    }
+
+
+#if 0
+    // Move assignment (emulation)
+    self& operator=(self src)
+    {
+	// Self-copy would be an indication of an error
+	assert(this != &src);
+
+	this->check_dim(src.num_rows(), src.num_cols());
+	if (this->category == memory_base::view || src.category == memory_base::view)
+	    matrix_copy(src, *this);
+	else
+	    memory_base::move_assignment(src);
+	return *this;
+    }
+#endif
+
+
+#ifdef MTL_WITH_MOVE
+    /// Move Assignment
+    self& operator=(self&& 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
+	    memory_base::move_assignment(src);
+	std::cout << "In dense2D::move_assignment\n";
+	return *this;
+    }
+
+    /// Copy Assignment
+    self& operator=(const self& src)
+    {
+	this->checked_change_dim(src.num_rows(), src.num_cols());
+	// this->check_dim(src.num_rows(), src.num_cols());
+	matrix_copy(src, *this);
+	return *this;
+    }
+#else
+    /// Copy assignment (with move emulation)
+    self& operator=(self src)
+    {
+	// Self-copy would be an indication of an error
+	assert(this != &src);
+
+	this->check_dim(src.num_rows(), src.num_cols());
+	if (this->category == memory_base::view || src.category == memory_base::view)
+	    matrix_copy(src, *this);
+	else
+	    memory_base::move_assignment(src);
+	return *this;
+    }
+#endif
+
+    using assign_base::operator=;
+
+
+    value_type operator() (key_type const& key) const
+    {
+	return this->data[key.dilated_row.dilated_value() + key.dilated_col.dilated_value()];
+    }
+
+    void operator()(key_type const& key, value_type const& value)
+    {
+	this->data[key.dilated_row.dilated_value() + key.dilated_col.dilated_value()]= value;
+    }
+
+    /// Constant reference to A[row][col]
+    const_reference operator() (size_type row, size_type col) const
+    {
+	MTL_DEBUG_THROW_IF(is_negative(row) || row >= this->num_rows() || is_negative(col) || col >= this->num_cols(), index_out_of_range());
+	return this->data[dilated_row_t(row).dilated_value() + dilated_col_t(col).dilated_value()];
+    }
+
+    /// Mutable reference to A[row][col]
+    value_type& operator() (size_type row, size_type col)
+    {
+	MTL_DEBUG_THROW_IF(is_negative(row) || row >= this->num_rows() || is_negative(col) || col >= this->num_cols(), index_out_of_range());
+	return this->data[dilated_row_t(row).dilated_value() + dilated_col_t(col).dilated_value()];
+    }
+
+    void crop() {} ///< Delete structural zeros, only dummy here
+
+  protected:
+    void set_nnz()
+    {
+      this->my_nnz = this->num_rows() * this->num_cols();
+    }
+
+    size_type memory_need(size_type rows, size_type cols)
+    {
+        dilated_row_t n_rows(rows - 1);
+        dilated_col_t n_cols(cols - 1);
+        return (n_rows.dilated_value() + n_cols.dilated_value() + 1);
+    }
+
+    /// Swap 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));
+    }
+
+    template <typename> friend struct sub_matrix_t;
+};
+
+
+// ================
+// Free functions
+// ================
+
+/// Number of rows
+template <typename Value, std::size_t Mask, typename Parameters>
+typename morton_dense<Value, Mask, Parameters>::size_type
+inline num_rows(const morton_dense<Value, Mask, Parameters>& matrix)
+{
+    return matrix.num_rows();
+}
+
+/// Number of columns
+template <typename Value, std::size_t Mask, typename Parameters>
+typename morton_dense<Value, Mask, Parameters>::size_type
+inline num_cols(const morton_dense<Value, Mask, Parameters>& matrix)
+{
+    return matrix.num_cols();
+}
+
+/// Matrix size, i.e. number of rows times columns
+template <typename Value, std::size_t Mask, typename Parameters>
+typename morton_dense<Value, Mask, Parameters>::size_type
+inline size(const morton_dense<Value, Mask, Parameters>& matrix)
+{
+    return matrix.num_cols() * matrix.num_rows();
+}
+
+}} // namespace mtl::matrix
+
+
+
+
+// ================
+// Range generators
+// ================
+
+namespace mtl { namespace traits {
+
+    // VC 8.0 finds ambiguity with mtl::tag::morton_dense (I wonder why)
+    using mtl::matrix::morton_dense;
+    using mtl::matrix::morton_dense_el_cursor;
+    using mtl::matrix::morton_dense_col_cursor;
+    using mtl::matrix::morton_dense_row_cursor;
+    using mtl::matrix::morton_dense_col_const_iterator;
+    using mtl::matrix::morton_dense_row_const_iterator;
+    using mtl::matrix::morton_dense_col_iterator;
+    using mtl::matrix::morton_dense_row_iterator;
+
+    // ===========
+    // For cursors
+    // ===========
+
+    template <class Elt, std::size_t BitMask, class Parameters>
+    struct range_generator<glas::tag::all, morton_dense<Elt, BitMask, Parameters> >
+    {
+	typedef morton_dense<Elt, BitMask, Parameters>        Matrix;
+	typedef complexity_classes::linear_cached             complexity;
+	static int const                                      level = 1;
+	typedef morton_dense_el_cursor<BitMask>  type;
+	type begin(Matrix const& matrix)
+	{
+	    return type(matrix.begin_row(), matrix.begin_col(), matrix.num_cols());
+	}
+	type end(Matrix const& matrix)
+	{
+	    return type(matrix.end_row(), matrix.begin_col(), matrix.num_cols());
+	}
+    };
+
+    template <class Elt, std::size_t BitMask, class Parameters>
+    struct range_generator<glas::tag::nz, morton_dense<Elt, BitMask, Parameters> >
+	: range_generator<glas::tag::all, morton_dense<Elt, BitMask, Parameters> >
+    {};
+
+    template <class Elt, std::size_t BitMask, class Parameters>
+    struct range_generator<glas::tag::row, morton_dense<Elt, BitMask, Parameters> >
+	: detail::all_rows_range_generator<morton_dense<Elt, BitMask, Parameters>, complexity_classes::linear_cached>
+    {};
+
+    // For a cursor pointing to some row give the range of elements in this row
+    template <class Elt, std::size_t BitMask, class Parameters>
+    struct range_generator<glas::tag::nz,
+			   detail::sub_matrix_cursor<morton_dense<Elt, BitMask, Parameters>, glas::tag::row, 2> >
+    {
+	typedef morton_dense<Elt, BitMask, Parameters>                   matrix;
+	typedef typename Collection<matrix>::size_type                   size_type;
+	typedef detail::sub_matrix_cursor<matrix, glas::tag::row, 2>     cursor;
+	typedef complexity_classes::linear_cached                        complexity;
+	static int const                                                 level = 1;
+	typedef morton_dense_col_cursor<BitMask>                         type;
+
+	type begin(cursor const& c) const
+	{
+	    return type(c.key, c.ref.begin_col());
+	}
+	type end(cursor const& c) const
+	{
+	    return type(c.key, c.ref.end_col());
+	}
+	type lower_bound(cursor const& c, size_type position) const
+	{
+	    return type(c.key, std::min(c.ref.end_col(), position));
+	}
+    };
+
+    template <class Elt, std::size_t BitMask, class Parameters>
+    struct range_generator<glas::tag::all,
+			   detail::sub_matrix_cursor<morton_dense<Elt, BitMask, Parameters>, glas::tag::row, 2> >
+        : range_generator<glas::tag::nz,
+			  detail::sub_matrix_cursor<morton_dense<Elt, BitMask, Parameters>, glas::tag::row, 2> >
+    {};
+
+    template <class Elt, std::size_t BitMask, class Parameters>
+    struct range_generator<glas::tag::col, morton_dense<Elt, BitMask, Parameters> >
+	: detail::all_cols_range_generator<morton_dense<Elt, BitMask, Parameters>, complexity_classes::linear_cached>
+    {};
+
+    // For a cursor pointing to some row give the range of elements in this row
+    template <class Elt, std::size_t BitMask, class Parameters>
+    struct range_generator<glas::tag::nz,
+			   detail::sub_matrix_cursor<morton_dense<Elt, BitMask, Parameters>, glas::tag::col, 2> >
+    {
+	typedef morton_dense<Elt, BitMask, Parameters>                   matrix;
+	typedef typename Collection<matrix>::size_type                   size_type;
+	typedef detail::sub_matrix_cursor<matrix, glas::tag::col, 2>     cursor;
+	typedef complexity_classes::linear_cached                        complexity;
+	static int const                                                 level = 1;
+	typedef morton_dense_row_cursor<BitMask>                         type;
+
+	type begin(cursor const& c)
+	{
+	    return type(c.ref.begin_row(), c.key);
+	}
+	type end(cursor const& c)
+	{
+	    return type(c.ref.end_row(), c.key);
+	}
+	type lower_bound(cursor const& c, size_type position) const
+	{
+	    return type(std::min(c.ref.end_row(), position), c.key);
+	}
+    };
+
+    template <class Elt, std::size_t BitMask, class Parameters>
+    struct range_generator<glas::tag::all,
+			   detail::sub_matrix_cursor<morton_dense<Elt, BitMask, Parameters>, glas::tag::col, 2> >
+        : range_generator<glas::tag::nz,
+			  detail::sub_matrix_cursor<morton_dense<Elt, BitMask, Parameters>, glas::tag::col, 2> >
+    {};
+
+
+// =============
+// For iterators
+// =============
+
+    namespace detail {
+
+        template <typename OuterTag, typename Matrix, bool is_const>
+        struct morton_dense_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;
+
+	    typedef complexity_classes::linear_cached                                     complexity;
+	    static int const                                                              level = 1;
+
+	    typedef typename boost::mpl::if_<
+		boost::is_same<OuterTag, glas::tag::row>
+	      , typename boost::mpl::if_c<
+    	            is_const
+		  , morton_dense_col_const_iterator<Matrix>
+		  , morton_dense_col_iterator<Matrix>
+		>::type
+	      , typename boost::mpl::if_c<
+    	            is_const
+		  , morton_dense_row_const_iterator<Matrix>
+		  , morton_dense_row_iterator<Matrix>
+		>::type
+	    >::type type;
+
+        private:
+
+	    typedef typename boost::mpl::if_c<is_const, const Matrix&, Matrix&>::type    mref_type;
+
+	    type begin_dispatch(cursor const& c, glas::tag::row)
+	    {
+		return type(const_cast<mref_type>(c.ref), c.key, c.ref.begin_col());
+	    }
+
+	    type end_dispatch(cursor const& c, glas::tag::row)
+	    {
+		return type(const_cast<mref_type>(c.ref), c.key, c.ref.end_col());
+	    }
+
+	    type begin_dispatch(cursor const& c, glas::tag::col)
+	    {
+		return type(const_cast<mref_type>(c.ref), c.ref.begin_row(), c.key);
+	    }
+
+	    type end_dispatch(cursor const& c, glas::tag::col)
+	    {
+		return type(const_cast<mref_type>(c.ref), c.ref.end_row(), c.key);
+	    }
+
+        public:
+
+	    type begin(cursor const& c)
+	    {
+		return begin_dispatch(c, OuterTag());
+	    }
+
+	    type end(cursor const& c)
+	    {
+		return end_dispatch(c, OuterTag());
+	    }
+        };
+
+    } // namespace detail
+
+
+    template <typename Value, std::size_t BitMask, typename Parameters, typename OuterTag>
+    struct range_generator<tag::iter::nz,
+			   detail::sub_matrix_cursor<morton_dense<Value, BitMask, Parameters>, OuterTag, 2> >
+      : public detail::morton_dense_iterator_range_generator<OuterTag, morton_dense<Value, BitMask, Parameters>, false>
+    {};
+
+    template <typename Value, std::size_t BitMask, typename Parameters, typename OuterTag>
+    struct range_generator<tag::iter::all,
+			   detail::sub_matrix_cursor<morton_dense<Value, BitMask, Parameters>, OuterTag, 2> >
+      : public detail::morton_dense_iterator_range_generator<OuterTag, morton_dense<Value, BitMask, Parameters>, false>
+    {};
+
+    template <typename Value, std::size_t BitMask, typename Parameters, typename OuterTag>
+    struct range_generator<tag::const_iter::nz,
+			   detail::sub_matrix_cursor<morton_dense<Value, BitMask, Parameters>, OuterTag, 2> >
+      : public detail::morton_dense_iterator_range_generator<OuterTag, morton_dense<Value, BitMask, Parameters>, true>
+    {};
+
+    template <typename Value, std::size_t BitMask, typename Parameters, typename OuterTag>
+    struct range_generator<tag::const_iter::all,
+			   detail::sub_matrix_cursor<morton_dense<Value, BitMask, Parameters>, OuterTag, 2> >
+      : public detail::morton_dense_iterator_range_generator<OuterTag, morton_dense<Value, BitMask, Parameters>, true>
+    {};
+
+
+}} // namespace mtl::traits
+
+
+namespace mtl { namespace matrix {
+
+    // ==========
+    // Sub matrix
+    // ==========
+
+    template <typename Value, std::size_t BitMask, typename Parameters>
+    struct sub_matrix_t<morton_dense<Value, BitMask, Parameters> >
+    {
+        typedef morton_dense<Value, BitMask, Parameters>    matrix_type;
+        typedef matrix_type                     sub_matrix_type;
+        typedef 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, morton_dense_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 {
+
+	using matrix::morton_dense;
+
+    // Enable cloning of dense matrices
+    template <typename Value, std::size_t BitMask, typename Parameters>
+    struct is_clonable< matrix::morton_dense<Value, BitMask, Parameters> > : boost::mpl::true_ {};
+
+} // namespace mtl
+
+#endif // MTL_MORTON_DENSE_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/multi_vector.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_MULTI_VECTOR_INCLUDE
+#define MTL_MATRIX_MULTI_VECTOR_INCLUDE
+
+#include <cassert>
+#include <boost/utility/enable_if.hpp>
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/matrix/parameter.hpp>
+#include <boost/numeric/mtl/matrix/crtp_base_matrix.hpp>
+#include <boost/numeric/mtl/matrix/mat_expr.hpp>
+#include <boost/numeric/mtl/matrix/multi_vector_range.hpp>
+#include <boost/numeric/mtl/matrix/make_fast_multi_vector_expr.hpp>
+#include <boost/numeric/mtl/utility/is_what.hpp>
+#include <boost/numeric/mtl/utility/is_composable_vector.hpp>
+#include <boost/numeric/mtl/utility/is_multi_vector_expr.hpp>
+#include <boost/numeric/mtl/utility/fast_multi_vector_expr.hpp>
+#include <boost/numeric/mtl/utility/static_assert.hpp>
+#include <boost/numeric/mtl/vector/parameter.hpp>
+
+
+namespace mtl { namespace matrix {
+
+
+// Might need to be defined later
+struct multi_vector_key {};
+
+/// Matrix constituting of set of column vectors (under development)
+template <typename Vector>
+class multi_vector
+  : public base_matrix<typename mtl::Collection<Vector>::value_type, parameters<> >,
+    public crtp_base_matrix< multi_vector<Vector>, typename Collection<Vector>::value_type,
+			     typename Collection<Vector>::size_type>,
+    public mat_expr< multi_vector<Vector> >
+{
+    typedef base_matrix<typename Collection<Vector>::value_type, parameters<> >           super;
+
+    // Vector must by column vector
+    MTL_STATIC_ASSERT((boost::is_same<typename OrientedCollection<Vector>::orientation,
+			              tag::col_major>::value),
+		      "Vector must be a column vector.");
+  public:
+    typedef multi_vector                             self;
+    // typedef mtl::matrix::parameters<>                parameters;
+    typedef Vector                                   vector_type;
+    typedef tag::col_major                           orientation;
+    typedef typename Collection<Vector>::value_type  value_type;
+    typedef typename Collection<Vector>::size_type   size_type;
+    typedef const value_type&                        const_reference;
+    typedef value_type&                              reference;
+    typedef multi_vector_key                         key_type;
+    typedef crtp_matrix_assign< self, value_type, size_type >    assign_base;
+
+    multi_vector() : super(non_fixed::dimensions(0, 0)), master(0) {}
+
+  private:
+    void setup_data(size_type num_rows, size_type num_cols, boost::mpl::false_)
+    {
+	for (size_type i= 0; i < num_cols; ++i)
+	    data[i]= Vector(num_rows);
+	master= 0;
+    }
+
+    void setup_data(size_type num_rows, size_type num_cols, boost::mpl::true_)
+    {
+	master= new Vector(num_rows * num_cols);
+	for (size_type i= 0; i < num_cols; ++i) {
+	    Vector tmp((*master)[irange(i * num_rows, (i+1) * num_rows)]);
+	    swap(data[i], tmp);
+	}
+    }
+
+  public:
+    /// Constructor by number of rows and columns
+    multi_vector(size_type num_rows, size_type num_cols)
+      : super(non_fixed::dimensions(num_rows, num_cols)), data(num_cols)
+    {
+	setup_data(num_rows, num_cols, mtl::traits::is_composable_vector<Vector>());
+	this->my_nnz= num_rows * num_cols;
+    }
+
+    /// Constructor column vector and number of columns (for easier initialization)
+    multi_vector(const Vector& v, size_type num_cols)
+      : super(non_fixed::dimensions(size(v), num_cols)), data(num_cols)
+    {
+	using mtl::vector::num_rows;
+	setup_data(num_rows(v), num_cols, mtl::traits::is_composable_vector<Vector>());
+	for (size_type i= 0; i < num_cols; ++i)
+	    data[i]= v;
+	this->my_nnz= num_cols * size(v);
+    }
+
+    ~multi_vector() { delete master; }
+
+  private:
+
+    void change_dim(size_type r, size_type c, boost::mpl::false_)
+    {
+	for (size_type i= 0; i < c; i++)
+	    data[i].change_dim(r);
+    }
+
+    void change_dim(size_type r, size_type c, boost::mpl::true_)
+    {
+	delete master;
+	setup_data(r, c, boost::mpl::true_());
+    }
+
+  public:
+    /// Change dimension, can keep old data
+    void change_dim(size_type r, size_type c)
+    {
+	if (r == super::num_rows() && c == super::num_cols()) return;
+	super::change_dim(r, c);
+	data.change_dim(c);
+	change_dim(r, c, mtl::traits::is_composable_vector<Vector>());
+    }
+
+    void self_assignment(const self& src, boost::mpl::false_)
+    {
+	for (size_type i= 0; i < super::num_cols(); i++)
+	    data[i]= src.data[i];
+    }
+
+    void self_assignment(const self& src, boost::mpl::true_)
+    {	*master= *src.master;    }
+
+    /// Copy constructor
+    /** Explicitly needed now. **/
+    self& operator=(const self& src)
+    {
+	assign_base::checked_change_dim(src.num_rows(), src.num_cols());
+	self_assignment(src, mtl::traits::is_composable_vector<Vector>());
+	return *this;
+    }
+
+    // Todo: multi_vector with other matrix expressions
+    /// Assign multi_vector and expressions thereof, general matrices currently not allowed
+    template <typename Src>
+    typename boost::enable_if_c<mtl::traits::is_multi_vector_expr<Src>::value
+				&& !mtl::traits::is_fast_multi_vector_expr<Src>::value, self&>::type
+    operator=(const Src& src)
+    {
+	MTL_THROW_IF((mtl::matrix::num_rows(src) != super::num_rows()
+		      || mtl::matrix::num_cols(src) != super::num_cols()), incompatible_size());
+	for (std::size_t i= 0, n= super::num_cols(); i < n; ++i)
+	    vector(i)= src.vector(i);
+	return *this;
+    }
+
+    /// Assign multi_vector and expressions thereof that are stored internally by a single vector
+    template <typename Src>
+    typename boost::enable_if<mtl::traits::is_fast_multi_vector_expr<Src>, self&>::type
+    operator=(const Src& src)
+    {
+	assert(master);
+	*master= make_fast_multi_vector_expr(src);
+	return *this;
+    }
+
+    template <typename Src>
+    typename boost::enable_if_c<mtl::traits::is_matrix<Src>::value
+				&& !mtl::traits::is_multi_vector_expr<Src>::value, self&>::type
+    operator=(const Src& src)
+    {
+	assign_base::operator=(src);
+	return *this;
+    }
+
+    /// Assign scalar
+    template <typename Src>
+    typename boost::enable_if<mtl::traits::is_scalar<Src>, self&>::type
+    operator=(const Src& src)
+    {
+	assign_base::operator=(src);
+	return *this;
+    }
+
+    const_reference operator() (size_type i, size_type j) const { return data[j][i]; }
+    reference operator() (size_type i, size_type j) { return data[j][i]; }
+
+    Vector& vector(size_type i) { return data[i]; }
+    const Vector& vector(size_type i) const { return data[i]; }
+
+    Vector& at(size_type i) { return data[i]; }
+    const Vector& at(size_type i) const { return data[i]; }
+
+    multi_vector_range<Vector> vector(irange const& r) const { return multi_vector_range<Vector>(*this, r); }
+
+    inline friend const Vector& make_fast_multi_vector_expr(const self& mv) { assert(mv.master); return *mv.master; }
+    inline friend Vector& make_fast_multi_vector_expr(self& mv) { assert(mv.master); return *mv.master; }
+
+  protected:
+    mtl::vector::dense_vector<Vector, mtl::vector::parameters<> >          data;
+    Vector*                                                                master;
+};
+
+/// Number of rows
+template< typename Vector >
+typename Collection< Vector >::size_type num_cols(const multi_vector< Vector >& A) { return A.num_cols(); }
+
+/// Number of columns
+template< typename Vector >
+typename Collection< Vector >::size_type num_rows(const multi_vector< Vector >& A) { return A.num_rows(); }
+
+/// Size as defined by number of rows times columns
+template< typename Vector >
+typename Collection< Vector >::size_type size(const multi_vector< Vector >& A) { return num_rows(A) * num_cols(A); }
+}} // namespace mtl::matrix
+
+namespace mtl {
+	using matrix::multi_vector;
+}
+
+#endif // MTL_MATRIX_MULTI_VECTOR_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/multi_vector_range.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_MULTI_VECTOR_RANGE_INCLUDE
+#define MTL_MATRIX_MULTI_VECTOR_RANGE_INCLUDE
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/matrix/multi_vector.hpp>
+
+namespace mtl { namespace matrix {
+
+// So far only const, might be refactored for mutability later
+template <typename Vector>
+class multi_vector_range
+{
+  public:
+    typedef multi_vector_range                       self;
+    typedef typename Collection<Vector>::size_type   size_type;
+    typedef typename Collection<Vector>::value_type  value_type;
+    typedef const value_type&                        const_reference;
+
+    multi_vector_range(const multi_vector<Vector>& ref, const irange& r) : ref(ref), r(r) {}
+
+    const_reference operator() (size_type i, size_type j) const { return ref[r.to_range(j)][i]; }
+    const Vector& vector(size_type i) const { return ref.vector(r.to_range(i)); }
+
+    /// Number of rows
+    friend size_type num_rows(const self& A) { return num_rows(A.ref); }
+
+    /// Number of columns
+    friend size_type num_cols(const self& A) { return A.r.size(); }
+
+    /// Size as defined by number of rows times columns
+    friend size_type size(const self& A) { return num_rows(A) * num_cols(A); }
+
+  protected:
+    const multi_vector<Vector>& ref;
+    const irange                r;
+};
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_MULTI_VECTOR_RANGE_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/operators.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_OPERATORS_INCLUDE
+#define MTL_MATRIX_OPERATORS_INCLUDE
+
+#include <boost/utility/enable_if.hpp>
+#include <boost/mpl/and.hpp>
+
+#include <boost/numeric/mtl/utility/ashape.hpp>
+#include <boost/numeric/mtl/utility/is_multi_vector_expr.hpp>
+#include <boost/numeric/mtl/utility/static_assert.hpp>
+#include <boost/numeric/mtl/matrix/all_mat_expr.hpp>
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+
+
+namespace mtl { namespace matrix {
+
+template <typename E1, typename E2>
+inline mat_mat_plus_expr<E1, E2>
+operator+ (const mat_expr<E1>& e1, const mat_expr<E2>& e2)
+{
+    // do not add matrices with inconsistent value types
+    MTL_STATIC_ASSERT((boost::is_same<typename ashape::ashape<E1>::type,
+				      typename ashape::ashape<E2>::type>::value), "Matrices have not consistent algebraic shape (i.e. nested types).");
+    return mat_mat_plus_expr<E1, E2>(static_cast<const E1&>(e1), static_cast<const E2&>(e2));
+}
+
+// if enabled it has priority over previous functions because that performs upcast
+template <typename E1, typename E2>
+typename boost::enable_if_c<mtl::traits::is_multi_vector_expr<E1>::value && mtl::traits::is_multi_vector_expr<E2>::value, mv_mv_plus_expr<E1, E2> >::type
+inline operator+(const E1& e1, const E2& e2)
+{
+    return mv_mv_plus_expr<E1, E2>(e1, e2);
+}
+
+// Specialization for multi_vector
+// template <typename V1, typename V2>
+// inline mv_mv_plus_expr<multi_vector<V1>, multi_vector<V2> >
+// operator+(const multi_vector<V1>& m1, const multi_vector<V2>& m2)
+// {
+//     return mv_mv_plus_expr<multi_vector<V1>, multi_vector<V2> >(m1, m2);
+// }
+
+#if 0
+// Planned for future optimizations on sums of dense matrix expressions
+template <typename E1, typename E2>
+inline dmat_dmat_plus_expr<E1, E2>
+operator+ (const dmat_expr<E1>& e1, const dmat_expr<E2>& e2)
+{
+    // do not add matrices with inconsistent value types
+    MTL_STATIC_ASSERT((boost::is_same<typename ashape::ashape<E1>::type,
+				      typename ashape::ashape<E2>::type>::value), "Matrices have not consistent algebraic shape (i.e. nested types).");
+    return dmat_dmat_plus_expr<E1, E2>(static_cast<const E1&>(e1), static_cast<const E2&>(e2));
+}
+#endif
+
+
+template <typename E1, typename E2>
+inline mat_mat_minus_expr<E1, E2>
+operator- (const mat_expr<E1>& e1, const mat_expr<E2>& e2)
+{
+    // do not add matrices with inconsistent value types
+    MTL_STATIC_ASSERT((boost::is_same<typename ashape::ashape<E1>::type,
+				      typename ashape::ashape<E2>::type>::value), "Matrices have not consistent algebraic shape (i.e. nested types).");
+    return mat_mat_minus_expr<E1, E2>(static_cast<const E1&>(e1), static_cast<const E2&>(e2));
+}
+
+// if enabled it has priority over previous functions because that performs upcast
+template <typename E1, typename E2>
+typename boost::enable_if_c<mtl::traits::is_multi_vector_expr<E1>::value && mtl::traits::is_multi_vector_expr<E2>::value, mv_mv_minus_expr<E1, E2> >::type
+inline operator-(const E1& e1, const E2& e2)
+{
+    return mv_mv_minus_expr<E1, E2>(e1, e2);
+}
+
+template <typename E1, typename E2>
+inline mat_mat_ele_times_expr<E1, E2>
+ele_prod(const mat_expr<E1>& e1, const mat_expr<E2>& e2)
+{
+    // do not multiply matrices element-wise with inconsistent value types
+    MTL_STATIC_ASSERT((boost::is_same<typename ashape::ashape<E1>::type,
+				      typename ashape::ashape<E2>::type>::value), "Matrices do not have consistent algebraic shape (i.e. nested types).");
+    return mat_mat_ele_times_expr<E1, E2>(static_cast<const E1&>(e1), static_cast<const E2&>(e2));
+}
+
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_OPERATORS_INCLUDE