![]() |
Eigen
3.3.3
|
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 VectorsType::Scalar Scalar; 00025 // const Index nbVecs = vectors.cols(); 00026 // eigen_assert(triFactor.rows() == nbVecs && triFactor.cols() == nbVecs && vectors.rows()>=nbVecs); 00027 // 00028 // for(Index i = 0; i < nbVecs; i++) 00029 // { 00030 // Index rs = vectors.rows() - i; 00031 // // Warning, note that hCoeffs may alias with vectors. 00032 // // It is then necessary to copy it before modifying vectors(i,i). 00033 // typename CoeffsType::Scalar h = hCoeffs(i); 00034 // // This hack permits to pass trough nested Block<> and Transpose<> expressions. 00035 // Scalar *Vii_ptr = const_cast<Scalar*>(vectors.data() + vectors.outerStride()*i + vectors.innerStride()*i); 00036 // Scalar Vii = *Vii_ptr; 00037 // *Vii_ptr = Scalar(1); 00038 // triFactor.col(i).head(i).noalias() = -h * vectors.block(i, 0, rs, i).adjoint() 00039 // * vectors.col(i).tail(rs); 00040 // *Vii_ptr = Vii; 00041 // // FIXME add .noalias() once the triangular product can work inplace 00042 // triFactor.col(i).head(i) = triFactor.block(0,0,i,i).template triangularView<Upper>() 00043 // * triFactor.col(i).head(i); 00044 // triFactor(i,i) = hCoeffs(i); 00045 // } 00046 // } 00047 00049 // This variant avoid modifications in vectors 00050 template<typename TriangularFactorType,typename VectorsType,typename CoeffsType> 00051 void make_block_householder_triangular_factor(TriangularFactorType& triFactor, const VectorsType& vectors, const CoeffsType& hCoeffs) 00052 { 00053 const Index nbVecs = vectors.cols(); 00054 eigen_assert(triFactor.rows() == nbVecs && triFactor.cols() == nbVecs && vectors.rows()>=nbVecs); 00055 00056 for(Index i = nbVecs-1; i >=0 ; --i) 00057 { 00058 Index rs = vectors.rows() - i - 1; 00059 Index rt = nbVecs-i-1; 00060 00061 if(rt>0) 00062 { 00063 triFactor.row(i).tail(rt).noalias() = -hCoeffs(i) * vectors.col(i).tail(rs).adjoint() 00064 * vectors.bottomRightCorner(rs, rt).template triangularView<UnitLower>(); 00065 00066 // FIXME add .noalias() once the triangular product can work inplace 00067 triFactor.row(i).tail(rt) = triFactor.row(i).tail(rt) * triFactor.bottomRightCorner(rt,rt).template triangularView<Upper>(); 00068 00069 } 00070 triFactor(i,i) = hCoeffs(i); 00071 } 00072 } 00073 00078 template<typename MatrixType,typename VectorsType,typename CoeffsType> 00079 void apply_block_householder_on_the_left(MatrixType& mat, const VectorsType& vectors, const CoeffsType& hCoeffs, bool forward) 00080 { 00081 enum { TFactorSize = MatrixType::ColsAtCompileTime }; 00082 Index nbVecs = vectors.cols(); 00083 Matrix<typename MatrixType::Scalar, TFactorSize, TFactorSize, RowMajor> T(nbVecs,nbVecs); 00084 00085 if(forward) make_block_householder_triangular_factor(T, vectors, hCoeffs); 00086 else make_block_householder_triangular_factor(T, vectors, hCoeffs.conjugate()); 00087 const TriangularView<const VectorsType, UnitLower> V(vectors); 00088 00089 // A -= V T V^* A 00090 Matrix<typename MatrixType::Scalar,VectorsType::ColsAtCompileTime,MatrixType::ColsAtCompileTime, 00091 (VectorsType::MaxColsAtCompileTime==1 && MatrixType::MaxColsAtCompileTime!=1)?RowMajor:ColMajor, 00092 VectorsType::MaxColsAtCompileTime,MatrixType::MaxColsAtCompileTime> tmp = V.adjoint() * mat; 00093 // FIXME add .noalias() once the triangular product can work inplace 00094 if(forward) tmp = T.template triangularView<Upper>() * tmp; 00095 else tmp = T.template triangularView<Upper>().adjoint() * tmp; 00096 mat.noalias() -= V * tmp; 00097 } 00098 00099 } // end namespace internal 00100 00101 } // end namespace Eigen 00102 00103 #endif // EIGEN_BLOCK_HOUSEHOLDER_H