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