Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #ifndef EIGEN_BLOCK_HOUSEHOLDER_H
00027 #define EIGEN_BLOCK_HOUSEHOLDER_H
00028
00029
00030
00031 namespace Eigen {
00032
00033 namespace internal {
00034
00036 template<typename TriangularFactorType,typename VectorsType,typename CoeffsType>
00037 void make_block_householder_triangular_factor(TriangularFactorType& triFactor, const VectorsType& vectors, const CoeffsType& hCoeffs)
00038 {
00039 typedef typename TriangularFactorType::Index Index;
00040 typedef typename VectorsType::Scalar Scalar;
00041 const Index nbVecs = vectors.cols();
00042 eigen_assert(triFactor.rows() == nbVecs && triFactor.cols() == nbVecs && vectors.rows()>=nbVecs);
00043
00044 for(Index i = 0; i < nbVecs; i++)
00045 {
00046 Index rs = vectors.rows() - i;
00047 Scalar Vii = vectors(i,i);
00048 vectors.const_cast_derived().coeffRef(i,i) = Scalar(1);
00049 triFactor.col(i).head(i).noalias() = -hCoeffs(i) * vectors.block(i, 0, rs, i).adjoint()
00050 * vectors.col(i).tail(rs);
00051 vectors.const_cast_derived().coeffRef(i, i) = Vii;
00052
00053 triFactor.col(i).head(i) = triFactor.block(0,0,i,i).template triangularView<Upper>()
00054 * triFactor.col(i).head(i);
00055 triFactor(i,i) = hCoeffs(i);
00056 }
00057 }
00058
00060 template<typename MatrixType,typename VectorsType,typename CoeffsType>
00061 void apply_block_householder_on_the_left(MatrixType& mat, const VectorsType& vectors, const CoeffsType& hCoeffs)
00062 {
00063 typedef typename MatrixType::Index Index;
00064 enum { TFactorSize = MatrixType::ColsAtCompileTime };
00065 Index nbVecs = vectors.cols();
00066 Matrix<typename MatrixType::Scalar, TFactorSize, TFactorSize> T(nbVecs,nbVecs);
00067 make_block_householder_triangular_factor(T, vectors, hCoeffs);
00068
00069 const TriangularView<const VectorsType, UnitLower>& V(vectors);
00070
00071
00072 Matrix<typename MatrixType::Scalar,VectorsType::ColsAtCompileTime,MatrixType::ColsAtCompileTime,0,
00073 VectorsType::MaxColsAtCompileTime,MatrixType::MaxColsAtCompileTime> tmp = V.adjoint() * mat;
00074
00075 tmp = T.template triangularView<Upper>().adjoint() * tmp;
00076 mat.noalias() -= V * tmp;
00077 }
00078
00079 }
00080
00081 }
00082
00083 #endif // EIGEN_BLOCK_HOUSEHOLDER_H