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

AMDiS/lib/mtl4/boost/numeric/mtl/mtl_fwd.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_FWD_INCLUDE
+#define MTL_MTL_FWD_INCLUDE
+
+/// Main name space for %Matrix Template Library
+namespace mtl {
+
+    template <typename T, typename U, typename Assign> struct lazy_assign;
+    template <typename VectorOut, typename Matrix, typename VectorIn, typename Assign> struct row_mat_cvec_index_evaluator;
+
+    /// Namespace for tags used for concept-free dispatching
+    namespace tag {
+	struct row_major;
+	struct col_major;
+
+	struct scalar;
+	struct vector;
+	struct matrix;
+
+	/// Namespace for constant iterator tags
+	namespace const_iter {}
+
+	/// Namespace for iterator tags
+	namespace iter {}
+    }
+    using tag::row_major;
+    using tag::col_major;
+
+    namespace index {
+	struct c_index;
+	struct f_index;
+    }
+
+    /// Namespace for compile-time parameters, e.g. %matrix dimensions
+    namespace fixed {
+	template <std::size_t Rows, std::size_t Cols> struct dimensions;
+    }
+
+    /// Namespace for run-time parameters, e.g. %matrix dimensions
+    namespace non_fixed {
+	struct dimensions;
+    }
+
+    /// Namespace for matrices and views and operations exclusively on matrices
+    namespace matrix {
+
+	template <typename Orientation, typename Index, typename Dimensions, bool OnStack, typename SizeType> struct parameters;
+
+        template <typename Value, typename Parameters> class dense2D;
+
+        template <typename Value, typename Parameters>
+        typename dense2D<Value, Parameters>::size_type num_cols(const dense2D<Value, Parameters>& matrix);
+        template <typename Value, typename Parameters>
+        typename dense2D<Value, Parameters>::size_type num_rows(const dense2D<Value, Parameters>& matrix);
+        template <typename Value, typename Parameters>
+        typename dense2D<Value, Parameters>::size_type size(const dense2D<Value, Parameters>& matrix);
+
+        template <typename Value, std::size_t Mask, typename Parameters> class morton_dense;
+
+#if !defined(_MSC_VER) || _MSC_VER != 1400 // Trouble in MSVC 2005
+        template <typename Value, std::size_t Mask, typename Parameters>
+        typename morton_dense<Value, Mask, Parameters>::size_type num_cols(const morton_dense<Value, Mask, Parameters>& matrix);
+        template <typename Value, std::size_t Mask, typename Parameters>
+        typename morton_dense<Value, Mask, Parameters>::size_type num_rows(const morton_dense<Value, Mask, Parameters>& matrix);
+        template <typename Value, std::size_t Mask, typename Parameters>
+        typename morton_dense<Value, Mask, Parameters>::size_type size(const morton_dense<Value, Mask, Parameters>& matrix);
+#endif
+
+        template <typename Value, typename Parameters> class compressed2D;
+
+        template <typename Value, typename Parameters>
+        typename compressed2D<Value, Parameters>::size_type num_cols(const compressed2D<Value, Parameters>& matrix);
+        template <typename Value, typename Parameters>
+        typename compressed2D<Value, Parameters>::size_type num_rows(const compressed2D<Value, Parameters>& matrix);
+        template <typename Value, typename Parameters>
+        // typename compressed2D<Value, Parameters>::size_type
+	std::size_t
+	size(const compressed2D<Value, Parameters>& matrix);
+
+        template <typename Value, typename Parameters, typename Updater> struct compressed2D_inserter;
+
+	template <typename T, typename Parameters> class coordinate2D;
+	template <typename Matrix, typename Updater> struct coordinate2D_inserter;
+
+	template <typename T, typename Parameters> class sparse_banded;
+	template <typename T, typename Parameters, typename Updater> struct sparse_banded_inserter;
+
+	template <typename Value, typename Parameters> class ell_matrix;
+	template <typename Value, typename Parameters, typename Updater> struct ell_matrix_inserter;
+
+	template <typename Matrix, typename Updater> struct inserter;
+	template <typename BaseInserter> class shifted_inserter;
+
+	template <typename Vector> class multi_vector;
+	template <typename Vector> class multi_vector_range;
+	template <typename Value> class element;
+	template <typename Value> class element_structure;
+	template <typename Functor> class implicit_dense;
+	template <typename Value> class ones_functor;
+	template <typename Value> class ones_matrix;
+	template <typename Value> class hilbert_functor;
+	template <typename Value> class hilbert_matrix;
+	template <typename Vector1, typename Vector2> class outer_product_functor;
+	template <typename Vector1, typename Vector2> class outer_product_matrix;
+
+        template <typename Matrix> struct transposed_view;
+
+	template <typename Matrix> struct mat_expr;
+	template <typename Matrix> struct dmat_expr;
+	template <typename Matrix> struct smat_expr;
+	template <typename M1, typename M2, typename SFunctor> struct mat_mat_op_expr;
+	template <typename M1, typename M2> struct mat_mat_plus_expr;
+	template <typename M1, typename M2> struct mv_mv_plus_expr;
+	template <typename M1, typename M2> struct mat_mat_minus_expr;
+	template <typename M1, typename M2> struct mv_mv_minus_expr;
+	template <typename M1, typename M2> struct mat_mat_ele_times_expr;
+	template <typename M1, typename M2> struct mat_mat_times_expr;
+	template <typename M1, typename M2> struct mat_mat_asgn_expr;
+
+	template <typename Matrix> struct mat_expr;
+	template <typename Functor, typename Matrix> struct map_view;
+	template <typename Scaling, typename Matrix> struct scaled_view;
+	template <typename Matrix, typename RScaling> struct rscaled_view; // added by Hui Li
+	template <typename Matrix, typename Divisor> struct divide_by_view; // added by Hui Li
+	template <typename Matrix>  struct conj_view;
+	template <typename Matrix>  struct negate_view;
+	template <typename Matrix>  struct imag_view;
+	template <typename Matrix>  struct real_view;
+	template <typename Matrix>  struct hermitian_view;
+	template <typename Matrix>  struct banded_view;
+	template <typename Matrix> struct indirect;
+
+	template <typename Matrix> class lu_solver;
+
+	template <typename Matrix> std::size_t size(const banded_view<Matrix>&);
+	template <class Matrix> std::size_t size(const transposed_view<Matrix>&);
+	template <typename Functor, typename Matrix> std::size_t size(const map_view<Functor, Matrix>&);
+	template <class Matrix> std::size_t size(const hermitian_view<Matrix>& );
+   }
+
+    //using matrix::dense2D;
+    //using matrix::morton_dense;
+    //using matrix::compressed2D;
+    //using matrix::coordinate2D;
+    //using matrix::multi_vector;
+    //using matrix::transposed_view;
+
+    template <typename E1, typename E2> struct mat_cvec_times_expr;
+
+    /// Namespace for vectors and views and %operations exclusively on vectors
+    namespace vector {
+	template <typename Vector> struct vec_expr;
+	template <typename Value, typename Parameters> class dense_vector;
+	template <typename Value, typename Parameters> class strided_vector_ref;
+	template <typename Value, typename Parameters> class sparse_vector;
+	template <typename Functor, typename Vector> struct map_view;
+	template <typename Vector>  struct conj_view;
+	template <typename Vector>  struct real_view;
+	template <typename Vector>  struct imag_view;
+	template <typename Vector>  struct negate_view;
+	template <typename Scaling, typename Vector> struct scaled_view;
+	template <typename Vector, typename RScaling> struct rscaled_view; // added by Hui Li
+	template <typename Vector, typename Divisor> struct divide_by_view; // added by Hui Li
+	template <class E1, class E2, typename SFunctor> struct vec_vec_op_expr;
+	template <class E1, class E2, typename SFunctor> struct vec_vec_pmop_expr;
+	template <class E1, class E2, typename SFunctor> struct vec_vec_aop_expr;
+	template <class E1, class E2, typename SFunctor> struct vec_vec_ele_prod_expr;
+	template <class E1, class E2, typename SFunctor> struct vec_scal_aop_expr;
+	template <class E1, class E2> struct vec_vec_asgn_expr;
+	template <class E1, class E2> struct vec_vec_plus_asgn_expr;
+	template <class E1, class E2> struct vec_vec_minus_asgn_expr;
+	template <class E1, class E2> struct vec_vec_times_asgn_expr; // is this really used???
+//	template <class E1, class E2> struct vec_vec_div_asgn_expr; // is this really used???
+	template <class E1, class E2> struct vec_scal_times_asgn_expr;
+	template <class E1, class E2> struct vec_scal_div_asgn_expr; // added by Hui Li
+	template <class E1, class E2> struct vec_scal_asgn_expr;
+	template <typename E1, typename E2> struct rvec_mat_times_expr;
+
+	template <typename Vector> struct vec_const_ref_expr;
+	template <unsigned BSize, typename Vector> class unrolled1;
+
+	template <typename Scalar, typename Vector, typename Functor, typename Assign> struct reduction_index_evaluator;
+	template <typename Scalar, typename Vector1, typename Vector2, typename ConjOpt, typename Assign> struct dot_index_evaluator;
+
+	template <typename Vector, typename Functor> struct lazy_reduction;
+
+	template <typename Matrix, typename VectorIn> struct mat_cvec_multiplier;
+
+	template <typename Value, typename Parameters, typename Value2>
+	inline void fill(dense_vector<Value, Parameters>& vector, const Value2& value);
+
+	template <typename Value, typename Parameters>
+	typename dense_vector<Value, Parameters>::size_type
+	inline size(const dense_vector<Value, Parameters>& vector);
+
+	template <typename Value, typename Parameters>
+	typename dense_vector<Value, Parameters>::size_type
+	inline num_rows(const dense_vector<Value, Parameters>& vector);
+
+	template <typename Value, typename Parameters>
+	typename dense_vector<Value, Parameters>::size_type
+	inline num_cols(const dense_vector<Value, Parameters>& vector);
+
+	template <typename Value, typename Parameters>
+	std::size_t size(const strided_vector_ref<Value, Parameters>& v);
+
+	template <typename Functor, typename Vector>
+	std::size_t size(const map_view<Functor, Vector>& v);
+
+	template <typename E1, typename E2, typename SFunctor>
+	std::size_t size(const vec_vec_op_expr<E1, E2, SFunctor>& v);
+
+	template <typename E1, typename E2, typename SFunctor>
+	std::size_t size(const vec_vec_aop_expr<E1, E2, SFunctor>& v);
+
+	template <typename E1, typename E2, typename SFunctor>
+	std::size_t size(const vec_vec_pmop_expr<E1, E2, SFunctor>& v);
+
+	template <typename E1, typename E2, typename SFunctor>
+	inline std::size_t size(const vec_scal_aop_expr<E1, E2, SFunctor>& v);
+
+	template <unsigned BSize, typename Vector>
+	inline std::size_t size(const unrolled1<BSize, Vector>& v);
+
+	template <typename E1, typename E2>
+	std::size_t inline size(const rvec_mat_times_expr<E1, E2>& x);
+
+	template <typename E1, typename E2>
+	std::size_t inline size(const mat_cvec_times_expr<E1, E2>& x);
+
+	template <typename Matrix, typename VectorIn>
+	std::size_t size(const mat_cvec_multiplier<Matrix, VectorIn>& m);
+
+	/// Namespace for fixed vector dimension types
+	namespace fixed {
+	    template <std::size_t Size> struct dimension;
+	}
+	/// Namespace for non-fixed vector dimension types, i.e. size dynamically determined at run-time
+	namespace non_fixed {
+	    struct dimension;
+	}
+
+    }
+
+    //using vector::dense_vector;
+
+    // Export free vector functions into mtl namespace
+    // It is also needed to handle STL vectors in MTL
+    using vector::fill;
+    using vector::size;
+    using vector::num_rows;
+    using vector::num_cols;
+
+    /// Namespace for %operations (if not defined in mtl)
+    namespace operations {
+	template <typename T> struct update_store;
+    }
+
+
+    namespace vector {
+
+	template <typename Vector, typename Updater = mtl::operations::update_store<typename Vector::value_type> > struct inserter;
+	template <typename Vector, typename Size> struct update_proxy;
+    }
+
+
+
+
+
+    /// Namespace for type %traits
+    namespace traits {
+	template <typename Value> struct category;
+	template <typename Value> struct algebraic_category;
+
+	template <typename Collection> struct value;
+	template <typename Collection> struct const_value;
+	template <typename Collection> struct row;
+	template <typename Collection> struct col;
+	template <class Matrix> struct offset;
+
+	template <class Vector> struct index;
+	template <typename Tag, typename Collection>  struct range_generator;
+
+	template <typename T> struct eval_dense;
+
+        template <typename Matrix> struct transposed_orientation;
+
+	template <typename Collection> struct root;
+
+	// for internal implementations
+	namespace detail {
+	    // needed collection.hpp (at least)
+	    template <typename Collection, typename Cursor, typename Complexity> struct dense_element_range_generator;
+	    template <typename Matrix, typename Cursor, typename Complexity>     struct all_offsets_range_generator;
+	    template <typename Matrix, typename Tag, int Level = 2>              struct sub_matrix_cursor;
+	    template <typename Matrix>                                           struct matrix_element_key;
+	    template <typename Matrix, int pos>                                  struct matrix_element_cursor;
+	    template <typename Matrix, typename Complexity, int Level = 2>       struct all_rows_range_generator;
+	    template <typename Cursor>                                           struct all_cols_in_row_range_generator;
+	    template <typename Matrix, typename Complexity, int Level = 2>       struct all_cols_range_generator;
+	    template <typename Cursor>                                           struct all_rows_in_col_range_generator;
+	    template <typename Collection, typename RangeGenerator>              struct referred_range_generator;
+	}
+    }
+
+    template <typename Matrix> struct ColumnInMatrix;
+    template <typename Matrix> struct RowInMatrix;
+
+    template <class Tag, class Collection> typename traits::range_generator<Tag, Collection>::type
+    begin(Collection const& c);
+
+    template <class Tag, class Collection> typename traits::range_generator<Tag, Collection>::type
+    end(Collection const& c);
+
+
+    /// Namespace for functors with application operator and fully typed parameters
+    namespace tfunctor {
+	/// Functor for scaling matrices, vectors and ordinary scalars
+	template <typename V1, typename V2, typename AlgebraicCategory = tag::scalar> struct scale;
+
+	/// Functor for scaling matrices, vectors and ordinary scalars
+	template <typename V1, typename V2, typename AlgebraicCategory = tag::scalar> struct rscale;
+
+	/// Functor for scaling matrices, vectors and ordinary scalars
+	template <typename V1, typename V2, typename AlgebraicCategory = tag::scalar> struct divide_by;
+    }
+
+    /// Namespace for functors with static function apply and fully typed parameters
+    namespace sfunctor {
+	template <typename Value, typename AlgebraicCategory = tag::scalar> struct conj_aux;
+	template <typename Value>                   struct conj;
+	template <typename Value>                   struct imag;
+	template <typename Value>                   struct real;
+	template <typename Value>                   struct negate;
+	template <typename Value1, typename Value2> struct plus;
+	template <typename Value1, typename Value2> struct minus;
+	template <typename Value1, typename Value2> struct times;
+	template <typename Value1, typename Value2> struct divide;
+	template <typename Value1, typename Value2> struct assign;
+	template <typename Value1, typename Value2> struct plus_assign;
+	template <typename Value1, typename Value2> struct minus_assign;
+	template <typename Value1, typename Value2> struct times_assign;
+	template <typename Value1, typename Value2> struct divide_assign;
+	template <typename Value>                   struct identity;
+	template <typename Value>                   struct abs;
+	template <typename Value>                   struct sqrt;
+	template <typename Value>                   struct square;
+	template <typename F, typename G>           struct compose;
+	template <typename F, typename G>           struct compose_first;
+	template <typename F, typename G>           struct compose_second;
+	template <typename F, typename G, typename H> struct compose_both;
+	template <typename F, typename G>           struct compose_binary;
+    }
+
+    // Namespace documentations
+
+    /// Namespace for static assignment functors
+    namespace assign {}
+
+    /// Namespace for complexity classes
+    namespace complexity_classes {}
+
+    /// Namespace for %operations (if not defined in mtl)
+    namespace operations {}
+
+    /// Namespace for recursive operations and types with recursive memory layout
+    namespace recursion {}
+
+    /// Namespace for %utilities
+    namespace utility {}
+
+    /// Namespace for implementations using recursators
+    namespace wrec {}
+
+    namespace matrix {
+	template <typename Matrix, typename ValueType, typename SizeType> struct crtp_matrix_assign;
+	template <typename Matrix, typename ValueType, typename SizeType> struct const_crtp_matrix_bracket;
+	template <typename Matrix, typename ValueType, typename SizeType> struct crtp_matrix_bracket;
+	template <typename Matrix, typename ValueType, typename SizeType> struct const_crtp_base_matrix;
+	template <typename Matrix, typename ValueType, typename SizeType> struct crtp_base_matrix;
+	template <typename Matrix> struct const_crtp_matrix_range_bracket;
+    }
+
+    namespace detail {
+	template <typename Value, bool OnStack, unsigned Size= 0> struct contiguous_memory_block;
+	template <typename Matrix, typename Updater> struct trivial_inserter;
+	template <typename Collection> struct with_format_t;
+    }
+
+    /// Free function defined for all matrix and vector types
+    template <typename Collection> void swap(Collection& c1, Collection& c2);
+
+    /// User registration that class has a clone constructor, otherwise use regular copy constructor.
+    template<typename T> struct is_clonable;
+
+    /// Helper type to define constructors that always copy
+    struct clone_ctor {};
+
+    /// Helper type to define constructors that refer to another object instead of copying it (e.g. transposed vectors)
+    struct refer_ctor {};
+
+    class irange;
+    class iset;
+    class srange;
+
+    template <typename T, typename U> struct fused_expr;
+    namespace vector {
+	template <typename T, typename U> struct fused_index_evaluator;
+	// template <typename T, typename U> size_t size(const fused_index_evaluator<T, U>&); // not needed currently
+    }
+
+    /// Namespace for I/O operations
+    namespace io {
+	class matrix_market_istream;
+	class matrix_market_ostream;
+
+	template <typename MatrixIStream, typename MatrixOStream> class matrix_file;
+	typedef matrix_file<matrix_market_istream, matrix_market_ostream> matrix_market;
+    }
+
+    // Multiplication functors
+    template <typename Assign, typename Backup> struct gen_cursor_dmat_dmat_mult_t;
+    template <typename Assign, typename Backup> struct gen_dmat_dmat_mult_t;
+    template <unsigned long Tiling1, unsigned long Tiling2, typename Assign, typename Backup> struct gen_tiling_dmat_dmat_mult_t;
+    template <typename Assign, typename Backup> struct gen_tiling_44_dmat_dmat_mult_t;
+    template <typename Assign, typename Backup> struct gen_tiling_22_dmat_dmat_mult_t;
+    template <typename BaseMult, typename BaseTest, typename Assign, typename Backup> struct gen_recursive_dmat_dmat_mult_t;
+    template <typename Assign, typename Backup> struct gen_platform_dmat_dmat_mult_t;
+    template <typename Assign, typename Backup> struct gen_blas_dmat_dmat_mult_t;
+    template <std::size_t SizeLimit, typename FunctorSmall, typename FunctorLarge> struct size_switch_dmat_dmat_mult_t;
+    template <bool IsStatic, typename FunctorStatic, typename FunctorDynamic> struct static_switch_dmat_dmat_mult_t;
+    template <typename Assign, typename Backup> struct fully_unroll_fixed_size_dmat_dmat_mult_t;
+
+} // namespace mtl
+
+#endif // MTL_MTL_FWD_INCLUDE
+
+

AMDiS/lib/mtl4/boost/numeric/mtl/operation/adjoint.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_ADJOINT_INCLUDE
+#define MTL_ADJOINT_INCLUDE
+
+#include <boost/numeric/mtl/operation/hermitian.hpp>
+
+namespace mtl { namespace matrix {
+
+namespace traits {
+
+    template <typename LinOp>
+    struct adjoint
+    {
+	typedef typename mtl::matrix::detail::hermitian<LinOp>::result_type  type;
+	type operator()(const LinOp& A)
+	{
+	    return hermitian(A);
+	}
+    };
+}
+
+/// Adjoint linear operator, typically Hermitian transposed
+template <typename LinOp>
+typename traits::adjoint<LinOp>::type
+inline adjoint(const LinOp& A)
+{
+    return traits::adjoint<LinOp>()(A);
+}
+
+}} // namespace mtl::matrix
+
+#endif // MTL_ADJOINT_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/adjust_cursor.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_DETAIL_ADJUST_CURSOR_INCLUDE
+#define MTL_DETAIL_ADJUST_CURSOR_INCLUDE
+
+namespace mtl { namespace detail {
+
+    template <typename Size, typename Cursor>
+    void inline adjust_cursor(Size diff, Cursor& c, tag::dense) { c+= diff; }
+
+    template <typename Size, typename Cursor>
+    void inline adjust_cursor(Size diff, Cursor& c, tag::sparse) {}
+
+
+}} // namespace mtl::detail
+
+#endif // MTL_DETAIL_ADJUST_CURSOR_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/assign_each_nonzero.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_ASSIGN_EACH_NONZERO_INCLUDE
+#define MTL_ASSIGN_EACH_NONZERO_INCLUDE
+
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/range_generator.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+
+namespace mtl {
+
+    namespace vector {
+
+	/// Assign the result of \p f(i) to each non-zero i in non-constant vector \p v
+	template <typename Vector, typename Functor>
+	inline void assign_each_nonzero(Vector& v, const Functor& f)
+	{
+	    vampir_trace<3008> tracer;
+	    typedef typename traits::range_generator<tag::iter::nz, Vector>::type iterator;
+	    for (iterator i= begin<tag::iter::nz>(v), iend= end<tag::iter::nz>(v); i != iend; ++i)
+		*i= f(*i);
+	}
+
+    } // namespace vector
+
+    namespace matrix {
+
+	/// Assign the result of \p f(i) to each non-zero i in non-constant matrix \p A
+	template <typename Matrix, typename Functor>
+	inline void assign_each_nonzero(Matrix& m, const Functor& f)
+	{
+	    vampir_trace<3008> tracer;
+	    typename mtl::traits::value<Matrix>::type     value(m);
+
+	    typedef typename mtl::traits::range_generator<tag::major, Matrix>::type     cursor_type;
+	    typedef typename mtl::traits::range_generator<tag::nz, cursor_type>::type   icursor_type;
+
+	    for (cursor_type cursor = begin<tag::major>(m), cend = end<tag::major>(m); cursor != cend; ++cursor)
+		for (icursor_type icursor = begin<tag::nz>(cursor), icend = end<tag::nz>(cursor);
+		     icursor != icend; ++icursor)
+			value(*icursor, f(value(*icursor)));
+	}
+
+    } // namespace matrix
+
+} // namespace mtl
+
+#endif // MTL_ASSIGN_EACH_NONZERO_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/assign_mode.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_ASSIGN_MODE_INCLUDE
+#define MTL_ASSIGN_MODE_INCLUDE
+
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+
+namespace mtl { namespace assign {
+
+/// Functor for assigning the result of a calculation (e.g. a sum) to a typically scalar value
+struct assign_sum
+{
+    static const bool init_to_zero= true;
+
+    template <typename T>
+    static void init(T& v)
+    {
+	using math::zero;
+	v= zero(v);
+    }
+
+    // Sets x with y when empty or not initialized
+    template <typename T, typename U>
+    static void set_empty(T& x, const U& y)
+    {
+	x= y;
+    }
+
+    // The first update sets the value and avoids the zeroing
+    template <typename T, typename U>
+    static void first_update(T& x, const U& y)
+    {
+	x= y;
+    }
+
+    template <typename T, typename U>
+    static void update(T& x, const U& y)
+    {
+	x+= y;
+    }
+
+    // To be used like sfunctor
+    template <typename T, typename U>
+    static T& apply(T& x, const U& y)
+    {
+	return x= y;
+    }
+
+};
+
+/// Functor for incrementing a typically scalar value with the result of a calculation (e.g. a sum)
+struct plus_sum
+{
+    static const bool init_to_zero= false;
+
+    template <typename T>
+    static void init(T&) {}
+
+    // Sets x with y when empty or not initialized
+    template <typename T, typename U>
+    static void set_empty(T& x, const U& y)
+    {
+	x= y;
+    }
+
+    template <typename T, typename U>
+    static void first_update(T& x, const U& y)
+    {
+	x+= y;
+    }
+
+    template <typename T, typename U>
+    static void update(T& x, const U& y)
+    {
+	x+= y;
+    }
+
+    // To be used like sfunctor
+    template <typename T, typename U>
+    static T& apply(T& x, const U& y)
+    {
+	return x+= y;
+    }
+};
+
+
+/// Functor for incrementing a typically scalar value with the result of a calculation (e.g. a sum)
+struct minus_sum
+{
+    static const bool init_to_zero= false;
+
+    template <typename T>
+    static void init(T&) {}
+
+    // Sets x with y when empty or not initialized
+    template <typename T, typename U>
+    static void set_empty(T& x, const U& y)
+    {
+	x= -y;
+    }
+
+    template <typename T, typename U>
+    static void first_update(T& x, const U& y)
+    {
+	x-= y;
+    }
+
+    template <typename T, typename U>
+    static void update(T& x, const U& y)
+    {
+	x-= y;
+    }
+
+    // To be used like sfunctor
+    template <typename T, typename U>
+    static T& apply(T& x, const U& y)
+    {
+	return x-= y;
+    }
+};
+
+
+
+/// Functor for multiplying a typically scalar value with the result of a calculation (e.g. a sum)
+struct times_sum
+{
+    static const bool init_to_zero= false;
+
+    template <typename T>
+    static void init(T&) {}
+
+    // Sets x with y when empty or not initialized
+    template <typename T, typename U>
+    static void set_empty(T& x, const U&)
+    {
+	x= T(0);
+    }
+
+    template <typename T, typename U>
+    static void first_update(T& x, const U& y)
+    {
+	x*= y;
+    }
+
+    template <typename T, typename U>
+    static void update(T& x, const U& y)
+    {
+	x*= y;
+    }
+
+    // To be used like sfunctor
+    template <typename T, typename U>
+    static T& apply(T& x, const U& y)
+    {
+	return x*= y;
+    }
+};
+
+/// Functor for dividing a typically scalar value with the result of a calculation (e.g. a sum)
+struct divide_sum
+{
+    static const bool init_to_zero= false;
+
+    template <typename T>
+    static void init(T&) {}
+
+    // Sets x with y when empty or not initialized
+    template <typename T, typename U>
+    static void set_empty(T& x, const U& y)
+    {
+        MTL_DEBUG_THROW_IF(y == U(0), division_by_zero());
+	x= T(0);
+    }
+
+    template <typename T, typename U>
+    static void first_update(T& x, const U& y)
+    {
+        MTL_DEBUG_THROW_IF(y == U(0), division_by_zero());
+	x/= y;
+    }
+
+    template <typename T, typename U>
+    static void update(T& x, const U& y)
+    {
+        MTL_DEBUG_THROW_IF(y == U(0), division_by_zero());
+	x/= y;
+    }
+
+    // To be used like sfunctor
+    template <typename T, typename U>
+    static T& apply(T& x, const U& y)
+    {
+	return x/= y;
+    }
+};
+
+}} // namespace mtl::assign
+
+#endif // MTL_ASSIGN_MODE_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/bin_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_BIN_OP_EXPR_INCLUDE
+#define MTL_BIN_OP_EXPR_INCLUDE
+
+namespace mtl {
+
+/// Minimalistic expression template for binary operation: keeps only references.
+template <typename E1, typename E2>
+struct bin_op_expr
+{
+    typedef bin_op_expr                          self;
+
+    typedef E1                                   first_argument_type ;
+    typedef E2                                   second_argument_type ;
+
+    bin_op_expr( first_argument_type const& v1, second_argument_type const& v2 )
+	: first( v1 ), second( v2 )
+    {}
+
+    first_argument_type const&  first ;
+    second_argument_type const& second ;
+};
+
+} // namespace mtl
+
+#endif // MTL_BIN_OP_EXPR_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/cholesky.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_CHOLESKY_INCLUDE
+#define MTL_CHOLESKY_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/recursion/matrix_recursator.hpp>
+#include <boost/numeric/mtl/utility/glas_tag.hpp>
+#include <boost/numeric/mtl/utility/range_generator.hpp>
+#include <boost/numeric/mtl/operation/dmat_dmat_mult.hpp>
+#include <boost/numeric/mtl/operation/assign_mode.hpp>
+#include <boost/numeric/mtl/matrix/transposed_view.hpp>
+#include <boost/numeric/mtl/recursion/base_case_cast.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace mtl { namespace matrix {
+
+namespace with_bracket {
+
+    // ============================================================================
+    // Generic Cholesky factorization and operands for Cholesky on with submatrices
+    // ============================================================================
+
+    template < typename Matrix >
+    void cholesky_base (Matrix & matrix)
+    {
+	vampir_trace<5001> tracer;
+	typedef typename  Collection<Matrix>::size_type    size_type;
+	for (size_type k = 0; k < matrix.num_cols(); k++) {
+	    matrix[k][k] = sqrt (matrix[k][k]);
+
+		for (size_type i = k + 1; i < matrix.num_rows(); i++) {
+		    matrix[i][k] /= matrix[k][k];
+		    typename Collection<Matrix>::value_type d = matrix[i][k];
+
+		    for (size_type j = k + 1; j <= i; j++)
+			matrix[i][j] -= d * matrix[j][k];
+		}
+	}
+    }
+
+
+    template < typename MatrixSW, typename MatrixNW >
+    void tri_solve_base(MatrixSW & SW, const MatrixNW & NW)
+    {
+	vampir_trace<5002> tracer;
+	typedef typename  Collection<MatrixSW>::size_type    size_type;
+	for (size_type k = 0; k < NW.num_rows (); k++) {
+
+	    for (size_type i = 0; i < SW.num_rows (); i++) {
+		SW[i][k] /= NW[k][k];
+		typename MatrixSW::value_type d = SW[i][k];
+
+		for (size_type j = k + 1; j < SW.num_cols (); j++)
+		    SW[i][j] -= d * NW[j][k];
+	    }
+	}
+    }
+
+
+    // Lower(SE) -= SW * SW^T
+    template < typename MatrixSE, typename MatrixSW >
+    void tri_schur_base(MatrixSE & SE, const MatrixSW & SW)
+    {
+	vampir_trace<5003> tracer;
+	typedef typename  Collection<MatrixSE>::size_type    size_type;
+	for (size_type k = 0; k < SW.num_cols (); k++)
+
+	    for (size_type i = 0; i < SE.num_rows (); i++) {
+		    typename MatrixSW::value_type d = SW[i][k];
+		    for (size_type j = 0; j <= i; j++)
+			SE[i][j] -= d * SW[j][k];
+	    }
+    }
+
+
+    template < typename MatrixNE, typename MatrixNW, typename MatrixSW >
+    void schur_update_base(MatrixNE & NE, const MatrixNW & NW, const MatrixSW & SW)
+    {
+	vampir_trace<5004> tracer;
+	typedef typename  Collection<MatrixNE>::size_type    size_type;
+	for (size_type k = 0; k < NW.num_cols (); k++)
+	    for (size_type i = 0; i < NE.num_rows (); i++) {
+		typename MatrixNW::value_type d = NW[i][k];
+		for (size_type j = 0; j < NE.num_cols (); j++)
+		    NE[i][j] -= d * SW[j][k];
+	    }
+    }
+
+
+    // ======================
+    // Corresponding functors
+    // ======================
+
+    struct cholesky_base_t
+    {
+	template < typename Matrix >
+	void operator() (Matrix & matrix)
+	{
+	    cholesky_base(matrix);
+	}
+    };
+
+    struct tri_solve_base_t
+    {
+	template < typename MatrixSW, typename MatrixNW >
+	void operator() (MatrixSW & SW, const MatrixNW & NW)
+	{
+	    tri_solve_base(SW, NW);
+	}
+    };
+
+    struct tri_schur_base_t
+    {
+	template < typename MatrixSE, typename MatrixSW >
+	void operator() (MatrixSE & SE, const MatrixSW & SW)
+	{
+	    tri_schur_base(SE, SW);
+	}
+    };
+
+    struct schur_update_base_t
+    {
+	template < typename MatrixNE, typename MatrixNW, typename MatrixSW >
+	void operator() (MatrixNE & NE, const MatrixNW & NW, const MatrixSW & SW)
+	{
+	    schur_update_base(NE, NW, SW);
+	}
+    };
+
+} // namespace with_bracket
+
+#if 0
+namespace with_iterator {
+
+    // ============================================================================
+    // Generic Cholesky factorization and operands for Cholesky on with submatrices
+    // ============================================================================
+    // CURRENTLY NOT SUPPORTED -- CAUSES SEGFAULT, e.g. with icc 11.0 in r8536 !!!!
+    // ============================================================================
+
+    template < typename Matrix >
+    void cholesky_base (Matrix& matrix)
+    {
+	vampir_trace<5001> tracer;
+	typedef typename  Collection<Matrix>::size_type    size_type;
+
+	using namespace glas::tag; using mtl::traits::range_generator;
+	typedef tag::iter::all    all_it;
+
+        typedef typename Collection<Matrix>::value_type                         value_type;
+        typedef typename range_generator<col, Matrix>::type       cur_type;
+        typedef typename range_generator<all_it, cur_type>::type    iter_type;
+
+	typedef typename range_generator<row, Matrix>::type       rcur_type;
+	typedef typename range_generator<all_it, rcur_type>::type   riter_type;
+
+	size_type k= 0;
+	for (cur_type kb= begin<col>(matrix), kend= end<col>(matrix); kb != kend; ++kb, ++k) {
+
+	    iter_type ib= begin<all_it>(kb), iend= end<all_it>(kb);
+	    ib+= k; // points now to matrix[k][k]
+
+	    value_type root= sqrt (*ib);
+	    *ib= root;
+
+	    ++ib; // points now to matrix[k+1][k]
+	    rcur_type rb= begin<row>(matrix); rb+= k+1; // to row k+1
+	    for (size_type i= k + 1; ib != iend; ++ib, ++rb, ++i) {
+		*ib = *ib / root;
+		typename Collection<Matrix>::value_type d = *ib;
+		riter_type it1= begin<all_it>(rb);    it1+= k+1;      // matrix[i][k+1]
+		riter_type it1end= begin<all_it>(rb); it1end+= i+1;   // matrix[i][i+1]
+		iter_type it2= begin<all_it>(kb);     it2+= k+1;      // matrix[k+1][k]
+		for (; it1 != it1end; ++it1, ++it2)
+		    *it1 = *it1 - d * *it2;
+	    }
+	}
+    }
+
+
+    template < typename MatrixSW, typename MatrixNW >
+    void tri_solve_base(MatrixSW & SW, const MatrixNW & NW)
+    {
+	vampir_trace<5002> tracer;
+	typedef typename  Collection<MatrixSW>::size_type    size_type;
+
+	using namespace glas::tag; using mtl::traits::range_generator;
+	typedef tag::iter::all        all_it;
+	typedef tag::const_iter::all  all_cit;
+
+        typedef typename range_generator<col, MatrixNW>::type       ccur_type;
+        typedef typename range_generator<all_cit, ccur_type>::type    citer_type;
+
+	typedef typename range_generator<row, MatrixSW>::type       rcur_type;
+	typedef typename range_generator<all_it, rcur_type>::type     riter_type;
+
+	for (size_type k = 0; k < NW.num_rows (); k++)
+	    for (size_type i = 0; i < SW.num_rows (); i++) {
+
+		typename MatrixSW::value_type d = SW[i][k] /= NW[k][k];
+
+		rcur_type sw_i= begin<row>(SW);     sw_i+= i;  // row i
+		riter_type it1= begin<all_it>(sw_i);  it1+= k+1; // SW[i][k+1]
+		riter_type it1end= end<all_it>(sw_i);
+
+		ccur_type nw_k= begin<col>(NW);     nw_k+= k;  // column k
+		citer_type it2= begin<all_cit>(nw_k); it2+= k+1; // NW[k+1][k]
+
+		for(; it1 != it1end; ++it1, ++it2)
+		    *it1 = *it1 - d * *it2;
+	    }
+    }
+
+
+    // Lower(SE) -= SW * SW^T
+    template < typename MatrixSE, typename MatrixSW >
+    void tri_schur_base(MatrixSE & SE, const MatrixSW & SW)
+    {
+	vampir_trace<5003> tracer;
+	typedef typename  Collection<MatrixSE>::size_type    size_type;
+
+	using namespace glas::tag; using mtl::traits::range_generator;
+	typedef tag::iter::all        all_it;
+	typedef tag::const_iter::all  all_cit;
+
+        typedef typename range_generator<col, MatrixSW>::type       ccur_type;
+        typedef typename range_generator<all_cit, ccur_type>::type    citer_type;
+
+	typedef typename range_generator<row, MatrixSE>::type       rcur_type;
+	typedef typename range_generator<all_it, rcur_type>::type     riter_type;
+
+	for (size_type k = 0; k < SW.num_cols (); k++)
+	    for (size_type i = 0; i < SE.num_rows (); i++) {
+		typename MatrixSW::value_type d = SW[i][k];
+
+		rcur_type se_i= begin<row>(SE);       se_i+= i;      // row i
+		riter_type it1= begin<all_it>(se_i);                   // SE[i][0]
+		riter_type it1end= begin<all_it>(se_i); it1end+= i+1;  // SE[i][i+i]
+
+		ccur_type sw_k= begin<col>(SW);     sw_k+= k;        // column k
+		citer_type it2= begin<all_cit>(sw_k);                  // SW[0][k]
+
+		for(; it1 != it1end; ++it1, ++it2)
+		    *it1 = *it1 - d * *it2;
+	    }
+    }
+
+
+    template < typename MatrixNE, typename MatrixNW, typename MatrixSW >
+    void schur_update_base(MatrixNE & NE, const MatrixNW & NW, const MatrixSW & SW)
+    {
+	vampir_trace<5004> tracer;
+	typedef typename  Collection<MatrixNE>::size_type    size_type;
+
+	using namespace glas::tag; using mtl::traits::range_generator;
+	typedef tag::iter::all        all_it;
+	typedef tag::const_iter::all  all_cit;
+
+        typedef typename range_generator<col, MatrixSW>::type       ccur_type;
+        typedef typename range_generator<all_cit, ccur_type>::type    citer_type;
+
+	typedef typename range_generator<row, MatrixNE>::type       rcur_type;
+	typedef typename range_generator<all_it, rcur_type>::type     riter_type;
+
+	for (size_type k = 0; k < NW.num_cols (); k++)
+	    for (size_type i = 0; i < NE.num_rows (); i++) {
+		typename MatrixNW::value_type d = NW[i][k];
+#if 0
+		rcur_type ne_i= begin<row>(NE);       ne_i+= i;      // row i
+		riter_type it1= begin<all_it>(ne_i);                   // NE[i][0]
+		riter_type it1end= end<all_it>(ne_i);                  // NE[i][num_col]
+
+		ccur_type sw_k= begin<col>(SW);     sw_k+= k;        // column k
+		citer_type it2= begin<all_cit>(sw_k);                  // SW[0][k]
+#endif
+		for (size_type j = 0; j < NE.num_cols (); j++)
+		    NE[i][j] -= d * SW[j][k];
+	    }
+    }
+
+
+    // ======================
+    // Corresponding functors
+    // ======================
+
+    struct cholesky_base_t
+    {
+	template < typename Matrix >
+	void operator() (Matrix & matrix)
+	{
+	    cholesky_base(matrix);
+	}
+    };
+
+    struct tri_solve_base_t
+    {
+	template < typename MatrixSW, typename MatrixNW >
+	void operator() (MatrixSW & SW, const MatrixNW & NW)
+	{
+	    tri_solve_base(SW, NW);
+	}
+    };
+
+    struct tri_schur_base_t
+    {
+	template < typename MatrixSE, typename MatrixSW >
+	void operator() (MatrixSE & SE, const MatrixSW & SW)
+	{
+	    tri_schur_base(SE, SW);
+	}
+    };
+
+    struct schur_update_base_t
+    {
+	template < typename MatrixNE, typename MatrixNW, typename MatrixSW >
+	void operator() (MatrixNE & NE, const MatrixNW & NW, const MatrixSW & SW)
+	{
+	    schur_update_base(NE, NW, SW);
+	}
+    };
+
+} // namespace with_iterator
+#endif
+
+
+// ==================================
+// Functor types for Cholesky visitor
+// ==================================
+
+
+template <typename BaseTest, typename CholeskyBase, typename TriSolveBase, typename TriSchur, typename SchurUpdate>
+struct recursive_cholesky_visitor_t
+{
+    typedef  BaseTest                   base_test;
+
+    template < typename Recursator >
+    bool is_base(const Recursator& recursator) const
+    {
+	return base_test()(recursator);
+    }
+
+    template < typename Matrix >
+    void cholesky_base(Matrix & matrix) const
+    {
+	CholeskyBase()(matrix);
+    }
+
+    template < typename MatrixSW, typename MatrixNW >
+    void tri_solve_base(MatrixSW & SW, const MatrixNW & NW) const
+    {
+	TriSolveBase()(SW, NW);
+    }
+
+    template < typename MatrixSE, typename MatrixSW >
+    void tri_schur_base(MatrixSE & SE, const MatrixSW & SW) const
+    {
+	TriSchur()(SE, SW);
+    }
+
+    template < typename MatrixNE, typename MatrixNW, typename MatrixSW >
+    void schur_update_base(MatrixNE & NE, const MatrixNW & NW, const MatrixSW & SW) const
+    {
+	SchurUpdate()(NE, NW, SW);
+    }
+};
+
+
+namespace detail {
+
+    // Compute schur update with external multiplication; must have Assign == minus_mult_assign_t !!!
+    template <typename MatrixMult>
+    struct mult_schur_update_t
+    {
+	template < typename MatrixNE, typename MatrixNW, typename MatrixSW >
+	void operator()(MatrixNE & NE, const MatrixNW & NW, const MatrixSW & SW)
+	{
+	    transposed_view<MatrixSW> trans_sw(const_cast<MatrixSW&>(SW));
+	    MatrixMult()(NW, trans_sw, NE);
+	}
+    };
+
+} // detail
+
+
+namespace with_bracket {
+    typedef recursive_cholesky_visitor_t<recursion::bound_test_static<64>, cholesky_base_t, tri_solve_base_t,
+					 tri_schur_base_t, schur_update_base_t >
+               recursive_cholesky_base_visitor_t;
+}
+
+#if 0
+namespace with_iterator {
+    typedef recursive_cholesky_visitor_t<recursion::bound_test_static<64>,
+					 cholesky_base_t, tri_solve_base_t, tri_schur_base_t, schur_update_base_t>
+               recursive_cholesky_base_visitor_t;
+}
+#endif
+
+typedef with_bracket::recursive_cholesky_base_visitor_t                    recursive_cholesky_default_visitor_t;
+
+
+
+
+
+
+namespace with_recursator {
+
+    template <typename Recursator, typename Visitor>
+    void schur_update(Recursator E, Recursator W, Recursator N, Visitor vis)
+    {
+	vampir_trace<5005> tracer;
+	using namespace recursion;
+
+	if (E.is_empty() || W.is_empty() || N.is_empty())
+	    return;
+
+	if (vis.is_base(E)) {
+	    typedef typename Visitor::base_test  base_test;
+	    typedef typename base_case_matrix<typename Recursator::matrix_type, base_test>::type matrix_type;
+
+	    matrix_type  base_E(base_case_cast<base_test>(E.get_value())),
+		         base_W(base_case_cast<base_test>(W.get_value())),
+		         base_N(base_case_cast<base_test>(N.get_value()));
+	    vis.schur_update_base(base_E, base_W, base_N);
+	} else{
+	    schur_update(     E.north_east(),W.north_west()     ,N.south_west()     , vis);
+	    schur_update(     E.north_east(),     W.north_east(),     N.south_east(), vis);
+
+	    schur_update(E.north_west()     ,     W.north_east(),     N.north_east(), vis);
+	    schur_update(E.north_west()     ,W.north_west()     ,N.north_west()     , vis);
+
+	    schur_update(E.south_west()     ,W.south_west()     ,N.north_west()     , vis);
+	    schur_update(E.south_west()     ,     W.south_east(),     N.north_east(), vis);
+
+	    schur_update(     E.south_east(),     W.south_east(),     N.south_east(), vis);
+	    schur_update(     E.south_east(),W.south_west()     ,N.south_west()     , vis);
+	}
+    }
+
+
+    template <typename Recursator, typename Visitor>
+    void tri_solve(Recursator S, Recursator N, Visitor vis)
+    {
+	using namespace recursion;
+	vampir_trace<5006> tracer;
+
+        if (S.is_empty())
+	    return;
+
+        if (vis.is_base(S)) {
+	    typedef typename Visitor::base_test  base_test;
+	    typedef typename base_case_matrix<typename Recursator::matrix_type, base_test>::type matrix_type;
+
+	    matrix_type  base_S(base_case_cast<base_test>(S.get_value())),
+		         base_N(base_case_cast<base_test>(N.get_value()));
+
+	    vis.tri_solve_base(base_S, base_N);
+        } else{
+
+	    tri_solve(S.north_west()     ,N.north_west(), vis);
+	    schur_update(  S.north_east(),S.north_west()     ,N.south_west(), vis);
+	    tri_solve(     S.north_east(),     N.south_east(), vis);
+	    tri_solve(S.south_west()     ,N.north_west()     , vis);
+	    schur_update(  S.south_east(),S.south_west()     ,N.south_west(), vis);
+	    tri_solve(     S.south_east(),     N.south_east(), vis);
+	}
+    }
+
+
+    template <typename Recursator, typename Visitor>
+    void tri_schur(Recursator E, Recursator W, Visitor vis)
+    {
+	using namespace recursion;
+	vampir_trace<5007> tracer;
+
+        if (E.is_empty() || W.is_empty())
+	    return;
+
+        if (vis.is_base(W)) {
+	    typedef typename Visitor::base_test  base_test;
+	    typedef typename base_case_matrix<typename Recursator::matrix_type, base_test>::type matrix_type;
+
+	    matrix_type  base_E(base_case_cast<base_test>(E.get_value())),
+               		 base_W(base_case_cast<base_test>(W.get_value()));
+	    vis.tri_schur_base(base_E, base_W);
+        } else{
+
+	    schur_update(E.south_west(),     W.south_west(),    W.north_west(), vis);
+	    schur_update(E.south_west(),     W.south_east(),    W.north_east(), vis);
+	    tri_schur(   E.south_east()     ,     W.south_east(), vis);
+	    tri_schur(   E.south_east()     ,W.south_west()     , vis);
+	    tri_schur(        E.north_west(),     W.north_east(), vis);
+	    tri_schur(        E.north_west(),W.north_west()     , vis);
+        }
+    }
+
+
+    template <typename Recursator, typename Visitor>
+    void cholesky(Recursator recursator, Visitor vis)
+    {
+	using namespace recursion;
+	vampir_trace<5008> tracer;
+
+        if (recursator.is_empty())
+	    return;
+
+        if (vis.is_base (recursator)){
+	    typedef typename Visitor::base_test  base_test;
+	    typedef typename base_case_matrix<typename Recursator::matrix_type, base_test>::type matrix_type;
+
+	    matrix_type  base_matrix(base_case_cast<base_test>(recursator.get_value()));
+	    vis.cholesky_base (base_matrix);
+        } else {
+	    cholesky(recursator.north_west(), vis);
+	    tri_solve(    recursator.south_west(),       recursator.north_west(), vis);
+	    tri_schur(    recursator.south_east(), recursator.south_west(), vis);
+	    cholesky(     recursator.south_east(), vis);
+        }
+    }
+
+} // namespace with_recursator
+
+
+
+template <typename Backup= with_bracket::cholesky_base_t>
+struct recursive_cholesky_t
+{
+    template <typename Matrix>
+    void operator()(Matrix& matrix)
+    {
+	(*this)(matrix, recursive_cholesky_default_visitor_t());
+    }
+
+    template <typename Matrix, typename Visitor>
+    void operator()(Matrix& matrix, Visitor vis)
+    {
+	apply(matrix, vis, typename mtl::traits::category<Matrix>::type());
+    }
+
+private:
+    // If the matrix is not sub-divisible then take backup function
+    template <typename Matrix, typename Visitor>
+    void apply(Matrix& matrix, Visitor, tag::universe)
+    {
+	Backup()(matrix);
+    }
+
+    // Only if matrix is sub-divisible, otherwise backup
+    template <typename Matrix, typename Visitor>
+    void apply(Matrix& matrix, Visitor vis, tag::qsub_divisible)
+    {
+	matrix::recursator<Matrix>  recursator(matrix);
+	with_recursator::cholesky(recursator, vis);
+    }
+};
+
+
+template <typename Matrix, typename Visitor>
+inline void recursive_cholesky(Matrix& matrix, Visitor vis)
+{
+    recursive_cholesky_t<>()(matrix, vis);
+}
+
+template <typename Matrix>
+inline void recursive_cholesky(Matrix& matrix)
+{
+    recursive_cholesky(matrix, recursive_cholesky_default_visitor_t());
+}
+
+
+
+
+
+template <typename Matrix>
+void fill_matrix_for_cholesky(Matrix& matrix)
+{
+    vampir_trace<5008> tracer;
+    typedef typename Collection<Matrix>::size_type   size_type;
+    typedef typename Collection<Matrix>::value_type  value_type;
+
+    value_type   x= 1.0;
+    for (size_type i= 0; i < num_rows(matrix); i++)
+       for (size_type j= 0; j <= i; j++)
+	   if (i != j) {
+	       matrix[i][j]= x; matrix[j][i]= x;
+	       x+= 1.0;
+	   }
+
+    for (size_type i= 0; i < num_rows(matrix); i++) {
+	value_type rowsum= 0.0;
+	for (size_type j=0; j<matrix.num_cols(); j++)
+	    if (i != j)
+		rowsum += matrix[i][j];
+	matrix[i][i]= rowsum * 2;
+    }
+}
+
+
+}} // namespace mtl::matrix
+
+
+
+
+#endif // MTL_CHOLESKY_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/clone.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_CLONE_INCLUDE
+#define MTL_CLONE_INCLUDE
+
+#include <boost/utility/enable_if.hpp>
+#include <boost/mpl/bool.hpp>
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace mtl {
+
+template<typename T>
+struct is_clonable : boost::mpl::false_
+{ };
+
+/// Move-semantics-related anti-dot: always copy in constructor.
+/** Some collections have referring semantics in copy constructors, e.g. sub-matrices.
+    That means
+    \code
+        Matrix B= sub_matrix(A, ...);
+    \endcode
+    creates a sub-matrix of A in B. As a consequence, changes in B modify A and vice versa
+    (unless it's outside the sub-matrix).
+    In contrast, clone forces the copy semantics
+    \code
+        Matrix B= clone(sub_matrix(A, ...));
+    \endcode
+    B now contains the values of A's sub-matrix but is an independent matrix.
+    Modifications to either A or B have no effect to each other.
+    Requires that type T is declared clonable in terms of
+    \code
+        is_clonable<T> : boost::mpl::true_ {};
+    \endcode
+**/
+template <typename T>
+typename boost::enable_if<is_clonable<T>, T>::type
+clone(const T& x)
+{
+    vampir_trace<3004> tracer;
+    // std::cout << "Cloning clone function.\n";
+    return T(x, clone_ctor());
+}
+
+
+template <typename T>
+typename boost::disable_if<is_clonable<T>, T>::type
+clone(const T& x)
+{
+    // std::cout << "Not cloning clone function.\n";
+    return x;
+}
+
+
+} // namespace mtl
+
+#endif // MTL_CLONE_INCLUDE