AMDiS/lib/mtl4/boost/numeric/mtl/operation/frobenius_norm.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_FROBENIUS_NORM_INCLUDE
+#define MTL_FROBENIUS_NORM_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/is_row_major.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/range_generator.hpp>
+#include <boost/numeric/mtl/utility/property_map.hpp>
+#include <boost/numeric/mtl/operation/max_of_sums.hpp>
+#include <boost/numeric/mtl/operation/squared_abs.hpp>
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+
+namespace mtl { namespace matrix {
+
+/// Frobenius norm, i.e. square root of sum of squares of all entries sqrt(sum_i sum_j(|a[i][j]|^2))
+template <typename Matrix>
+typename RealMagnitude<typename Collection<Matrix>::value_type>::type
+inline frobenius_norm(const Matrix& matrix)
+{
+    vampir_trace<3010> tracer;
+    using std::sqrt; using std::abs; using math::zero;
+    namespace traits = mtl::traits;
+    typename traits::const_value<Matrix>::type     value(matrix);
+
+    typedef typename Collection<Matrix>::value_type   value_type;
+    typedef typename RealMagnitude<value_type>::type  real_type;
+    real_type ref, sum= zero(ref);
+
+    typedef typename traits::range_generator<tag::major, Matrix>::type     cursor_type;
+    typedef typename traits::range_generator<tag::nz, cursor_type>::type   icursor_type;
+
+    for (cursor_type cursor = begin<tag::major>(matrix), cend = end<tag::major>(matrix); cursor != cend; ++cursor)
+	for (icursor_type icursor = begin<tag::nz>(cursor), icend = end<tag::nz>(cursor); icursor != icend; ++icursor)
+	    sum+= squared_abs(value(*icursor));
+    return sqrt(sum);
+}
+
+}} // namespace mtl::matrix
+
+#endif // MTL_FROBENIUS_NORM_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/fuse.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_FUSE_INCLUDE
+#define MTL_FUSE_INCLUDE
+
+#include <boost/utility/enable_if.hpp>
+#include <boost/mpl/and.hpp>
+
+#include <boost/numeric/mtl/utility/is_lazy.hpp>
+#include <boost/numeric/mtl/operation/fused_expr.hpp>
+
+namespace mtl {
+
+/// Fuse two lazy expressions (operator notation of fuse)
+template <typename T, typename U>
+typename boost::enable_if<boost::mpl::and_<traits::is_lazy<T>, traits::is_lazy<U> >, fused_expr<T, U> >::type
+operator||(const T& x, const U& y)
+{
+    return fused_expr<T, U>(const_cast<T&>(x), const_cast<U&>(y));
+}
+
+/// Fuse two lazy expressions
+template <typename T, typename U>
+typename boost::enable_if<boost::mpl::and_<traits::is_lazy<T>, traits::is_lazy<U> >, fused_expr<T, U> >::type
+fuse(const T& x, const U& y)
+{
+    return fused_expr<T, U>(const_cast<T&>(x), const_cast<U&>(y));
+}
+
+
+
+} // namespace mtl
+
+#endif // MTL_FUSE_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/fused_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_FUSED_EXPR_INCLUDE
+#define MTL_FUSED_EXPR_INCLUDE
+
+#include <boost/mpl/bool.hpp>
+#include <boost/mpl/and.hpp>
+#include <boost/numeric/mtl/operation/index_evaluator.hpp>
+#include <boost/numeric/mtl/utility/index_evaluatable.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace mtl {
+
+/// Expression template for fusing other expression
+template <typename T, typename U>
+struct fused_expr
+{
+    typedef boost::mpl::and_<traits::forward_index_evaluatable<T>, traits::forward_index_evaluatable<U> >   forward;
+    typedef boost::mpl::and_<traits::backward_index_evaluatable<T>, traits::backward_index_evaluatable<U> > backward;
+
+    fused_expr(T& first, U& second)
+      : first(first), second(second), delayed_assign(false)
+    {
+	first.delay_assign(); second.delay_assign();
+    }
+
+    ~fused_expr() { if (!delayed_assign) eval(forward(), backward()); }
+
+    void delay_assign() const { delayed_assign= true; }
+
+    template <typename TT, typename UU>
+    void forward_eval_loop(const TT& const_first_eval, const UU& const_second_eval, boost::mpl::false_)
+    {
+	vampir_trace<6001> tracer;
+	// hope there is a more elegant way; copying the arguments causes errors due to double destructor evaluation
+	TT& first_eval= const_cast<TT&>(const_first_eval);
+	UU& second_eval= const_cast<UU&>(const_second_eval);
+	MTL_DEBUG_THROW_IF(mtl::vector::size(first_eval) != mtl::vector::size(second_eval), incompatible_size());
+
+	for (std::size_t i= 0, s= size(first_eval); i < s; i++) {
+	    first_eval(i); second_eval(i);
+	}
+    }
+
+    template <typename TT, typename UU>
+    void forward_eval_loop(const TT& const_first_eval, const UU& const_second_eval, boost::mpl::true_)
+    {
+	vampir_trace<6002> tracer;
+	// hope there is a more elegant way; copying the arguments causes errors due to double destructor evaluation
+	TT& first_eval= const_cast<TT&>(const_first_eval);
+	UU& second_eval= const_cast<UU&>(const_second_eval);
+	MTL_DEBUG_THROW_IF(mtl::vector::size(first_eval) != mtl::vector::size(second_eval), incompatible_size());
+
+	const std::size_t s= size(first_eval), sb= s >> 2 << 2;
+
+	for (std::size_t i= 0; i < sb; i+= 4) {
+	    first_eval.template at<0>(i); second_eval.template at<0>(i);
+	    first_eval.template at<1>(i); second_eval.template at<1>(i);
+	    first_eval.template at<2>(i); second_eval.template at<2>(i);
+	    first_eval.template at<3>(i); second_eval.template at<3>(i);
+	}
+
+	for (std::size_t i= sb; i < s; i++) {
+	    first_eval(i); second_eval(i);
+	}
+    }
+
+    // Forward evaluation dominates backward
+    template <bool B2>
+    void eval(boost::mpl::true_, boost::mpl::bool_<B2>)
+    {
+#ifdef MTL_LAZY_LOOP_WO_UNROLL
+	typedef boost::mpl::false_                                                                              to_unroll;
+#else
+	typedef boost::mpl::and_<traits::unrolled_index_evaluatable<T>, traits::unrolled_index_evaluatable<U> > to_unroll;
+#endif
+	// Currently lazy evaluation is only available on vector expressions, might change in the future
+	// std::cout << "Forward evaluation\n";
+	forward_eval_loop(index_evaluator(first), index_evaluator(second), to_unroll());
+    }
+
+    template <typename TT, typename UU>
+    void backward_eval_loop(const TT& const_first_eval, const UU& const_second_eval, boost::mpl::false_)
+    {
+	vampir_trace<6003> tracer;
+	// hope there is a more elegant way; copying the arguments causes errors due to double destructor evaluation
+	TT& first_eval= const_cast<TT&>(const_first_eval);
+	UU& second_eval= const_cast<UU&>(const_second_eval);
+	MTL_DEBUG_THROW_IF(mtl::vector::size(first_eval) != mtl::vector::size(second_eval), incompatible_size());
+
+	for (std::size_t i= size(first_eval); i-- > 0; ) {
+	    // std::cout << "i is " << i << "\n";
+	    first_eval(i); second_eval(i);
+	}
+    }
+
+    template <typename TT, typename UU>
+    void backward_eval_loop(const TT& const_first_eval, const UU& const_second_eval, boost::mpl::true_)
+    {
+	vampir_trace<6004> tracer;
+	// hope there is a more elegant way; copying the arguments causes errors due to double destructor evaluation
+	TT& first_eval= const_cast<TT&>(const_first_eval);
+	UU& second_eval= const_cast<UU&>(const_second_eval);
+	MTL_DEBUG_THROW_IF(mtl::vector::size(first_eval) != mtl::vector::size(second_eval), incompatible_size());
+
+	std::size_t s= size(first_eval), i= s-1, m= s % 4;
+	for (; m; m--) {
+	    // std::cout << "i is " << i << "\n";
+	    first_eval(i); second_eval(i--);
+	}
+	for(long j= i - 3; j >= 0; j-= 4) {
+	    // std::cout << "i is " << j+3 << ".." << j << "\n";
+	    first_eval.template at<3>(j); second_eval.template at<3>(j);
+	    first_eval.template at<2>(j); second_eval.template at<2>(j);
+	    first_eval.template at<1>(j); second_eval.template at<1>(j);
+	    first_eval.template at<0>(j); second_eval.template at<0>(j);
+	}
+    }
+
+    // Backward evaluation if forward isn't possible
+    void eval(boost::mpl::false_, boost::mpl::true_)
+    {
+	// std::cout << "Backward evaluation\n";
+	backward_eval_loop(index_evaluator(first), index_evaluator(second), boost::mpl::true_());
+    }
+
+    // Sequential evaluation
+    void eval(boost::mpl::false_, boost::mpl::false_)
+    {
+	// std::cout << "Non-fused evaluation\n";
+	evaluate_lazy(first); evaluate_lazy(second);
+    }
+
+    T& first;
+    U& second;
+    mutable bool                        delayed_assign;
+};
+
+
+} // namespace mtl
+
+#endif // MTL_FUSED_EXPR_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/fused_index_evaluator.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_FUSED_INDEX_EVALUATOR_INCLUDE
+#define MTL_FUSED_INDEX_EVALUATOR_INCLUDE
+
+#include <boost/numeric/mtl/utility/index_evaluator.hpp>
+
+namespace mtl { namespace vector {
+
+template <typename T, typename U>
+struct fused_index_evaluator
+{
+    fused_index_evaluator(T& first, U& second)
+      : first(index_evaluator(first)), second(index_evaluator(second)) {}
+
+    template <unsigned Offset>
+    void at(std::size_t i)
+    {
+	first.template at<Offset>(i);
+	second.template at<Offset>(i);
+    }
+
+    void operator() (std::size_t i) { at<0>(i); }
+    void operator[] (std::size_t i) { at<0>(i); }
+
+    typename traits::index_evaluator<T>::type first;
+    typename traits::index_evaluator<U>::type second;
+};
+
+template <typename T, typename U>
+inline size_t size(const fused_index_evaluator<T, U>& expr) { return size(expr.first); }
+
+}} // namespace mtl::vector
+
+#endif // MTL_FUSED_INDEX_EVALUATOR_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/givens.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.
+
+// With contributions from Cornelius Steinhardt
+
+#ifndef MTL_MATRIX_GIVENS_INCLUDE
+#define MTL_MATRIX_GIVENS_INCLUDE
+
+#include <cmath>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/irange.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/operation/householder.hpp>
+#include <boost/numeric/mtl/operation/rank_one_update.hpp>
+#include <boost/numeric/mtl/operation/trans.hpp>
+
+namespace mtl { namespace matrix {
+
+/// Given's transformator
+/** Requires Hessenberg form, i.e. for transformations near the diagonal.
+    For general form use qr_givens.
+    \sa qr_givens. **/
+template <typename Matrix>
+class givens
+{
+    typedef typename Collection<Matrix>::value_type   value_type;
+    typedef typename Collection<Matrix>::size_type    size_type;
+
+  public:
+    /// Re-set the rotation parameters \p a and \p b
+    void set_rotation(value_type a, value_type b)
+    {
+	using std::abs;
+	value_type zero= math::zero(a), one= math::one(b), t;
+
+	if ( b == zero ) {
+	    c= one; s= zero;
+	} else if ( abs(b) > abs(a) ) {
+	    t= -a / b;
+	    s= one / sqrt(one + t*t);
+	    c= s * t;
+	} else {
+	    t= -b / a;
+	    c= one / sqrt(one + t*t);
+	    s= c * t;
+	}
+	G= c, s,
+	  -s, c;
+    }
+
+
+    /// Constructor takes %matrix \p H to be transformed and the rotation parameters \p a and \p b
+    givens(Matrix& H, value_type a, value_type b) : H(H), G(2, 2)
+    {	set_rotation(a, b);    }
+
+    /// Given's transformation of \p H with \p G regarding column \p k
+    Matrix& trafo(const Matrix& G, size_type k)
+    {
+	    irange r(k,k+2);
+	    // trans(H[r][ind])*= G; H[ind][r]*= G; // most compact form but does not work yet
+
+	    Matrix col_block(H[r][iall]), col_perm(trans(G) * col_block);
+	    H[r][iall]= col_perm;
+	    Matrix row_perm(H[iall][r] * G);
+	    H[iall][r]= row_perm;
+
+	    return H;
+    }
+
+    /// Given's transformation of \p H regarding column \p k
+    Matrix& trafo(size_type k)
+    {
+	return trafo(G, k);
+    }
+
+  private:
+    Matrix&    H, G;
+    value_type c, s;
+};
+
+}// namespace matrix
+
+
+namespace vector {
+
+/// Given's transformator on %vector (swap a*line(k) with b*line(k+1) )
+template <typename Vector>
+class givens
+{
+    typedef typename Collection<Vector>::value_type   value_type;
+    typedef typename Collection<Vector>::size_type    size_type;
+
+  public:
+    /// Constructor takes %vector \p v to be transformed and the rotation parameters \p a and \p b
+    givens(Vector& v, value_type a, value_type b) : v(v), a(a), b(b)
+    {  }
+
+    /// Given's transformation of \p v with \p a and \p b regarding column \p k
+    Vector& trafo(size_type k)
+    {
+	    value_type w1(0), w2(0);
+	    w1= a*v[k] - b*v[k+1]; //given's rotation on solution
+            w2= b*v[k] + a*v[k+1]; //rotation on vector
+            v[k]= w1;
+            v[k+1]= w2;
+
+	    return v;
+    }
+
+  private:
+    Vector&    v;
+    value_type a, b;
+};
+
+}// namespace vector
+
+
+} // namespace mtl
+
+#endif // MTL_MATRIX_GIVENS_INCLUDE

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

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_HERMITIAN_INCLUDE
+#define MTL_HERMITIAN_INCLUDE
+
+#include <boost/numeric/mtl/matrix/hermitian_view.hpp>
+#include <boost/numeric/mtl/matrix/view_ref.hpp>
+#include <boost/numeric/mtl/utility/enable_if.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>
+
+namespace mtl {
+
+    // vector version to be done
+
+    namespace matrix {
+
+	namespace detail {
+
+	    template <typename Matrix>
+	    struct hermitian
+	    {
+		static const unsigned code= mtl::traits::view_toggle_hermitian<mtl::traits::view_code<Matrix> >::value;
+		typedef typename mtl::traits::compose_view<code, typename mtl::traits::viewed_collection<Matrix>::type>::type result_type;
+
+		static inline result_type apply(const Matrix& A)
+		{
+		    return result_type(view_ref(A));
+		}
+	    };
+
+	} // namespace detail
+
+	/// Return hermitian of matrix A
+	template <typename Matrix>
+	typename mtl::traits::enable_if_matrix<Matrix, typename detail::hermitian<Matrix>::result_type >::type
+	inline hermitian(const Matrix& A)
+	{
+	    return detail::hermitian<Matrix>::apply(A);
+	}
+    }
+
+    using matrix::hermitian;
+
+} // namespace mtl
+
+#endif // MTL_HERMITIAN_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/hessenberg.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.
+
+// With contributions from Cornelius Steinhardt
+
+#ifndef MTL_MATRIX_HESSENBERG_INCLUDE
+#define MTL_MATRIX_HESSENBERG_INCLUDE
+
+#include <cmath>
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/mtl/vector/parameter.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/irange.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/concept/magnitude.hpp>
+#include <boost/numeric/mtl/operation/householder.hpp>
+#include <boost/numeric/mtl/operation/rank_one_update.hpp>
+#include <boost/numeric/mtl/operation/trans.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+
+namespace mtl { namespace matrix {
+
+/// Hessenberg-Factorization of matrix A with householder-vectors
+/** Return Hessenberg matrix and tril(B,-2) are Householder-vectors **/
+template <typename Matrix>
+Matrix inline hessenberg_factors(const Matrix& A)
+{
+    vampir_trace<5014> tracer;
+    if (num_rows(A) < 3)
+	return A;
+
+    using mtl::imax;
+    typedef typename Collection<Matrix>::value_type   value_type;
+    typedef typename Collection<Matrix>::size_type    size_type;
+    size_type        ncols = num_cols(A), nrows = num_rows(A);
+    value_type       zero= math::zero(A[0][0]), beta;
+    Matrix           B(clone(A));
+
+    for(size_type i= 0; i < ncols-2; i++){
+	// mtl::vector::dense_vector<value_type>  v(B[irange(i+1, imax)][i]);
+		mtl::vector::dense_vector<value_type, vector::parameters<> >  v(nrows-i-1), w(nrows);
+        for (size_type j = 0; j < size(v); j++)
+            v[j]= B[j+i+1][i];
+        beta= householder(v).second;
+        v= householder(v).first;
+;
+	if( beta != zero){
+            w= beta * B[irange(0,imax)][irange(i+1,imax)] * v;
+	    //rank_one_update(B[irange(0,imax)][irange(i+1,imax)],-w,v);
+            for(size_type row = 0; row < nrows; row++){
+                for(size_type col = i+1; col < ncols; col++){
+                    B[row][col] -= w[row] * v[col-i-1];
+                }
+            }
+            //vector*Matrix
+            for(size_type k=0; k < size(w); k++){
+                w[k]= zero;
+                for(size_type j = 0; j < size(v); j++){
+                    w[k] += beta * v[j] * B[j+i+1][k];
+                }
+            }
+            //rank_one_update(A[irange(i+1,imax)][irange(0,imax)],-v,w);
+            for(size_type row = i+1; row < nrows; row++){
+                for(size_type col = 0; col < ncols; col++){
+                    B[row][col] -= v[row-i-1] * w[col];
+                }
+            }
+	    // B[irange(i+2, imax)][i]= v[irange(1, nrows-i-1)];
+            for(size_type row = i+2; row < nrows; row++){
+                B[row][i] = v[row-i-1];
+            }
+        }
+    }
+    return B;
+}
+
+
+/// Extract Householder vectors from Hessenberg factorization H of some A
+template <typename Matrix>
+Matrix inline extract_householder_hessenberg(const Matrix& H)
+{
+    vampir_trace<5015> tracer;
+    return Matrix(tril(H, -2));
+}
+
+/// Compute Householder vectors from Hessenberg factorization of A
+template <typename Matrix>
+Matrix inline householder_hessenberg(const Matrix& A)
+{
+    vampir_trace<5016> tracer;
+    return Matrix(tril(hessenberg_factors(A), -2));
+}
+
+
+/// Extract Hessenberg form from factorization H of some A
+template <typename Matrix>
+Matrix inline extract_hessenberg(const Matrix& H)
+{
+    vampir_trace<5017> tracer;
+    return Matrix(triu(H, -1));
+}
+
+/// Hessenberg form of A
+template <typename Matrix>
+Matrix inline hessenberg(const Matrix& A)
+{
+    vampir_trace<5018> tracer;
+    // return triu(hessenberg_factors(A), -1);
+    MTL_THROW_IF(num_rows(A) < 3, matrix_too_small());
+
+    typedef typename Collection<Matrix>::value_type   value_type;
+    // typedef typename Magnitude<value_type>::type      magnitude_type; // to multiply with 2 not 2+0i
+    typedef typename Collection<Matrix>::size_type    size_type;
+    size_type        ncols = num_cols(A), nrows = num_rows(A);
+    value_type       zero= math::zero(A[0][0]);
+    Matrix           H(nrows,ncols);
+
+    H= hessenberg_factors(A);
+
+    // H= bands(hessenberg_factors(A), -nrows, -1);
+    // set (doubly) strict lower triangle to zero
+    for(size_type row = 2; row < nrows; row++){
+        for(size_type col = 0; col < row-1; col++){
+            H[row][col]= zero;
+        }
+    }
+
+    return H;
+}
+
+
+/// Return Q where Q'*A*Q == hessenberg(A)
+template <typename Matrix>
+Matrix inline hessenberg_q(const Matrix& A)
+{
+    vampir_trace<5013> tracer;
+    using std::abs;
+    typedef typename Collection<Matrix>::value_type   value_type;
+    typedef typename Magnitude<value_type>::type      magnitude_type; // to multiply with 2 not 2+0i
+    typedef typename Collection<Matrix>::size_type    size_type;
+    size_type        ncols = num_cols(A), nrows = num_rows(A), mini;
+    value_type       zero= math::zero(A[0][0]), one= math::one(A[0][0]);
+    const magnitude_type two(2);
+    Matrix           Q(nrows,ncols);
+
+    MTL_THROW_IF(num_rows(A) < 3, matrix_too_small());
+
+    Q= one;
+
+    //Extract Q
+    for(size_type i = 0; i < nrows-2; i++){
+		mtl::vector::dense_vector<value_type, vector::parameters<> >   v(nrows-1), w(nrows);
+        v[0]= one;
+// 	std::cout<< "v=" << v << "\n";
+// 	std::cout<< "w=" << w << "\n";
+        for(size_type k = 1; k < size(v); k++)
+            v[k]= A[nrows-k][i];
+
+	magnitude_type beta= two / abs(dot(v, v)); // abs: x+0i -> x
+        if (beta != two) {
+            //trans(Vector)*Matrix
+            for(size_type k = 0; k < size(w); k++) {
+                w[k]= zero;
+                for(size_type j = 0; j < size(v); j++){
+// 		     std::cout<< "k=" << k << "  j=" << j << "\n";
+                    w[k]+= beta * v[j] * Q[j+i+1][k];
+		}
+            }
+            //rank_one_update(Q[irange(i+1,imax)][irange(0,imax)],-v,w);
+            for(size_type row = i+1; row < nrows; row++)
+                for(size_type col = 0; col < ncols; col++){
+// 		    std::cout<< "row=" << row << "  col=" << col << "\n";
+                    Q[row][col] -= v[row-1] * w[col];
+		}
+        }
+    }
+    return Q;
+}
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_HESSENBERG_INCLUDE
+

AMDiS/lib/mtl4/boost/numeric/mtl/operation/householder.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.
+
+// With contributions from Cornelius Steinhardt
+
+#ifndef MTL_MATRIX_HOUSEHOLDER_INCLUDE
+#define MTL_MATRIX_HOUSEHOLDER_INCLUDE
+
+#include <cmath>
+#include <cassert>
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/vector/parameter.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+
+namespace mtl { namespace vector {
+
+
+/// Computes Householder vector v and scalar b for vector \p y
+/** such that identity_matrix(size(y))-b*v*v' projects the vector y
+    to a positive multiple of the first unit vector. **/
+template <typename Vector>
+std::pair<typename mtl::vector::dense_vector<typename Collection<Vector>::value_type, parameters<> >, typename Collection<Vector>::value_type>
+inline householder(Vector& y)
+{
+    vampir_trace<2004> tracer;
+    assert(size(y) > 0);
+    typedef typename  Collection<Vector>::value_type   value_type;
+    // typedef typename  Collection<Vector>::size_type    size_type;
+    const value_type  zero= math::zero(y[0]), one= math::one(y[0]);
+
+    Vector            v(y);
+    v[0]= one;
+    irange            tail(1, imax);
+    value_type        s( dot(v[tail], v[tail]) ), b, v0;
+
+    //evaluation of v and b
+    if (s == zero)
+        b= zero;
+    else {
+	value_type mu= sqrt(y[0] * y[0] + s);
+	v0= v[0]= y[0] <= zero ? y[0] - mu : -s / (y[0] + mu); // complex < zero????
+	b= 2 * v0 * v0 / (s + v0 * v0);                        // 2* complex???????
+	v/= v0;                                                // normalization of the first entry
+    }
+
+    return std::make_pair(v,b);
+}
+
+/// Computes Householder vector for vector \p y
+/** More stable Householder transformation, also for non-square matrices. **/
+template <typename Vector>
+typename mtl::vector::dense_vector<typename Collection<Vector>::value_type, parameters<> >
+inline householder_s(Vector& y)
+{
+    vampir_trace<2005> tracer;
+    typedef typename  Collection<Vector>::value_type   value_type;
+    const value_type  zero= math::zero(y[0]);
+
+    Vector            u(y);
+    value_type        nu(sqrt( dot(u, u) )), s;
+
+    if (nu != zero){
+	if(u[0] < 0){
+		s= -1;
+	} else {
+		s= 1;
+	}
+	u[0]= u[0] + s * nu;
+	u/= sqrt( dot(u, u) );
+    }
+
+    return u;
+}
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_HOUSEHOLDER_INCLUDE
+

AMDiS/lib/mtl4/boost/numeric/mtl/operation/imag.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_IMAG_INCLUDE
+#define MTL_IMAG_INCLUDE
+
+#include <complex>
+#include <boost/numeric/linear_algebra/identity.hpp>
+
+namespace mtl {
+
+namespace sfunctor {
+
+    template <typename Value>
+    struct imag
+    {
+	typedef Value result_type;
+
+	static inline result_type apply(const Value& v)
+	{
+	    using math::zero;
+	    return zero(v);
+	}
+	result_type operator()(const Value& v) const
+	{
+	    using math::zero;
+	    return zero(v);
+	}
+    };
+
+    template <typename Value>
+    struct imag<std::complex<Value> >
+    {
+	typedef Value result_type;
+
+	static inline result_type apply(const std::complex<Value>& v)
+	{
+	    return std::imag(v);
+	}
+	result_type operator()(const std::complex<Value>& v) const
+	{
+	    return std::imag(v);
+	}
+    };
+
+}
+
+/// Imaginary part of scalars (including non-complex)
+template <typename Value>
+typename mtl::traits::enable_if_scalar<Value, typename sfunctor::imag<Value>::result_type>::type
+inline imag(const Value& v)
+{
+    return sfunctor::imag<Value>::apply(v);
+}
+
+namespace vector {
+
+    /// Imaginary part of an vector
+    template <typename Vector>
+    typename mtl::traits::enable_if_vector<Vector, imag_view<Vector> >::type
+    inline imag(const Vector& v)
+    {
+	return imag_view<Vector>(v);
+    }
+}
+
+namespace matrix {
+
+    /// Imaginary part of an vector
+    template <typename Matrix>
+    typename mtl::traits::enable_if_matrix<Matrix, imag_view<Matrix> >::type
+    inline imag(const Matrix& v)
+    {
+	return imag_view<Matrix>(v);
+    }
+}
+
+
+} // namespace mtl
+
+#endif // MTL_IMAG_INCLUDE