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

AMDiS/lib/mtl4/boost/numeric/mtl/operation/column_in_matrix.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_COLUMN_IN_MATRIX_INCLUDE
+#define MTL_COLUMN_IN_MATRIX_INCLUDE
+
+#include <boost/mpl/if.hpp>
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/vector/parameter.hpp>
+#include <boost/numeric/mtl/utility/irange.hpp>
+
+namespace mtl {
+
+
+/// Type of column in matrix as vector if accessible
+template <typename Matrix>
+struct ColumnInMatrix
+{
+    static const bool exists= false;
+};
+
+template <typename Value, typename Parameters>
+struct ColumnInMatrix<mtl::matrix::dense2D<Value, Parameters> >
+{
+    typedef mtl::matrix::dense2D<Value, Parameters>         ref_type;
+    typedef typename ref_type::size_type       size_type;
+    typedef typename ref_type::value_type      value_type;
+
+    static const bool aligned= !boost::is_same<typename Parameters::orientation, row_major>::value;
+    static const bool exists= true;
+
+    typedef typename boost::mpl::if_c<
+	aligned
+      , vector::dense_vector<Value, vector::parameters<> >
+      , vector::strided_vector_ref<Value, vector::parameters<> >
+    >::type type;
+
+    static inline type apply(ref_type& A, const irange& row_range, size_type col)
+    {
+	return dispatch<type, ref_type>(A, row_range, col, boost::mpl::bool_<aligned>());
+    }
+
+    template <typename Ref>
+    static inline size_type vector_size(const Ref& A, const irange& row_range)
+    {
+	using std::min;
+	size_type finish= min(row_range.finish(), num_rows(A));
+	return row_range.start() < finish ? finish - row_range.start() : 0;
+    }
+
+    template <typename Return, typename Ref>
+    static inline Return dispatch(Ref& A, const irange& row_range, size_type col, boost::mpl::true_)
+    {
+	return Return(vector_size(A, row_range), const_cast<value_type*>(&A[row_range.start()][col])); // TODO make work without const cast
+    }
+
+    template <typename Return, typename Ref>
+    static inline Return dispatch(Ref& A, const irange& row_range, size_type col, boost::mpl::false_)
+    {
+	return Return(vector_size(A, row_range), &A[row_range.start()][col], num_cols(A));
+    }
+};
+
+template <typename Value, typename Parameters>
+struct ColumnInMatrix<const mtl::matrix::dense2D<Value, Parameters> >
+{
+    typedef mtl::matrix::dense2D<Value, Parameters> const   ref_type;
+    typedef mtl::matrix::dense2D<Value, Parameters>         ref2_type;
+    typedef typename ref2_type::size_type      size_type;
+
+    static const bool aligned= !boost::is_same<typename Parameters::orientation, row_major>::value;
+    static const bool exists= true;
+
+    typedef typename boost::mpl::if_c<
+	aligned
+      , vector::dense_vector<Value, vector::parameters<> > // TODO needs constification !!!
+      , vector::strided_vector_ref<const Value, vector::parameters<> >
+    >::type type;
+
+    static inline type apply(ref_type& A, const irange& row_range, size_type col)
+    {
+	return ColumnInMatrix<ref2_type>::template dispatch<type, ref_type>(A, row_range, col, boost::mpl::bool_<aligned>());
+    }
+};
+
+} // namespace mtl
+
+#endif // MTL_COLUMN_IN_MATRIX_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/compute_factors.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_COMPUTE_FACTORS_INCLUDE
+#define MTL_COMPUTE_FACTORS_INCLUDE
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+
+namespace mtl { namespace operation {
+
+// Default is to just refer to the expression
+template <typename Result, typename Expr>
+struct compute_one_factor
+{
+    typedef const Expr&           type;
+    typedef const Expr&           const_reference;
+
+    compute_one_factor(type src) : value(src) {}
+
+    type value;
+};
+
+template <typename Result, typename E1, typename E2>
+struct compute_one_factor<Result, matrix::mat_mat_times_expr<E1, E2> >
+{
+    typedef Result                type;
+    typedef const Result&         const_reference;
+
+    compute_one_factor(const matrix::mat_mat_times_expr<E1, E2>& src)
+	: value(src.first * src.second) {}
+
+    type value;
+};
+
+template <typename Result, typename E1, typename E2>
+struct compute_one_factor<Result, matrix::mat_mat_ele_times_expr<E1, E2> >
+{
+    typedef Result                type;
+    typedef const Result&         const_reference;
+
+    compute_one_factor(const matrix::mat_mat_ele_times_expr<E1, E2>& src)
+	: value(ele_prod(src.first, src.second)) {}
+
+    type value;
+};
+
+
+// Only defined for matrix::mat_mat_times_expr and matrix::mat_mat_ele_times_expr
+template <typename Result, typename Expr>
+struct compute_factors {};
+
+template <typename Result, typename E1, typename E2>
+struct compute_factors<Result, matrix::mat_mat_times_expr<E1, E2> >
+{
+    compute_factors(const matrix::mat_mat_times_expr<E1, E2>& src)
+	: first_factor(src.first), second_factor(src.second),
+	  first(first_factor.value), second(second_factor.value)
+    {}
+
+    compute_one_factor<Result, E1> first_factor;
+    compute_one_factor<Result, E2> second_factor;
+
+    typename compute_one_factor<Result, E1>::const_reference first;
+    typename compute_one_factor<Result, E2>::const_reference second;
+};
+
+// First factor is handled implicitly in the evaluation
+template <typename Result, typename E1, typename E2>
+struct compute_factors<Result, matrix::mat_mat_ele_times_expr<E1, E2> >
+{
+    compute_factors(const matrix::mat_mat_ele_times_expr<E1, E2>& src)
+	: first(src.first), second_factor(src.second),
+	  second(second_factor.value)
+    {}
+
+    const E1&                                                first;
+    compute_one_factor<Result, E2>                           second_factor;
+    typename compute_one_factor<Result, E2>::const_reference second;
+};
+
+}} // namespace mtl::operation
+
+#endif // MTL_COMPUTE_FACTORS_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/compute_summand.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__OPERATION_COMPUTE_SUMMAND_INCLUDE
+#define MTL__OPERATION_COMPUTE_SUMMAND_INCLUDE
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/utility/copy_expression_const_ref_container.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace mtl { namespace operation {
+
+/// Compute a summand in an expression
+/** For instance matrix vector products are transformed into mult function calls
+    when assigned to a column vector. Adding such a vector to other vector expressions
+    requires to compute (evaluate) the summand first and then add the resulting
+    vector. The current implementation assumes that the result of the operation
+    can be represented by the multiplied column vector. This is for instance wrong
+    when the matrix is complex and the vector real. To handle this is signifantly
+    more complicated and is planned for the future.
+**/
+template <typename Expr>
+struct compute_summand
+{
+    typedef Expr type;
+    compute_summand(const Expr& expr) : value(expr) {}
+    // value is a const& if Expr is a true vector and a copy if it is an expression
+    typename mtl::traits::copy_expression_const_ref_container<Expr>::type value;
+};
+
+
+/// Specialization for matrix vector products
+template <typename Matrix, typename CVector>
+struct compute_summand< mat_cvec_times_expr<Matrix, CVector> >
+{
+    typedef CVector    type;
+
+    compute_summand(const mat_cvec_times_expr<Matrix, CVector>& expr)
+      : value(num_rows(expr.first))
+#ifndef NDEBUG // might be helpful, e.g. for generating AST
+      , first(expr.first), second(expr.second)
+#endif
+    {
+	vampir_trace<3005> tracer;
+	value= expr.first * expr.second;
+    }
+
+    CVector value;
+#ifndef NDEBUG
+    const Matrix&  first;
+    const CVector& second;
+#endif
+};
+
+/// Specialization for matrix vector products
+template <typename Matrix, typename CVector>
+struct compute_summand< vector::mat_cvec_multiplier<Matrix, CVector> >
+{
+    typedef CVector    type;
+
+    compute_summand(const vector::mat_cvec_multiplier<Matrix, CVector>& expr)
+      : value(num_rows(expr.A))
+    {
+	vampir_trace<3005> tracer;
+	expr.assign_to(value);
+    }
+
+    CVector value;
+};
+
+}} // namespace mtl::operation
+
+#endif // MTL__OPERATION_COMPUTE_SUMMAND_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/conj.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_CONJ_INCLUDE
+#define MTL_CONJ_INCLUDE
+
+#include <complex>
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/utility/enable_if.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/mtl/utility/algebraic_category.hpp>
+#include <boost/numeric/mtl/utility/is_what.hpp>
+#include <boost/numeric/mtl/utility/view_code.hpp>
+#include <boost/numeric/mtl/utility/viewed_collection.hpp>
+#include <boost/numeric/mtl/utility/compose_view.hpp>
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/mtl/matrix/view_ref.hpp>
+
+namespace mtl {
+
+namespace sfunctor {
+
+    template <typename Value, typename AlgebraicCategory>
+    struct conj_aux
+    {
+	typedef Value result_type;
+
+	static inline result_type apply(const Value& v)
+	{
+	    return v;
+	}
+
+	result_type operator() (const Value& v) const
+	{
+	    return v;
+	}
+    };
+
+
+    template <typename Value, typename AlgebraicCategory>
+    struct conj_aux<std::complex<Value>, AlgebraicCategory>
+    {
+	typedef std::complex<Value> result_type;
+
+	static inline result_type apply(const std::complex<Value>& v)
+	{
+	    return std::conj(v);
+	}
+
+	result_type operator() (const std::complex<Value>& v) const
+	{
+	    return std::conj(v);
+	}
+    };
+
+    // Only declarations here, definitions in matrix::map_view (vector::map_view)
+    template <typename Matrix>
+    struct conj_aux<Matrix, tag::matrix>;
+
+    template <typename Vector>
+    struct conj_aux<Vector, tag::vector>;
+
+    // Short cut for result type
+    template <typename Value>
+    struct conj
+	: public conj_aux<Value, typename mtl::traits::algebraic_category<Value>::type>
+    {};
+
+} // namespace sfunctor
+
+    namespace vector {
+
+	/// Conjugate of an vector
+	template <typename Vector>
+	typename mtl::traits::enable_if_vector<Vector, conj_view<Vector> >::type
+	inline conj(const Vector& v)
+	{
+	    return conj_view<Vector>(v);
+	}
+    }
+
+    namespace matrix {
+
+	namespace detail {
+
+	    template <typename Matrix>
+	    struct conj_trait
+	    {
+		static const unsigned code= mtl::traits::view_toggle_conj<mtl::traits::view_code<Matrix> >::value;
+		typedef typename mtl::traits::compose_view<code, typename mtl::traits::viewed_collection<Matrix>::type>::type type;
+
+		static inline type apply(const Matrix& A)
+		{
+		    return type(view_ref(A));
+		}
+	    };
+
+	}
+
+	/// Conjugate of a matrix
+	template <typename Matrix>
+	typename mtl::traits::lazy_enable_if_matrix<Matrix, detail::conj_trait<Matrix> >::type
+	inline conj(const Matrix& A)
+	{
+	    return detail::conj_trait<Matrix>::apply(A);
+	}
+    }
+
+    namespace scalar {
+
+	// Only scalar values remain here
+	template <typename Value>
+	typename mtl::traits::enable_if_scalar<
+	    Value
+	  , typename sfunctor::conj<Value>::result_type
+	>::type
+	inline conj(const Value& v)
+	{
+	    return mtl::sfunctor::conj<Value>::apply(v);
+	}
+
+	float inline conj(float v) { return v; }
+	double inline conj(double v) { return v; }
+	long double inline conj(long double v) { return v; }
+    }
+
+    /// Conjugate of vector, matrix, or scalar
+    using vector::conj;
+    using matrix::conj;
+    using scalar::conj;
+    // using std::conj;
+
+} // namespace mtl
+
+#endif // MTL_CONJ_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/copy.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_COPY_INCLUDE
+#define MTL_COPY_INCLUDE
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/detail/index.hpp>
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/mtl/utility/flatcat.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/is_row_major.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/range_generator.hpp>
+#include <boost/numeric/mtl/utility/ashape.hpp>
+#include <boost/numeric/mtl/utility/property_map.hpp>
+#include <boost/numeric/mtl/utility/static_assert.hpp>
+#include <boost/numeric/mtl/utility/updater_to_assigner.hpp>
+#include <boost/numeric/mtl/matrix/inserter.hpp>
+#include <boost/numeric/mtl/operation/set_to_zero.hpp>
+#include <boost/numeric/mtl/operation/update.hpp>
+#include <boost/numeric/mtl/operation/print.hpp>
+#include <boost/numeric/mtl/operation/crop.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+#include <boost/type_traits/is_same.hpp>
+#include <boost/utility/enable_if.hpp>
+#include <iostream>
+#include <limits>
+
+namespace mtl {
+
+    namespace detail {
+
+	// Set Destination matrix to zero when source is sparse
+	// (otherwise everything is overwritten anyway)
+	template <typename MatrixDest>
+	inline void zero_with_sparse_src(MatrixDest& dest, tag::flat<tag::sparse>)
+	{
+	    set_to_zero(dest);
+	}
+
+	template <typename MatrixDest>
+	inline void zero_with_sparse_src(MatrixDest&, tag::universe) {}
+
+	// Adapt inserter size to operation
+	template <typename Updater> struct copy_inserter_size {};
+
+	// Specialization for store
+	template <typename Value>
+	struct copy_inserter_size< operations::update_store<Value> >
+	{
+	    template <typename MatrixSrc, typename MatrixDest>
+	    static inline int apply(const MatrixSrc& src, const MatrixDest& dest)
+	    {
+		// std::cout << "nnz = " << src.nnz() << ", dim1 = " << dest.dim1() << "\n";
+		return int(src.nnz() * 1.2 / dest.dim1());
+	    }
+	};
+
+	struct sum_of_sizes
+	{
+	    template <typename MatrixSrc, typename MatrixDest>
+	    static inline int apply(const MatrixSrc& src, const MatrixDest& dest)
+	    {	return int((src.nnz() + dest.nnz()) * 1.2 / dest.dim1()); }
+	};
+
+	// Specialization for plus and minus
+	template <typename Value> struct copy_inserter_size< operations::update_plus<Value> > : sum_of_sizes {};
+	template <typename Value> struct copy_inserter_size< operations::update_minus<Value> > : sum_of_sizes {};
+
+    } // namespace detail
+
+
+    template <typename Updater, typename MatrixSrc, typename MatrixDest>
+    inline void gen_matrix_copy(const MatrixSrc& src, MatrixDest& dest, bool with_reset)
+    {
+	vampir_trace<3002> tracer;
+	MTL_THROW_IF(num_rows(src) != num_rows(dest) || num_cols(src) != num_cols(dest), incompatible_size());
+
+	if (with_reset)
+	    detail::zero_with_sparse_src(dest, traits::sparsity_flatcat<MatrixSrc>());
+
+	typename traits::row<MatrixSrc>::type             row(src);
+	typename traits::col<MatrixSrc>::type             col(src);
+	typename traits::const_value<MatrixSrc>::type     value(src);
+	typedef typename traits::range_generator<tag::major, MatrixSrc>::type  cursor_type;
+
+	// std::cout << "Slot size is " << detail::copy_inserter_size<Updater>::apply(src, dest) << "\n";
+	matrix::inserter<MatrixDest, Updater>   ins(dest, detail::copy_inserter_size<Updater>::apply(src, dest));
+	for (cursor_type cursor = mtl::begin<tag::major>(src), cend = mtl::end<tag::major>(src);
+	     cursor != cend; ++cursor) {
+	    // std::cout << dest << '\n';
+
+	    typedef typename traits::range_generator<tag::nz, cursor_type>::type icursor_type;
+	    for (icursor_type icursor = mtl::begin<tag::nz>(cursor), icend = mtl::end<tag::nz>(cursor);
+		 icursor != icend; ++icursor) {
+		//std::cout << "in " << row(*icursor) << ", " << col(*icursor) << " insert " << value(*icursor) << '\n';
+		ins(row(*icursor), col(*icursor)) << value(*icursor); }
+	}
+    }
+
+    // Specialization for multi_vector
+    template <typename Updater, typename MatrixSrc, typename Vector>
+    inline void gen_matrix_copy(const MatrixSrc& src, mtl::matrix::multi_vector<Vector>& dest, bool)
+    {
+	MTL_THROW_IF(num_rows(src) != num_rows(dest) || num_cols(src) != num_cols(dest), incompatible_size());
+	typedef typename mtl::traits::updater_to_assigner<Updater>::type Assigner;
+
+	for (std::size_t i= 0, n= num_cols(src); i < n; ++i)
+	    Assigner::first_update(dest.vector(i), src.vector(i));
+    }
+
+    namespace {
+#    ifdef __clang__
+#     pragma clang diagnostic ignored "-Wunneeded-internal-declaration"
+#     pragma clang diagnostic ignored "-Wunused-function"
+#    endif
+	inline long inc_wo_over(long i) { return i == std::numeric_limits<long>::max() ? i : i+1; }
+	inline long negate_wo_over(long i) { return i == std::numeric_limits<long>::min() ? std::numeric_limits<long>::max() : -i; }
+    }
+
+    template <typename Updater, typename ValueSrc, typename Para, typename ValueDest>
+    typename boost::enable_if<boost::is_same<Updater, operations::update_store<ValueDest> > >::type
+    inline gen_matrix_copy(const matrix::banded_view<mtl::matrix::compressed2D<ValueSrc, Para> >& src, mtl::matrix::compressed2D<ValueDest, Para>& dest, bool)
+    {
+	vampir_trace<3061> tracer;
+	typedef typename Para::size_type size_type;
+	dest.change_dim(num_rows(src), num_cols(src)); // contains make_empty
+	set_to_zero(dest);
+	const mtl::matrix::compressed2D<ValueSrc, Para>  &sref= src.ref;
+	const std::vector<size_type>        &sstarts= sref.ref_major(), &sindices= sref.ref_minor();
+	long first, last;
+	if (traits::is_row_major<Para>::value) {
+	    first= src.get_begin();
+	    last= src.get_end();
+	} else {
+	    first= inc_wo_over(negate_wo_over(src.get_end()));
+	    last=  inc_wo_over(negate_wo_over(src.get_begin()));
+	}
+
+	long jd= 0, j_end= sstarts[0];
+	for (long i= 0, i_end= src.dim1(), f= first, l= last; i < i_end; ++i) {
+	    dest.ref_major()[i]= jd;
+	    long j= j_end;
+	    j_end= sstarts[i+1];
+	    while (j < j_end && long(sindices[j]) < f) j++;
+	    while (j < j_end && long(sindices[j]) < l) jd++, j++;
+	    f= inc_wo_over(f);
+	    l= inc_wo_over(l);
+	}
+	dest.ref_major()[src.dim1()]= jd;
+	dest.set_nnz(jd); // resizes indices and data
+
+	for (long i= 0, i_end= src.dim1(), jd= 0, j_end= sstarts[0]; i < i_end; ++i) {
+	    dest.ref_major()[i]= jd;
+	    long j= j_end;
+	    j_end= sstarts[i+1];
+	    while (j < j_end && long(sindices[j]) < first) j++;
+	    while (j < j_end && long(sindices[j]) < last) {
+		dest.ref_minor()[jd]= sindices[j];
+		dest.data[jd++]= sref.data[j++];
+	    }
+	    first= inc_wo_over(first);
+	    last= inc_wo_over(last);
+	}
+    }
+
+
+    /// Copy matrix \p src into matrix \p dest
+    template <typename MatrixSrc, typename MatrixDest>
+    inline void matrix_copy(const MatrixSrc& src, MatrixDest& dest)
+    {
+	gen_matrix_copy< operations::update_store<typename MatrixDest::value_type> >(src, dest, true);
+    }
+
+
+    /// Add matrix \p src to matrix \p dest in copy-like style
+    template <typename MatrixSrc, typename MatrixDest>
+    inline void matrix_copy_plus(const MatrixSrc& src, MatrixDest& dest)
+    {
+	gen_matrix_copy< operations::update_plus<typename MatrixDest::value_type> >(src, dest, false);
+    }
+
+    /// Subtract matrix \p src from matrix \p dest in copy-like style
+    template <typename MatrixSrc, typename MatrixDest>
+    inline void matrix_copy_minus(const MatrixSrc& src, MatrixDest& dest)
+    {
+	gen_matrix_copy< operations::update_minus<typename MatrixDest::value_type> >(src, dest, false);
+    }
+
+    /// Multiply matrix \p src element-wise with matrix \p dest in copy-like style
+    template <typename MatrixSrc, typename MatrixDest>
+    inline void matrix_copy_ele_times(const MatrixSrc& src, MatrixDest& dest)
+    {
+	vampir_trace<3001> tracer;
+	MTL_THROW_IF(num_rows(src) != num_rows(dest) || num_cols(src) != num_cols(dest), incompatible_size());
+
+	typename traits::row<MatrixDest>::type             row(dest);
+	typename traits::col<MatrixDest>::type             col(dest);
+	typename traits::value<MatrixDest>::type           value(dest);
+	typedef typename traits::range_generator<tag::major, MatrixDest>::type  cursor_type;
+	typedef typename traits::range_generator<tag::nz, cursor_type>::type icursor_type;
+
+	for (cursor_type cursor = begin<tag::major>(dest), cend = end<tag::major>(dest); cursor != cend; ++cursor)
+	    for (icursor_type icursor = begin<tag::nz>(cursor), icend = end<tag::nz>(cursor); icursor != icend; ++icursor)
+		value(*icursor, value(*icursor) * src[row(*icursor)][col(*icursor)]);
+#if 0   // copy would result in a*0 = a and 0*b = b!!!!
+	gen_matrix_copy< operations::update_times<typename MatrixDest::value_type> >(src, dest, false);
+#endif
+	crop(dest);
+    }
+
+
+    template <typename MatrixSrc, typename MatrixDest>
+    inline void copy(const MatrixSrc& src, tag::flat<tag::matrix>, MatrixDest& dest, tag::flat<tag::matrix>)
+	// inline void copy(const MatrixSrc& src, tag::matrix_expr, MatrixDest& dest, tag::matrix)
+    {
+	return matrix_copy(src, dest);
+    }
+
+
+
+    template <typename Updater, typename VectorSrc, typename VectorDest>
+    inline void gen_vector_copy(const VectorSrc& src, VectorDest& dest, bool with_reset)
+    {
+	// Works only with dense vectors as dest !!!!! (source could be sparse)
+	// Needs vector inserter
+	vampir_trace<2001> tracer;
+
+	MTL_STATIC_ASSERT((boost::is_same<typename ashape::ashape<VectorSrc>::type,
+					  typename ashape::ashape<VectorDest>::type>::value), "Source and target must have the same algebraic shape.");
+
+	MTL_THROW_IF(size(src) != size(dest), incompatible_size());
+
+	if (with_reset)
+	    detail::zero_with_sparse_src(dest, typename traits::category<VectorSrc>::type());
+
+	typename traits::index<VectorSrc>::type           index(src);
+	typename traits::const_value<VectorSrc>::type     value(src);
+
+	typedef typename traits::range_generator<tag::nz, VectorSrc>::type  cursor_type;
+	for (cursor_type cursor = begin<tag::nz>(src), cend = end<tag::nz>(src);
+	     cursor != cend; ++cursor)
+	    Updater()(dest[index(*cursor)], value(*cursor));
+    }
+
+    /// Copy vector \p src into vector \p dest
+    template <typename VectorSrc, typename VectorDest>
+    inline void vector_copy(const VectorSrc& src, VectorDest& dest)
+    {
+	gen_vector_copy< operations::update_store<typename VectorDest::value_type> >(src, dest, true);
+    }
+
+
+    /// Add vector \p src to vector \p dest in copy-like style
+    template <typename VectorSrc, typename VectorDest>
+    inline void vector_copy_plus(const VectorSrc& src, VectorDest& dest)
+    {
+	gen_vector_copy< operations::update_plus<typename VectorDest::value_type> >(src, dest, false);
+    }
+
+    /// Subtract vector \p src from vector \p dest in copy-like style
+    template <typename VectorSrc, typename VectorDest>
+    inline void vector_copy_minus(const VectorSrc& src, VectorDest& dest)
+    {
+	gen_vector_copy< operations::update_minus<typename VectorDest::value_type> >(src, dest, false);
+    }
+
+
+
+    template <typename VectorSrc, typename VectorDest>
+    inline void copy(const VectorSrc& src, tag::flat<tag::vector>, VectorDest& dest, tag::flat<tag::vector>)
+    {
+	return vector_copy(src, dest);
+    }
+
+
+    template <typename CollSrc, typename CollDest>
+    inline void copy(const CollSrc& src, CollDest& dest)
+    {
+	vampir_trace<3003> tracer;
+	return copy(src, traits::flatcat2<CollSrc, tag::matrix, tag::vector>(),
+		    dest, traits::flatcat2<CollDest, tag::matrix, tag::vector>());
+    }
+
+
+} // namespace mtl
+
+#endif // MTL_COPY_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/copysign.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_COPYSIGN_INCLUDE
+#define MTL_COPYSIGN_INCLUDE
+
+#include <complex>
+#include <cmath>
+#include <boost/numeric/mtl/operation/real.hpp>
+#include <boost/numeric/mtl/operation/imag.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace mtl {
+
+namespace sfunctor {
+
+    template <typename Value1, typename Value2>
+    struct copysign
+    {
+	static inline Value1 apply(const Value1& v, const Value2& s)
+	{
+	    using math::zero; using std::abs;
+
+	    Value1 a(abs(v));
+	    return s < zero(s) ? -a : a;
+	}
+    };
+
+    // This specialization is questionable and thus subject to elimination
+    template <typename Value1, typename Value2>
+    struct copysign<std::complex<Value1>, Value2>
+    {
+	static inline Value1 apply(const std::complex<Value1>& v, const Value2& s)
+	{
+	    using mtl::real; using mtl::imag;
+	    return std::complex<Value1>(copysign<Value1, Value2>::apply(real(v), s),
+					copysign<Value1, Value2>::apply(imag(v), s));
+	}
+    };
+#if 0 // ndef _MSC_VER doesn't work on BigRed with gcc 4.2.2 (????) and isn't used anyway
+    template <>
+    struct copysign<float, float>
+    {
+	static inline float apply(float v, float s)
+	{ return ::copysignf(v, s); }
+    };
+
+    template <>
+    struct copysign<double, double>
+    {
+	static inline double apply(double v, double s)
+	{ return ::copysign(v, s); }
+    };
+
+    template <>
+    struct copysign<long double, long double>
+    {
+	static inline long double apply(long double v, long double s)
+	{ return copysignl(v, s); }
+    };
+#endif // _MSC_VER
+}
+
+/// sign of scalars; for complex numbers sign of real part
+template <typename Value1, typename Value2>
+inline Value1 copysign(const Value1& v, const Value2& s)
+{
+    vampir_trace<1001> tracer;
+    return sfunctor::copysign<Value1, Value2>::apply(v, s);
+}
+
+
+
+
+} // namespace mtl
+
+#endif // MTL_COPYSIGN_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/crop.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_CROP_INCLUDE
+#define MTL_CROP_INCLUDE
+
+#include <boost/numeric/mtl/utility/enable_if.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace mtl {
+
+    namespace vector {
+
+	/// Remove all zero entries from a collection
+	/** Does nothing for dense collections **/
+	template <typename T>
+	typename mtl::traits::enable_if_vector<T, T&>::type inline crop(T& x)
+	{
+	    vampir_trace<3006> tracer;
+	    x.crop(); return x;
+	}
+    }
+
+    namespace matrix {
+
+	/// Remove all zero entries from a collection
+	/** Does nothing for dense collections **/
+	template <typename T>
+	typename mtl::traits::enable_if_matrix<T, T&>::type inline crop(T& x)
+	{
+	    vampir_trace<3006> tracer;
+	    x.crop(); return x;
+	}
+    }
+
+    using vector::crop;
+    using matrix::crop;
+
+} // namespace mtl
+
+#endif // MTL_CROP_INCLUDE