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

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/parameter.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_PARAMETERS_INCLUDE
+#define MTL_MATRIX_PARAMETERS_INCLUDE
+
+#include <boost/numeric/mtl/utility/static_assert.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/detail/index.hpp>
+#include <boost/numeric/mtl/matrix/dimension.hpp>
+#include <boost/numeric/mtl/utility/is_static.hpp>
+
+namespace mtl { namespace matrix {
+
+/// Type for bundling template parameters of common matrix types
+/** OnStack = true can only be used with fixed::dimensions.
+    \sa \ref matrix_parameters
+    \sa \ref tuning_fsize
+    \sa \ref tuning_sizetype **/
+template <typename Orientation= row_major,
+	  typename Index= index::c_index,
+	  typename Dimensions= mtl::non_fixed::dimensions,
+	  bool OnStack= mtl::traits::is_static<Dimensions>::value,
+	  typename SizeType= std::size_t>
+struct parameters
+{
+    typedef Orientation orientation;
+    typedef Index       index;
+    typedef Dimensions  dimensions;
+    static bool const   on_stack= OnStack;
+    typedef SizeType    size_type;
+
+    // Matrix dimensions must be known at compile time to be on the stack
+    MTL_STATIC_ASSERT(( !on_stack || dimensions::is_static ), "Types to be stored on stack must provide static size.");
+};
+
+/// Short-cut to define parameters with unsigned and defaults otherwise
+typedef parameters<row_major, index::c_index, mtl::non_fixed::dimensions, false, unsigned> unsigned_parameters;
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_PARAMETERS_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/permutation.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_PERMUTATION_INCLUDE
+#define MTL_MATRIX_PERMUTATION_INCLUDE
+
+#include <boost/numeric/mtl/operation/size.hpp>
+#include <boost/numeric/mtl/matrix/reorder.hpp>
+
+namespace mtl { namespace matrix {
+
+
+namespace traits {
+
+    //\ Return type of mtl::matrix::permutation
+    // Only for completeness
+    template <typename Value= short>
+    struct permutation
+    {
+	typedef typename reorder<Value>::type  type;
+    };
+}
+
+template <typename Value, typename PermutationVector>
+typename traits::permutation<Value>::type
+inline permutation(const PermutationVector& v)
+{
+    using mtl::size;
+    return reorder(v, size(v));
+}
+
+/// Computes permutation matrix from corresponding vector
+template <typename PermutationVector>
+typename traits::permutation<>::type
+inline permutation(const PermutationVector& v)
+{
+    return permutation<short>(v);
+}
+
+
+}} // namespace mtl::matrix
+
+namespace mtl { namespace vector {
+
+    /// Import into vector namespace; see \ref mtl::matrix::permutation
+    using mtl::matrix::permutation;
+
+}} // namespace mtl::vector
+
+#endif // MTL_MATRIX_PERMUTATION_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/poisson2D_dirichlet.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_POISSON2D_DIRICHLET_INCLUDE
+#define MTL_MATRIX_POISSON2D_DIRICHLET_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/ashape.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/vector/mat_cvec_multiplier.hpp>
+
+namespace mtl { namespace matrix {
+
+/// Matrix-free linear operator for a Poisson equation on a rectangular domain of \p m by \p n with Dirichlet boundary conditions
+struct poisson2D_dirichlet
+{
+    /// Constructor
+    poisson2D_dirichlet(int m, int n) : m(m), n(n), s(m * n) {}
+
+    /// Member function that realizes the multiplication
+    template <typename VectorIn, typename VectorOut, typename Assign>
+    void mult(const VectorIn& v, VectorOut& w, Assign) const
+    {
+	MTL_DEBUG_THROW_IF(int(size(v)) != s, incompatible_size());
+	MTL_DEBUG_THROW_IF(size(v) != size(w), incompatible_size());
+
+	const int nb = n < 3 ? 1 : (n - 2) / 4 * 4 + 1;
+
+	// Inner domain
+	for (int i= 1; i < m-1; i++) {
+	    int kmax= i * n + nb;
+	    for (int k= i * n + 1; k < kmax; k+= 4) {
+		typename Collection<VectorIn>::value_type const v0= v[k], v1= v[k+1], v2= v[k+2], v3= v[k+3];
+		Assign::apply(w[k], 4 * v0 - v[k-n] - v[k+n] - v[k-1] - v1);
+		Assign::apply(w[k+1], 4 * v1 - v[k-n+1] - v[k+n+1] - v0 - v2);
+		Assign::apply(w[k+2], 4 * v2 - v[k-n+2] - v[k+n+2] - v1 - v3);
+		Assign::apply(w[k+3], 4 * v3 - v[k-n+3] - v[k+n+3] - v2 - v[k+4]);
+	    }
+	    for (int j= nb, k= i * n + j; j < n-1; j++, k++)
+		Assign::apply(w[k], 4 * v[k] - v[k-n] - v[k+n] - v[k-1] - v[k+1]);
+	}
+
+	// Upper border
+	for (int j= 1; j < n-1; j++)
+	    Assign::apply(w[j], 4 * v[j] - v[j+n] - v[j-1] - v[j+1]);
+
+	// Lower border
+	for (int j= 1, k= (m-1) * n + j; j < n-1; j++, k++)
+	    Assign::apply(w[k], 4 * v[k] - v[k-n] - v[k-1] - v[k+1]);
+
+	// Left border
+	for (int i= 1, k= n; i < m-1; i++, k+= n)
+	    Assign::apply(w[k], 4 * v[k] - v[k-n] - v[k+n] - v[k+1]);
+
+	// Right border
+	for (int i= 1, k= n+n-1; i < m-1; i++, k+= n)
+	    Assign::apply(w[k], 4 * v[k] - v[k-n] - v[k+n] - v[k-1]);
+
+	// Corners
+	Assign::apply(w[0], 4 * v[0] - v[1] - v[n]);
+	Assign::apply(w[n-1], 4 * v[n-1] - v[n-2] - v[2*n - 1]);
+	Assign::apply(w[(m-1)*n], 4 * v[(m-1)*n] - v[(m-2)*n] - v[(m-1)*n+1]);
+	Assign::apply(w[m*n-1], 4 * v[m*n-1] - v[m*n-2] - v[m*n-n-1]);
+    }
+
+    /// Multiplication is procastinated until we know where the product goes
+    template <typename VectorIn>
+    vector::mat_cvec_multiplier<poisson2D_dirichlet, VectorIn> operator*(const VectorIn& v) const
+    {	return vector::mat_cvec_multiplier<poisson2D_dirichlet, VectorIn>(*this, v);    }
+
+    int m, n, s;
+};
+
+inline std::size_t size(const poisson2D_dirichlet& A) { return A.s * A.s; } ///< Matrix size
+inline std::size_t num_rows(const poisson2D_dirichlet& A) { return A.s; } ///< Number of rows
+inline std::size_t num_cols(const poisson2D_dirichlet& A) { return A.s; } ///< Number of columns
+
+}} // namespace mtl::matrix
+
+namespace mtl {
+
+    template <>
+    struct Collection<matrix::poisson2D_dirichlet>
+    {
+	typedef double value_type;
+	typedef int    size_type;
+    };
+
+    namespace ashape {
+	template <> struct ashape_aux<matrix::poisson2D_dirichlet>
+	{	typedef nonscal type;    };
+    }
+}
+
+#endif // MTL_MATRIX_POISSON2D_DIRICHLET_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/reorder.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_REORDER_INCLUDE
+#define MTL_MATRIX_REORDER_INCLUDE
+
+#include <algorithm>
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/matrix/parameter.hpp>
+#include <boost/numeric/mtl/matrix/compressed2D.hpp>
+#include <boost/numeric/mtl/matrix/reorder_ref.hpp>
+
+namespace mtl { namespace matrix {
+
+
+namespace traits {
+
+    /// Return type of mtl::matrix::reorder
+    template <typename Value= short>
+    struct reorder
+    {
+	typedef mtl::matrix::compressed2D<Value, parameters<> >  type;
+    };
+}
+
+
+template <typename Value, typename ReorderVector>
+typename traits::reorder<Value>::type
+reorder(const ReorderVector& v, std::size_t cols= 0)
+{
+    typename traits::reorder<Value>::type A;
+    reorder_ref(v, A, cols);
+    return A;
+}
+
+
+/// Computes reordering matrix from corresponding vector
+template <typename ReorderVector>
+typename traits::reorder<>::type
+inline reorder(const ReorderVector& v, std::size_t cols= 0)
+{
+    return reorder<short>(v, cols);
+}
+
+
+}} // namespace mtl::matrix
+
+namespace mtl { namespace vector {
+
+    /// Import into vector namespace; see \ref mtl::matrix::reorder
+    using mtl::matrix::reorder;
+
+}} // namespace mtl::vector
+
+#endif // MTL_MATRIX_REORDER_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/reorder_matrix_rows.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_REORDER_MATRIX_ROWS_INCLUDE
+#define MTL_MATRIX_REORDER_MATRIX_ROWS_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/range_generator.hpp>
+#include <boost/numeric/mtl/utility/property_map.hpp>
+#include <boost/numeric/mtl/matrix/inserter.hpp>
+#include <boost/numeric/mtl/operation/size.hpp>
+
+namespace mtl { namespace matrix {
+
+/// Reorder the rows of a matrix without creating a reorder matrix
+/** This is less elegant but avoids cyclic dependencies. Does not work with CCS matrices. **/
+template <typename ReorderVector, typename Matrix>
+Matrix reorder_matrix_rows(const ReorderVector& v, const Matrix& A)
+{
+    using mtl::size;
+
+    typename mtl::traits::col<Matrix>::type                                   col(A);
+    typename mtl::traits::const_value<Matrix>::type                           value(A);
+    typedef typename mtl::traits::range_generator<tag::row, Matrix>::type     cursor_type;
+    typedef typename mtl::traits::range_generator<tag::nz, cursor_type>::type icursor_type;
+    typedef typename mtl::Collection<Matrix>::size_type                       size_type;
+
+    Matrix B(size(v), num_cols(A));
+
+    inserter<Matrix>      ins(B, size_type(B.nnz() / num_cols(B) * 1.2));
+
+    for (std::size_t i= 0; i < size(v); i++) {
+	cursor_type cursor(v[i], A);   // go to row given by reorder
+	for (icursor_type icursor = begin<tag::nz>(cursor), icend = end<tag::nz>(cursor); icursor != icend; ++icursor)
+	    ins[i][col(*icursor)] << value(*icursor);
+    }
+    return B;
+}
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_REORDER_MATRIX_ROWS_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/reorder_ref.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_REORDER_REF_INCLUDE
+#define MTL_MATRIX_REORDER_REF_INCLUDE
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/operation/size.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/matrix/inserter.hpp>
+#include <boost/numeric/linear_algebra/identity.hpp>
+
+namespace mtl { namespace matrix {
+
+
+template <typename ReorderVector, typename Matrix>
+void reorder_ref(const ReorderVector& v, Matrix& A, std::size_t cols= 0)
+{
+    using math::one; using mtl::size;
+    typedef typename Collection<Matrix>::value_type value_type;
+
+    if (mtl::size(v) == 0) {
+	A.change_dim(0, cols); return; }
+
+    // Find maximal entry (don't use mtl::max to allow for arrays and others)
+    std::size_t  s= mtl::size(v),
+	         my_max= std::size_t(*std::max_element(&v[0], &v[0] + s)) + 1;
+
+    if (cols == 0)
+	cols= my_max;
+    else
+	MTL_THROW_IF(my_max > cols, range_error("Too large value in reorder vector"));
+
+    A.change_dim(s, cols);
+    inserter<Matrix>      ins(A, 1);
+    for (std::size_t i= 0; i < s; i++)
+	ins[i][v[i]] << one(value_type());
+}
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_REORDER_REF_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/shifted_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_SHIFTED_INSERTER_INCLUDE
+#define MTL_MATRIX_SHIFTED_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>
+class shifted_inserter
+{
+  public:
+    typedef shifted_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 slot size and the offsets
+    shifted_inserter(matrix_type& A, size_type slot_size= 0,
+		     size_type row_offset= 0, size_type col_offset= 0)
+      : ins(A, slot_size), row_offset(row_offset), col_offset(col_offset) {}
+
+    void set_row_offset(size_type ro) { row_offset= ro; } ///< Change row offset
+    void set_col_offset(size_type co) { col_offset= co; } ///< Change column offset
+
+    size_type get_row_offset() const { return row_offset; } ///< Get row offset
+    size_type get_col_offset() const { return col_offset; } ///< Get column offset
+
+  private:
+    struct bracket_proxy
+    {
+	bracket_proxy(BaseInserter& ref, size_type row, size_type col_offset) : ref(ref), row(row), col_offset(col_offset) {}
+
+	proxy_type operator[](size_type col)
+	{   return proxy_type(ref, row, col+col_offset); }
+
+	BaseInserter&   ref;
+	size_type       row, col_offset;
+    };
+
+  public:
+    /// To be used in ins[r][c] << value;
+    bracket_proxy operator[] (size_type row)
+    {	return bracket_proxy(ins, row+row_offset, col_offset);    }
+
+    /// To be used in ins(r, c) << value;
+    proxy_type operator() (size_type row, size_type col)
+    {	return proxy_type(ins, row+row_offset, col+col_offset);    }
+
+    // update, modify and operator<< are used from BaseInserter
+
+  private:
+    BaseInserter ins;
+    size_type    row_offset, col_offset;
+};
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_SHIFTED_INSERTER_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/sparse_banded.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_SPARSE_BANDED_INCLUDE
+#define MTL_MATRIX_SPARSE_BANDED_INCLUDE
+
+#include <algorithm>
+#include <cassert>
+#include <ostream>
+#include <vector>
+#include <boost/type_traits/make_signed.hpp>
+
+#include <boost/numeric/mtl/matrix/dimension.hpp>
+#include <boost/numeric/mtl/matrix/parameter.hpp>
+#include <boost/numeric/mtl/matrix/base_matrix.hpp>
+#include <boost/numeric/mtl/matrix/crtp_base_matrix.hpp>
+#include <boost/numeric/mtl/matrix/mat_expr.hpp>
+#include <boost/numeric/mtl/operation/is_negative.hpp>
+#include <boost/numeric/mtl/operation/update.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/is_row_major.hpp>
+#include <boost/numeric/mtl/utility/maybe.hpp>
+#include <boost/numeric/mtl/utility/static_assert.hpp>
+
+namespace mtl { namespace matrix {
+
+/// Sparse banded matrix class
+template <typename Value, typename Parameters = matrix::parameters<> >
+class sparse_banded
+  : public base_matrix<Value, Parameters>,
+    public const_crtp_base_matrix< sparse_banded<Value, Parameters>, Value, typename Parameters::size_type >,
+    public crtp_matrix_assign< sparse_banded<Value, Parameters>, Value, typename Parameters::size_type >,
+    public mat_expr< sparse_banded<Value, Parameters> >
+{
+    MTL_STATIC_ASSERT((mtl::traits::is_row_major<Parameters>::value), "Only row-major sparse banded matrices supported so far.");
+
+    typedef std::size_t                                size_t;
+    typedef base_matrix<Value, Parameters>             super;
+    typedef sparse_banded<Value, Parameters>           self;
+
+  public:
+    typedef Value                                      value_type;
+    typedef typename Parameters::size_type             size_type;
+    typedef typename boost::make_signed<size_type>::type  band_size_type;
+    typedef crtp_matrix_assign<self, Value, size_type> assign_base;
+
+    /// Construct matrix of dimension \p nr by \p nc
+    explicit sparse_banded(size_type nr= 0, size_type nc= 0)
+      : super(non_fixed::dimensions(nr, nc)), data(0), inserting(false)
+    {}
+
+    /// Copy from other types
+    template <typename MatrixSrc>
+    explicit sparse_banded(const MatrixSrc& src) : data(0), inserting(false)
+    {	*this= src;    }
+
+
+    ~sparse_banded() { delete[] data; }
+    void check() const { MTL_DEBUG_THROW_IF(inserting, access_during_insertion()); }
+    void check(size_type MTL_DEBUG_ARG(r), size_type MTL_DEBUG_ARG(c)) const
+    {
+	check();
+	MTL_DEBUG_THROW_IF(is_negative(r) || r >= this->num_rows()
+			   || is_negative(c) || c >= this->num_cols(), index_out_of_range());
+    }
+
+    using assign_base::operator=;
+
+    void make_empty() ///< Delete all entries
+    {
+	delete[] data; data= 0;
+	bands.resize(0);
+    }
+
+    /// Change dimension to \p r by \p c; deletes all entries
+    void change_dim(size_type r, size_type c)
+    {
+	make_empty();
+	super::change_dim(r, c);
+    }
+
+    /// Offset of entry r,c if in matrix (otherwise offset of next entry)
+    utilities::maybe<size_type> offset(size_type r, size_type c) const
+    {
+	check(r, c);
+	band_size_type dia= band_size_type(c) - band_size_type(r);
+	size_type      b= find_dia(dia);
+	return utilities::maybe<size_type>(r * bands.size() + b, size_t(b) < bands.size() && bands[b] == dia);
+    }
+
+    /// Value of matrix entry
+    value_type operator()(size_type r, size_type c) const
+    {
+	using math::zero;
+	utilities::maybe<size_type> o= offset(r, c);
+	return o ? data[o.value()] : zero(value_type());
+    }
+
+    /// L-value reference of stored matrix entry
+    /** To be used with care; in debug mode, exception is thrown if entry is not found **/
+    value_type& lvalue(size_type r, size_type c)
+    {
+	utilities::maybe<size_type> o= offset(r, c);
+	MTL_DEBUG_THROW_IF(!o, runtime_error("Entry doesn't exist."));
+	return data[o.value()];
+    }
+
+    /// Print matrix on \p os
+    friend std::ostream& operator<<(std::ostream& os, const self& A)
+    {
+	const size_type bs= A.bands.size(), nc= num_cols(A), s= bs * num_rows(A);
+	print_size_max p;
+	for (size_type i= 0; i < s; i++)
+	    p(A.data[i]);
+
+	for (size_type r= 0; r < num_rows(A); r++) {
+	    os << '[';
+	    size_type b= 0;
+	    while (b < bs && -A.bands[b] > long(r)) b++;
+	    for (size_type c= 0; c < nc; c++) {
+		os.width(p.max);
+		if (b == bs || long(c) - long(r) < A.bands[b])
+		    os << Value(0);
+		else
+		    os << A.data[r * bs + b++];
+		os << (c + 1 == nc ? "]\n" : " ");
+	    }
+	}
+	return os;
+    }
+
+    /// Number of structural non-zeros (i.e. stored entries) which is the the number of bands times rows
+    size_type nnz() const { return bands.size() * this->num_rows(); }
+    size_type find_major(size_type offset) const { return offset % bands.size(); }
+
+    friend size_t num_rows(const self& A) { return A.num_rows(); }
+    friend size_t num_cols(const self& A) { return A.num_cols(); }
+
+    std::vector<band_size_type> const& ref_bands() const { return bands; } ///< Refer bands vector [advanced]
+    const value_type*                  ref_data() const { return data; } ///< Refer data pointer [advanced]
+
+  private:
+    size_type find_dia(band_size_type dia) const
+    {
+	size_type i= 0;
+	while (i < size_type(bands.size()) && dia > bands[i]) i++;
+	return i;
+    }
+
+    template <typename, typename, typename> friend struct sparse_banded_inserter;
+
+    std::vector<band_size_type>    bands;
+    value_type*                    data;
+    bool                           inserting;
+};
+
+/// Inserter for sparse banded matrix
+template <typename Value, typename Parameters, typename Updater = mtl::operations::update_store<Value> >
+struct sparse_banded_inserter
+{
+    typedef Value                                                            value_type;
+    typedef typename Parameters::size_type                                   size_type;
+    typedef sparse_banded<Value, Parameters>                                 matrix_type;
+    typedef sparse_banded_inserter                                           self;
+    typedef typename matrix_type::band_size_type                             band_size_type;
+    typedef operations::update_proxy<self, size_type>                        proxy_type;
+
+  private:
+    struct bracket_proxy
+    {
+	bracket_proxy(self& ref, size_type row) : ref(ref), row(row) {}
+
+	proxy_type operator[](size_type col) { return proxy_type(ref, row, col); }
+
+	self&      ref;
+	size_type  row;
+    };
+
+    void insert_dia(size_type dia_band, band_size_type dia)
+    {
+	using std::swap;
+	// empty vector and entry at the end
+	diagonals.push_back(std::vector<Value>(num_rows(A), Value(0)));
+	A.bands.push_back(dia);
+
+	for (size_type i= diagonals.size() - 1; i > dia_band; i--) {
+	    swap(diagonals[i-1], diagonals[i]);
+	    swap(A.bands[i-1], A.bands[i]);
+	}
+    }
+
+  public:
+    /// Construct inserter for matrix \p A; second argument for slot_size ignored
+    sparse_banded_inserter(matrix_type& A, size_type) : A(A)
+    {
+	MTL_THROW_IF(A.inserting, runtime_error("Two inserters on same matrix"));
+	MTL_THROW_IF(A.data, runtime_error("Inserter does not support modifications yet."));
+	A.inserting= true;
+    }
+
+    ~sparse_banded_inserter()
+    {
+	const size_type bs= A.bands.size();
+	Value* p= A.data= new Value[A.bands.size() * num_rows(A)];
+	for (size_type r= 0; r < num_rows(A); r++)
+	    for (size_type b= 0; b < bs; b++)
+		*p++= diagonals[b][r];
+	assert(p - A.data == long(A.bands.size() * num_rows(A)));
+	A.inserting= false;
+    }
+
+    /// Proxy to insert into A[row][col]
+    bracket_proxy operator[] (size_type row)
+    {
+	return bracket_proxy(*this, row);
+    }
+
+    /// Proxy to insert into A[row][col]
+    proxy_type operator() (size_type row, size_type col)
+    {
+	return proxy_type(*this, row, col);
+    }
+
+    /// Modify A[row][col] with \p val using \p Modifier
+    template <typename Modifier>
+    void modify(size_type row, size_type col, Value val)
+    {
+	MTL_DEBUG_THROW_IF(is_negative(row) || row >= num_rows(A) || is_negative(col) || col >= num_cols(A),
+			   index_out_of_range());
+	const band_size_type dia= col - row;
+	const size_type      dia_band= A.find_dia(dia);
+
+	if (dia_band == size_type(A.bands.size()) || dia != A.bands[dia_band])
+	    insert_dia(dia_band, dia);
+
+	Modifier()(diagonals[dia_band][row], val);
+    }
+
+    /// Modify A[row][col] with \p val using the class' updater
+    void update(size_type row, size_type col, Value val)
+    {
+	using math::zero;
+	modify<Updater>(row, col, val);
+    }
+
+  private:
+    matrix_type&                      A;
+    std::vector<std::vector<Value> >  diagonals;
+};
+
+
+template <typename Value, typename Parameters>
+struct sparse_banded_key
+{
+    typedef std::size_t                                   size_t;
+    typedef sparse_banded_key                             self;
+    typedef mtl::matrix::sparse_banded<Value, Parameters> matrix_type;
+
+    explicit sparse_banded_key(const matrix_type& A, size_t offset)
+      : A(A), offset(offset) {}
+
+    bool operator== (sparse_banded_key const& other) const { return offset == other.offset; }
+    bool operator!= (sparse_banded_key const& other) const { return !(*this == other); }
+
+    const matrix_type&   A;
+    size_t               offset;
+};
+
+// Cursor over every element
+template <typename Value, typename Parameters>
+struct sparse_banded_cursor
+  : public sparse_banded_key<Value, Parameters>
+{
+    typedef sparse_banded_key<Value, Parameters>          base;
+    typedef mtl::matrix::sparse_banded<Value, Parameters> matrix_type;
+    typedef sparse_banded_cursor                          self;
+
+    explicit sparse_banded_cursor(const matrix_type& A, size_t offset)
+      : base(A, offset)    {}
+
+    self& operator++() { ++this->offset; return *this;  }
+    const base& operator* () const { return *this; }
+};
+
+// ================
+// Free functions
+// ================
+
+template <typename Value, typename Parameters>
+inline std::size_t size(const sparse_banded<Value, Parameters>& matrix)
+{
+    return std::size_t(matrix.num_cols()) * std::size_t(matrix.num_rows());
+}
+
+}} // namespace mtl::matrix
+
+namespace mtl {
+    using matrix::sparse_banded;
+}
+
+
+// ================
+// Range generators
+// ================
+
+namespace mtl { namespace traits {
+
+    // Cursor over all rows
+    // Supported if row major matrix
+    template <typename Value, typename Parameters>
+    struct range_generator<glas::tag::row, sparse_banded<Value, Parameters> >
+      : detail::all_rows_range_generator<sparse_banded<Value, Parameters>, complexity_classes::linear_cached>
+    {};
+
+    template <class Value, class Parameters>
+    struct range_generator<glas::tag::nz,
+			   detail::sub_matrix_cursor<sparse_banded<Value, Parameters>, glas::tag::row, 2> >
+    {
+	typedef typename Parameters::size_type                                                 size_type;
+	typedef detail::sub_matrix_cursor<sparse_banded<Value, Parameters>, glas::tag::row, 2> cursor_type;
+	typedef complexity_classes::linear_cached                                              complexity;
+	typedef mtl::matrix::sparse_banded_cursor<Value, Parameters>                           type;
+	static int const                                                                       level = 1;
+	MTL_STATIC_ASSERT((is_row_major<Parameters>::value), "Only row-major sparse banded matrices supported so far.");
+
+	type begin(cursor_type const& cursor) const
+	{
+	    size_type r= cursor.key, b= 0, bs= cursor.ref.ref_bands().size();
+	    while (b < bs && -cursor.ref.ref_bands()[b] > long(r)) b++;
+	    return type(cursor.ref, r * bs + b);
+	}
+	type end(cursor_type const& cursor) const
+	{
+	    size_type r= cursor.key, bs= cursor.ref.ref_bands().size(), b= bs, nc= cursor.ref.num_cols();
+	    while (b > 0 && cursor.ref.ref_bands()[b-1] + long(r) >= long(nc)) b--;
+	    return type(cursor.ref, r * bs + b);
+	}
+
+	// type lower_bound(cursor_type const& cursor, size_type position) const {}
+    };
+
+
+    // Cursor over all rows or columns, depending which one is major
+    template <typename Value, typename Parameters>
+    struct range_generator<glas::tag::major, sparse_banded<Value, Parameters> >
+      : range_generator<glas::tag::row, sparse_banded<Value, Parameters> >
+    {};
+
+}} // namespace mtl::traits
+
+#endif // MTL_MATRIX_SPARSE_BANDED_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/strict_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_STRICT_LOWER_INCLUDE
+#define MTL_MATRIX_STRICT_LOWER_INCLUDE
+
+namespace mtl { namespace matrix {
+
+namespace traits {
+
+    template <typename Matrix>
+    struct strict_lower
+    {
+	typedef typename traits::bands<Matrix>::type type;
+    };
+}
+
+///  Strict lower triangular matrix
+template <typename Matrix>
+typename traits::strict_lower<Matrix>::type
+inline strict_lower(const Matrix& A)
+{
+    return bands(A, std::numeric_limits<long>::min(), 0);
+}
+
+/// Triangle-lower starting at off-diagonoal \p d (for compatibility with matlib)
+template <typename Matrix>
+typename traits::strict_lower<Matrix>::type
+inline tril(const Matrix& A, long d= 0)
+{
+    return bands(A, std::numeric_limits<long>::min(), d+1);
+}
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_STRICT_LOWER_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/strict_upper.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_STRICT_UPPER_INCLUDE
+#define MTL_MATRIX_STRICT_UPPER_INCLUDE
+
+#include <boost/numeric/mtl/matrix/bands.hpp>
+
+namespace mtl { namespace matrix {
+
+namespace traits {
+
+    template <typename Matrix>
+    struct strict_upper
+    {
+	typedef typename bands<Matrix>::type type;
+    };
+}
+
+///  Strict upper triangle matrix
+template <typename Matrix>
+typename traits::strict_upper<Matrix>::type
+inline strict_upper(const Matrix& A)
+{
+    return bands(A, 1, std::numeric_limits<long>::max());
+}
+
+/// Triangle-upper starting at off-diagonoal \p d (for compatibility with matlib)
+template <typename Matrix>
+typename traits::strict_upper<Matrix>::type
+inline triu(const Matrix& A, long d= 0)
+{
+    return bands(A, d, std::numeric_limits<long>::max());
+}
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_STRICT_UPPER_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/transposed_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_TRANSPOSED_VIEW_INCLUDE
+#define MTL_TRANSPOSED_VIEW_INCLUDE
+
+#include <utility>
+#include <boost/shared_ptr.hpp>
+#include <boost/mpl/if.hpp>
+#include <boost/type_traits/is_const.hpp>
+#include <boost/type_traits/is_same.hpp>
+#include <boost/type_traits/remove_const.hpp>
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/utility/transposed_orientation.hpp>
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/mtl/utility/property_map_impl.hpp>
+#include <boost/numeric/mtl/matrix/crtp_base_matrix.hpp>
+#include <boost/numeric/mtl/operation/sub_matrix.hpp>
+#include <boost/numeric/mtl/matrix/mat_expr.hpp>
+
+namespace mtl { namespace matrix {
+
+
+template <class Matrix>
+struct transposed_view
+  : public boost::mpl::if_<
+          boost::is_const<Matrix>
+        , const_crtp_base_matrix< const transposed_view<const Matrix>,
+				  typename Matrix::value_type, typename Matrix::size_type >
+        , crtp_base_matrix< transposed_view<Matrix>,
+			    typename Matrix::value_type, typename Matrix::size_type >
+      >::type,
+    public matrix::mat_expr< transposed_view<Matrix> >
+{
+    typedef transposed_view               self;
+    typedef mat_expr< self >              expr_base;
+    typedef Matrix                        other;
+
+    typedef typename mtl::traits::transposed_orientation<typename Matrix::orientation>::type orientation;
+    typedef typename Matrix::index_type                index_type;
+    typedef typename Matrix::value_type                value_type;
+    typedef typename Matrix::const_reference           const_reference;
+    typedef typename Matrix::key_type                  key_type;
+    typedef typename Matrix::size_type                 size_type;
+    typedef typename Matrix::dim_type::transposed_type dim_type;
+    typedef matrix::parameters<orientation, index_type, dim_type> parameters;
+
+    typedef typename boost::mpl::if_<boost::is_const<Matrix>,
+				     const_reference,
+				     value_type&
+				    >::type                  access_type;
+
+    typedef typename boost::mpl::if_<boost::is_const<Matrix>,
+				     const Matrix&,
+				     Matrix&
+				    >::type                  ref_type;
+    typedef const Matrix&                                    const_ref_type;
+
+
+    transposed_view (ref_type ref) : ref(ref) {}
+
+    transposed_view (const boost::shared_ptr<Matrix>& p) : my_copy(p), ref(*p) {}
+
+#ifdef MTL_WITH_CPP11_MOVE
+    transposed_view (self&& that) : my_copy(std::move(that.my_copy)), ref(that.ref) {}
+    transposed_view (const self& that) : ref(that.ref) { assert(that.my_copy.use_count() == 0); }
+#endif
+
+    const_reference operator() (size_type r, size_type c) const
+    {         return ref(c, r);     }
+
+    access_type operator() (size_type r, size_type c)
+    {         return ref(c, r);     }
+
+    size_type dim1() const  {  return ref.dim2(); }
+    size_type dim2() const  {  return ref.dim1(); }
+
+    dim_type dimensions() const
+    {        return ref.dimensions().transpose();  }
+
+    size_type begin_row() const { return ref.begin_col(); }
+    size_type end_row() const   { return ref.end_col();   }
+    size_type begin_col() const { return ref.begin_row(); }
+    size_type end_col() const   { return ref.end_row();   }
+    size_type nnz() const       { return ref.nnz();       }
+
+    friend size_type inline num_rows(const self& A)
+    { 	using mtl::matrix::num_cols; return num_cols(A.ref);     }
+    friend size_type inline num_cols(const self& A)
+    { 	using mtl::matrix::num_rows; return num_rows(A.ref);     }
+  protected:
+    boost::shared_ptr<Matrix>           my_copy;
+  public:
+    ref_type                            ref;
+};
+
+template <class Matrix>
+inline std::size_t size(const transposed_view<Matrix>& A)
+{
+    return num_rows(A.ref) * num_rows(A.ref);
+}
+
+// ==========
+// Sub matrix
+// ==========
+
+template <typename Matrix>
+struct sub_matrix_t< transposed_view<Matrix> >
+{
+    typedef transposed_view<Matrix>                                               matrix_type;
+    typedef typename boost::remove_const<Matrix>::type                            tmp_type;
+
+    // Transposed of submatrix type
+    typedef transposed_view<typename sub_matrix_t<tmp_type>::sub_matrix_type>       sub_matrix_type;
+    typedef transposed_view<typename sub_matrix_t<tmp_type>::const_sub_matrix_type> const_sub_matrix_type;
+    typedef typename matrix_type::size_type                                       size_type;
+
+    sub_matrix_type operator()(matrix_type& A, size_type begin_r, size_type end_r, size_type begin_c, size_type end_c)
+    {
+	typedef typename sub_matrix_t<Matrix>::sub_matrix_type   ref_sub_type;
+	typedef boost::shared_ptr<ref_sub_type>                  pointer_type;
+
+	// Submatrix of referred matrix, colums and rows interchanged
+	// Create a submatrix, whos address will be kept by transposed_view
+	pointer_type p(new ref_sub_type(sub_matrix(A.ref, begin_c, end_c, begin_r, end_r)));
+	return sub_matrix_type(p);
+    }
+
+    const_sub_matrix_type operator()(matrix_type const& A, size_type begin_r, size_type end_r,
+				     size_type begin_c, size_type end_c)
+    {
+	typedef typename sub_matrix_t<Matrix>::const_sub_matrix_type   ref_sub_type;
+	typedef boost::shared_ptr<ref_sub_type>                        pointer_type;
+
+	// Submatrix of referred matrix, colums and rows interchanged
+	// Create a submatrix, whos address will be kept by transposed_view
+	pointer_type p(new ref_sub_type(sub_matrix(A.ref, begin_c, end_c, begin_r, end_r)));
+	return const_sub_matrix_type(p);
+    }
+
+};
+
+}} // mtl::matrix
+
+namespace mtl { namespace traits {
+
+
+    namespace detail {
+
+	template <class Matrix, class Ref>
+	struct transposed_row
+	{
+	    typedef typename Matrix::key_type   key_type;
+	    typedef typename Matrix::size_type  size_type;
+
+	    transposed_row(mtl::matrix::transposed_view<Ref> const& transposed_matrix)
+	      : its_col(transposed_matrix.ref) {}
+
+	    size_type operator() (key_type const& key) const
+	    {
+		return its_col(key);
+	    }
+
+	  protected:
+	    typename col<typename boost::remove_const<Matrix>::type>::type  its_col;
+        };
+
+
+        template <class Matrix, class Ref>
+        struct transposed_col
+        {
+	    typedef typename Matrix::key_type   key_type;
+	    typedef typename Matrix::size_type  size_type;
+
+		transposed_col(mtl::matrix::transposed_view<Ref> const& transposed_matrix)
+	      : its_row(transposed_matrix.ref) {}
+
+	    size_type operator() (key_type const& key) const
+	    {
+		return its_row(key);
+	    }
+
+          protected:
+	    typename row<typename boost::remove_const<Matrix>::type>::type  its_row;
+        };
+
+    } // namespace detail
+
+    template <class Matrix>
+    struct row<matrix::transposed_view<Matrix> >
+    {
+	typedef detail::transposed_row<Matrix, Matrix>  type;
+    };
+
+    template <class Matrix>
+    struct col<matrix::transposed_view<Matrix> >
+    {
+	typedef detail::transposed_col<Matrix, Matrix>  type;
+    };
+
+    template <class Matrix>
+    struct const_value<mtl::matrix::transposed_view<Matrix> >
+    {
+	typedef mtl::detail::const_value_from_other<mtl::matrix::transposed_view<Matrix> > type;
+    };
+
+
+    template <class Matrix>
+    struct value<mtl::matrix::transposed_view<Matrix> >
+    {
+	typedef mtl::detail::value_from_other<mtl::matrix::transposed_view<Matrix> > type;
+    };
+
+
+// ================
+// Range generators
+// ================
+
+    namespace detail
+    {
+	template <class UseTag, class Matrix>
+	struct range_transposer_impl
+	{
+	    typedef range_generator<UseTag, typename boost::remove_const<Matrix>::type>  generator;
+	    // typedef range_generator<UseTag, Matrix>  generator;
+	    typedef typename generator::complexity   complexity;
+	    typedef typename generator::type         type;
+	    static int const                         level = generator::level;
+		type begin(mtl::matrix::transposed_view<Matrix> const& m)
+	    {
+		return generator().begin(m.ref);
+	    }
+		type end(mtl::matrix::transposed_view<Matrix> const& m)
+	    {
+		return generator().end(m.ref);
+	    }
+	};
+
+	// If considered range_generator for Matrix isn't supported, i.e. has infinite complexity
+	// then define as unsupported for transposed view
+	// (range_transposer_impl wouldn't compile in this case)
+	template <class UseTag, class Matrix>
+	struct range_transposer
+	    : boost::mpl::if_<
+	          boost::is_same<typename range_generator<UseTag, typename boost::remove_const<Matrix>::type>::complexity,
+				 complexity_classes::infinite>
+	        , range_generator<tag::unsupported, Matrix>
+	        , range_transposer_impl<UseTag, Matrix>
+	      >::type {};
+    }
+
+    // Row and column cursors are interchanged
+    template <class Matrix>
+    struct range_generator<glas::tag::col, matrix::transposed_view<Matrix> >
+      : detail::range_transposer<glas::tag::row, Matrix>
+    {};
+
+    template <class Matrix>
+    struct range_generator<glas::tag::row, matrix::transposed_view<Matrix> >
+      : detail::range_transposer<glas::tag::col, Matrix>
+    {};
+
+
+    // To traverse the major dimension refer to the Matrix
+    template <class Matrix>
+    struct range_generator<tag::major, matrix::transposed_view<Matrix> >
+      : detail::range_transposer<tag::major, Matrix>
+    {};
+
+
+    // Other cursors still use the same tag, e.g. elements
+    template <class Tag, class Matrix>
+    struct range_generator<Tag, matrix::transposed_view<Matrix> >
+      : detail::range_transposer<Tag, Matrix>
+    {};
+
+
+
+}} // namespace mtl::traits
+
+namespace mtl {
+	using matrix::transposed_view;
+}
+
+#endif // MTL_TRANSPOSED_VIEW_INCLUDE
+
+
+

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/upper.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_UPPER_INCLUDE
+#define MTL_MATRIX_UPPER_INCLUDE
+
+#include <boost/numeric/mtl/matrix/bands.hpp>
+#include <limits>
+
+namespace mtl { namespace matrix {
+
+namespace traits {
+
+    template <typename Matrix>
+    struct upper
+    {
+	typedef typename traits::bands<Matrix>::type type;
+    };
+}
+
+/// Upper triangular matrix
+template <typename Matrix>
+typename traits::upper<Matrix>::type
+inline upper(const Matrix& A)
+{
+    return bands(A, 0, std::numeric_limits<long>::max());
+}
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_UPPER_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/matrix/view_ref.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_VIEW_REF_INCLUDE
+#define MTL_MATRIX_VIEW_REF_INCLUDE
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/matrix/map_view.hpp>
+#include <boost/numeric/mtl/matrix/transposed_view.hpp>
+#include <boost/numeric/mtl/matrix/hermitian_view.hpp>
+
+
+namespace mtl { namespace matrix {
+
+template <typename Matrix>
+inline Matrix& view_ref(Matrix& A)
+{    return A; }
+
+template <typename Matrix>
+inline const Matrix& view_ref(const Matrix& A)
+{    return A; }
+
+template <typename Matrix>
+inline Matrix& view_ref(transposed_view<Matrix>& A)
+{    return A.ref; }
+
+template <typename Matrix>
+inline const Matrix& view_ref(transposed_view<const Matrix>& A)
+{    return A.ref; }
+
+template <typename Matrix>
+inline const Matrix& view_ref(const transposed_view<Matrix>& A)
+{    return A.ref; }
+
+template <typename Matrix>
+inline const Matrix& view_ref(const conj_view<Matrix>& A)
+{    return A.ref; }
+
+template <typename Matrix>
+inline const Matrix& view_ref(const hermitian_view<Matrix>& A)
+{    return A.const_ref(); }
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_VIEW_REF_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/mtl.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_MTL_INCLUDE
+#define MTL_MTL_INCLUDE
+
+#include <boost/numeric/mtl/types.hpp>
+#include <boost/numeric/mtl/operations.hpp>
+#include <boost/numeric/mtl/interfaces.hpp>
+
+#endif // MTL_MTL_INCLUDE