$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2010 Vincent Lejeune 00005 // Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr> 00006 // 00007 // This Source Code Form is subject to the terms of the Mozilla 00008 // Public License v. 2.0. If a copy of the MPL was not distributed 00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00010 00011 #ifndef EIGEN_BLOCK_HOUSEHOLDER_H 00012 #define EIGEN_BLOCK_HOUSEHOLDER_H 00013 00014 // This file contains some helper function to deal with block householder reflectors 00015 00016 namespace Eigen { 00017 00018 namespace internal { 00019 00021 template<typename TriangularFactorType,typename VectorsType,typename CoeffsType> 00022 void make_block_householder_triangular_factor(TriangularFactorType& triFactor, const VectorsType& vectors, const CoeffsType& hCoeffs) 00023 { 00024 typedef typename TriangularFactorType::Index Index; 00025 typedef typename VectorsType::Scalar Scalar; 00026 const Index nbVecs = vectors.cols(); 00027 eigen_assert(triFactor.rows() == nbVecs && triFactor.cols() == nbVecs && vectors.rows()>=nbVecs); 00028 00029 for(Index i = 0; i < nbVecs; i++) 00030 { 00031 Index rs = vectors.rows() - i; 00032 Scalar Vii = vectors(i,i); 00033 vectors.const_cast_derived().coeffRef(i,i) = Scalar(1); 00034 triFactor.col(i).head(i).noalias() = -hCoeffs(i) * vectors.block(i, 0, rs, i).adjoint() 00035 * vectors.col(i).tail(rs); 00036 vectors.const_cast_derived().coeffRef(i, i) = Vii; 00037 // FIXME add .noalias() once the triangular product can work inplace 00038 triFactor.col(i).head(i) = triFactor.block(0,0,i,i).template triangularView<Upper>() 00039 * triFactor.col(i).head(i); 00040 triFactor(i,i) = hCoeffs(i); 00041 } 00042 } 00043 00045 template<typename MatrixType,typename VectorsType,typename CoeffsType> 00046 void apply_block_householder_on_the_left(MatrixType& mat, const VectorsType& vectors, const CoeffsType& hCoeffs) 00047 { 00048 typedef typename MatrixType::Index Index; 00049 enum { TFactorSize = MatrixType::ColsAtCompileTime }; 00050 Index nbVecs = vectors.cols(); 00051 Matrix<typename MatrixType::Scalar, TFactorSize, TFactorSize, ColMajor> T(nbVecs,nbVecs); 00052 make_block_householder_triangular_factor(T, vectors, hCoeffs); 00053 00054 const TriangularView<const VectorsType, UnitLower>& V(vectors); 00055 00056 // A -= V T V^* A 00057 Matrix<typename MatrixType::Scalar,VectorsType::ColsAtCompileTime,MatrixType::ColsAtCompileTime,0, 00058 VectorsType::MaxColsAtCompileTime,MatrixType::MaxColsAtCompileTime> tmp = V.adjoint() * mat; 00059 // FIXME add .noalias() once the triangular product can work inplace 00060 tmp = T.template triangularView<Upper>().adjoint() * tmp; 00061 mat.noalias() -= V * tmp; 00062 } 00063 00064 } // end namespace internal 00065 00066 } // end namespace Eigen 00067 00068 #endif // EIGEN_BLOCK_HOUSEHOLDER_H