1f49d36a214375ac9a735ab9af9e1b3a4b9d2a3cdSergey Vasilinets// This file is part of Eigen, a lightweight C++ template library 2f49d36a214375ac9a735ab9af9e1b3a4b9d2a3cdSergey Vasilinets// for linear algebra. 3f49d36a214375ac9a735ab9af9e1b3a4b9d2a3cdSergey Vasilinets// 4f49d36a214375ac9a735ab9af9e1b3a4b9d2a3cdSergey Vasilinets// Copyright (C) 2010 Vincent Lejeune 5f49d36a214375ac9a735ab9af9e1b3a4b9d2a3cdSergey Vasilinets// Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr> 6f49d36a214375ac9a735ab9af9e1b3a4b9d2a3cdSergey Vasilinets// 7f49d36a214375ac9a735ab9af9e1b3a4b9d2a3cdSergey Vasilinets// This Source Code Form is subject to the terms of the Mozilla 8f49d36a214375ac9a735ab9af9e1b3a4b9d2a3cdSergey Vasilinets// Public License v. 2.0. If a copy of the MPL was not distributed 9f49d36a214375ac9a735ab9af9e1b3a4b9d2a3cdSergey Vasilinets// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10f49d36a214375ac9a735ab9af9e1b3a4b9d2a3cdSergey Vasilinets 11f49d36a214375ac9a735ab9af9e1b3a4b9d2a3cdSergey Vasilinets#ifndef EIGEN_BLOCK_HOUSEHOLDER_H 12f49d36a214375ac9a735ab9af9e1b3a4b9d2a3cdSergey Vasilinets#define EIGEN_BLOCK_HOUSEHOLDER_H 13f49d36a214375ac9a735ab9af9e1b3a4b9d2a3cdSergey Vasilinets 14f49d36a214375ac9a735ab9af9e1b3a4b9d2a3cdSergey Vasilinets// This file contains some helper function to deal with block householder reflectors 15f49d36a214375ac9a735ab9af9e1b3a4b9d2a3cdSergey Vasilinets 16f49d36a214375ac9a735ab9af9e1b3a4b9d2a3cdSergey Vasilinetsnamespace Eigen { 171503d52153986fdcfe7e744795010708b7410892Ian Lake 18f49d36a214375ac9a735ab9af9e1b3a4b9d2a3cdSergey Vasilinetsnamespace internal { 19f49d36a214375ac9a735ab9af9e1b3a4b9d2a3cdSergey Vasilinets 20f49d36a214375ac9a735ab9af9e1b3a4b9d2a3cdSergey Vasilinets/** \internal */ 21f49d36a214375ac9a735ab9af9e1b3a4b9d2a3cdSergey Vasilinets// 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 VectorsType::Scalar Scalar; 25// const Index nbVecs = vectors.cols(); 26// eigen_assert(triFactor.rows() == nbVecs && triFactor.cols() == nbVecs && vectors.rows()>=nbVecs); 27// 28// for(Index i = 0; i < nbVecs; i++) 29// { 30// Index rs = vectors.rows() - i; 31// // Warning, note that hCoeffs may alias with vectors. 32// // It is then necessary to copy it before modifying vectors(i,i). 33// typename CoeffsType::Scalar h = hCoeffs(i); 34// // This hack permits to pass trough nested Block<> and Transpose<> expressions. 35// Scalar *Vii_ptr = const_cast<Scalar*>(vectors.data() + vectors.outerStride()*i + vectors.innerStride()*i); 36// Scalar Vii = *Vii_ptr; 37// *Vii_ptr = Scalar(1); 38// triFactor.col(i).head(i).noalias() = -h * vectors.block(i, 0, rs, i).adjoint() 39// * vectors.col(i).tail(rs); 40// *Vii_ptr = Vii; 41// // FIXME add .noalias() once the triangular product can work inplace 42// triFactor.col(i).head(i) = triFactor.block(0,0,i,i).template triangularView<Upper>() 43// * triFactor.col(i).head(i); 44// triFactor(i,i) = hCoeffs(i); 45// } 46// } 47 48/** \internal */ 49// This variant avoid modifications in vectors 50template<typename TriangularFactorType,typename VectorsType,typename CoeffsType> 51void make_block_householder_triangular_factor(TriangularFactorType& triFactor, const VectorsType& vectors, const CoeffsType& hCoeffs) 52{ 53 const Index nbVecs = vectors.cols(); 54 eigen_assert(triFactor.rows() == nbVecs && triFactor.cols() == nbVecs && vectors.rows()>=nbVecs); 55 56 for(Index i = nbVecs-1; i >=0 ; --i) 57 { 58 Index rs = vectors.rows() - i - 1; 59 Index rt = nbVecs-i-1; 60 61 if(rt>0) 62 { 63 triFactor.row(i).tail(rt).noalias() = -hCoeffs(i) * vectors.col(i).tail(rs).adjoint() 64 * vectors.bottomRightCorner(rs, rt).template triangularView<UnitLower>(); 65 66 // FIXME add .noalias() once the triangular product can work inplace 67 triFactor.row(i).tail(rt) = triFactor.row(i).tail(rt) * triFactor.bottomRightCorner(rt,rt).template triangularView<Upper>(); 68 69 } 70 triFactor(i,i) = hCoeffs(i); 71 } 72} 73 74/** \internal 75 * if forward then perform mat = H0 * H1 * H2 * mat 76 * otherwise perform mat = H2 * H1 * H0 * mat 77 */ 78template<typename MatrixType,typename VectorsType,typename CoeffsType> 79void apply_block_householder_on_the_left(MatrixType& mat, const VectorsType& vectors, const CoeffsType& hCoeffs, bool forward) 80{ 81 enum { TFactorSize = MatrixType::ColsAtCompileTime }; 82 Index nbVecs = vectors.cols(); 83 Matrix<typename MatrixType::Scalar, TFactorSize, TFactorSize, RowMajor> T(nbVecs,nbVecs); 84 85 if(forward) make_block_householder_triangular_factor(T, vectors, hCoeffs); 86 else make_block_householder_triangular_factor(T, vectors, hCoeffs.conjugate()); 87 const TriangularView<const VectorsType, UnitLower> V(vectors); 88 89 // A -= V T V^* A 90 Matrix<typename MatrixType::Scalar,VectorsType::ColsAtCompileTime,MatrixType::ColsAtCompileTime, 91 (VectorsType::MaxColsAtCompileTime==1 && MatrixType::MaxColsAtCompileTime!=1)?RowMajor:ColMajor, 92 VectorsType::MaxColsAtCompileTime,MatrixType::MaxColsAtCompileTime> tmp = V.adjoint() * mat; 93 // FIXME add .noalias() once the triangular product can work inplace 94 if(forward) tmp = T.template triangularView<Upper>() * tmp; 95 else tmp = T.template triangularView<Upper>().adjoint() * tmp; 96 mat.noalias() -= V * tmp; 97} 98 99} // end namespace internal 100 101} // end namespace Eigen 102 103#endif // EIGEN_BLOCK_HOUSEHOLDER_H 104