householder.hpp 2.67 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
// 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 {


29
30
/// Computes Householder vector v and scalar b for vector \p y
/** such that identity_matrix(size(y))-b*v*v' projects the vector y
31
32
33
34
35
36
37
38
39
40
41
42
43
    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;
44
    irange            tail(1, imax);
45
46
47
48
49
50
51
52
53
54
55
    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
    }
56

57
58
59
    return std::make_pair(v,b);
}

60
/// Computes Householder vector for vector \p y
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
/** 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 {
77
		s= 1;
78
79
80
81
82
83
84
85
86
87
88
89
90
	}
	u[0]= u[0] + s * nu;
	u/= sqrt( dot(u, u) );
    }

    return u;
}


}} // namespace mtl::matrix

#endif // MTL_MATRIX_HOUSEHOLDER_INCLUDE