1 | // This file is part of Eigen, a lightweight C++ template library
|
---|
2 | // for linear algebra.
|
---|
3 | //
|
---|
4 | // Copyright (C) 2010 Vincent Lejeune
|
---|
5 | // Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr>
|
---|
6 | //
|
---|
7 | // This Source Code Form is subject to the terms of the Mozilla
|
---|
8 | // Public License v. 2.0. If a copy of the MPL was not distributed
|
---|
9 | // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
---|
10 |
|
---|
11 | #ifndef EIGEN_BLOCK_HOUSEHOLDER_H
|
---|
12 | #define EIGEN_BLOCK_HOUSEHOLDER_H
|
---|
13 |
|
---|
14 | // This file contains some helper function to deal with block householder reflectors
|
---|
15 |
|
---|
16 | namespace Eigen {
|
---|
17 |
|
---|
18 | namespace internal {
|
---|
19 |
|
---|
20 | /** \internal */
|
---|
21 | template<typename TriangularFactorType,typename VectorsType,typename CoeffsType>
|
---|
22 | void make_block_householder_triangular_factor(TriangularFactorType& triFactor, const VectorsType& vectors, const CoeffsType& hCoeffs)
|
---|
23 | {
|
---|
24 | typedef typename TriangularFactorType::Index Index;
|
---|
25 | typedef typename VectorsType::Scalar Scalar;
|
---|
26 | const Index nbVecs = vectors.cols();
|
---|
27 | eigen_assert(triFactor.rows() == nbVecs && triFactor.cols() == nbVecs && vectors.rows()>=nbVecs);
|
---|
28 |
|
---|
29 | for(Index i = 0; i < nbVecs; i++)
|
---|
30 | {
|
---|
31 | Index rs = vectors.rows() - i;
|
---|
32 | Scalar Vii = vectors(i,i);
|
---|
33 | vectors.const_cast_derived().coeffRef(i,i) = Scalar(1);
|
---|
34 | triFactor.col(i).head(i).noalias() = -hCoeffs(i) * vectors.block(i, 0, rs, i).adjoint()
|
---|
35 | * vectors.col(i).tail(rs);
|
---|
36 | vectors.const_cast_derived().coeffRef(i, i) = Vii;
|
---|
37 | // FIXME add .noalias() once the triangular product can work inplace
|
---|
38 | triFactor.col(i).head(i) = triFactor.block(0,0,i,i).template triangularView<Upper>()
|
---|
39 | * triFactor.col(i).head(i);
|
---|
40 | triFactor(i,i) = hCoeffs(i);
|
---|
41 | }
|
---|
42 | }
|
---|
43 |
|
---|
44 | /** \internal */
|
---|
45 | template<typename MatrixType,typename VectorsType,typename CoeffsType>
|
---|
46 | void apply_block_householder_on_the_left(MatrixType& mat, const VectorsType& vectors, const CoeffsType& hCoeffs)
|
---|
47 | {
|
---|
48 | typedef typename MatrixType::Index Index;
|
---|
49 | enum { TFactorSize = MatrixType::ColsAtCompileTime };
|
---|
50 | Index nbVecs = vectors.cols();
|
---|
51 | Matrix<typename MatrixType::Scalar, TFactorSize, TFactorSize, ColMajor> T(nbVecs,nbVecs);
|
---|
52 | make_block_householder_triangular_factor(T, vectors, hCoeffs);
|
---|
53 |
|
---|
54 | const TriangularView<const VectorsType, UnitLower>& V(vectors);
|
---|
55 |
|
---|
56 | // A -= V T V^* A
|
---|
57 | Matrix<typename MatrixType::Scalar,VectorsType::ColsAtCompileTime,MatrixType::ColsAtCompileTime,0,
|
---|
58 | VectorsType::MaxColsAtCompileTime,MatrixType::MaxColsAtCompileTime> tmp = V.adjoint() * mat;
|
---|
59 | // FIXME add .noalias() once the triangular product can work inplace
|
---|
60 | tmp = T.template triangularView<Upper>().adjoint() * tmp;
|
---|
61 | mat.noalias() -= V * tmp;
|
---|
62 | }
|
---|
63 |
|
---|
64 | } // end namespace internal
|
---|
65 |
|
---|
66 | } // end namespace Eigen
|
---|
67 |
|
---|
68 | #endif // EIGEN_BLOCK_HOUSEHOLDER_H
|
---|