| 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> 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 |
|---|
