From 35f7829af10c61e33dd2e2a7a015058e11a11ea0 Mon Sep 17 00:00:00 2001 From: Stanislaw Halik Date: Sat, 25 Mar 2017 14:17:07 +0100 Subject: update --- eigen/Eigen/src/Householder/BlockHouseholder.h | 77 ++++++++++++++++------- eigen/Eigen/src/Householder/CMakeLists.txt | 6 -- eigen/Eigen/src/Householder/Householder.h | 4 +- eigen/Eigen/src/Householder/HouseholderSequence.h | 56 ++++++++++++----- 4 files changed, 100 insertions(+), 43 deletions(-) delete mode 100644 eigen/Eigen/src/Householder/CMakeLists.txt (limited to 'eigen/Eigen/src/Householder') diff --git a/eigen/Eigen/src/Householder/BlockHouseholder.h b/eigen/Eigen/src/Householder/BlockHouseholder.h index 60dbea5..01a7ed1 100644 --- a/eigen/Eigen/src/Householder/BlockHouseholder.h +++ b/eigen/Eigen/src/Householder/BlockHouseholder.h @@ -16,48 +16,83 @@ namespace Eigen { namespace internal { + +/** \internal */ +// template +// void make_block_householder_triangular_factor(TriangularFactorType& triFactor, const VectorsType& vectors, const CoeffsType& hCoeffs) +// { +// typedef typename VectorsType::Scalar Scalar; +// const Index nbVecs = vectors.cols(); +// eigen_assert(triFactor.rows() == nbVecs && triFactor.cols() == nbVecs && vectors.rows()>=nbVecs); +// +// for(Index i = 0; i < nbVecs; i++) +// { +// Index rs = vectors.rows() - i; +// // Warning, note that hCoeffs may alias with vectors. +// // It is then necessary to copy it before modifying vectors(i,i). +// typename CoeffsType::Scalar h = hCoeffs(i); +// // This hack permits to pass trough nested Block<> and Transpose<> expressions. +// Scalar *Vii_ptr = const_cast(vectors.data() + vectors.outerStride()*i + vectors.innerStride()*i); +// Scalar Vii = *Vii_ptr; +// *Vii_ptr = Scalar(1); +// triFactor.col(i).head(i).noalias() = -h * vectors.block(i, 0, rs, i).adjoint() +// * vectors.col(i).tail(rs); +// *Vii_ptr = Vii; +// // FIXME add .noalias() once the triangular product can work inplace +// triFactor.col(i).head(i) = triFactor.block(0,0,i,i).template triangularView() +// * triFactor.col(i).head(i); +// triFactor(i,i) = hCoeffs(i); +// } +// } /** \internal */ +// This variant avoid modifications in vectors template void make_block_householder_triangular_factor(TriangularFactorType& triFactor, const VectorsType& vectors, const CoeffsType& hCoeffs) { - typedef typename TriangularFactorType::Index Index; - typedef typename VectorsType::Scalar Scalar; const Index nbVecs = vectors.cols(); eigen_assert(triFactor.rows() == nbVecs && triFactor.cols() == nbVecs && vectors.rows()>=nbVecs); - for(Index i = 0; i < nbVecs; i++) + for(Index i = nbVecs-1; i >=0 ; --i) { - Index rs = vectors.rows() - i; - Scalar Vii = vectors(i,i); - vectors.const_cast_derived().coeffRef(i,i) = Scalar(1); - triFactor.col(i).head(i).noalias() = -hCoeffs(i) * vectors.block(i, 0, rs, i).adjoint() - * vectors.col(i).tail(rs); - vectors.const_cast_derived().coeffRef(i, i) = Vii; - // FIXME add .noalias() once the triangular product can work inplace - triFactor.col(i).head(i) = triFactor.block(0,0,i,i).template triangularView() - * triFactor.col(i).head(i); + Index rs = vectors.rows() - i - 1; + Index rt = nbVecs-i-1; + + if(rt>0) + { + triFactor.row(i).tail(rt).noalias() = -hCoeffs(i) * vectors.col(i).tail(rs).adjoint() + * vectors.bottomRightCorner(rs, rt).template triangularView(); + + // FIXME add .noalias() once the triangular product can work inplace + triFactor.row(i).tail(rt) = triFactor.row(i).tail(rt) * triFactor.bottomRightCorner(rt,rt).template triangularView(); + + } triFactor(i,i) = hCoeffs(i); } } -/** \internal */ +/** \internal + * if forward then perform mat = H0 * H1 * H2 * mat + * otherwise perform mat = H2 * H1 * H0 * mat + */ template -void apply_block_householder_on_the_left(MatrixType& mat, const VectorsType& vectors, const CoeffsType& hCoeffs) +void apply_block_householder_on_the_left(MatrixType& mat, const VectorsType& vectors, const CoeffsType& hCoeffs, bool forward) { - typedef typename MatrixType::Index Index; enum { TFactorSize = MatrixType::ColsAtCompileTime }; Index nbVecs = vectors.cols(); - Matrix T(nbVecs,nbVecs); - make_block_householder_triangular_factor(T, vectors, hCoeffs); - - const TriangularView& V(vectors); + Matrix T(nbVecs,nbVecs); + + if(forward) make_block_householder_triangular_factor(T, vectors, hCoeffs); + else make_block_householder_triangular_factor(T, vectors, hCoeffs.conjugate()); + const TriangularView V(vectors); // A -= V T V^* A - Matrix tmp = V.adjoint() * mat; // FIXME add .noalias() once the triangular product can work inplace - tmp = T.template triangularView().adjoint() * tmp; + if(forward) tmp = T.template triangularView() * tmp; + else tmp = T.template triangularView().adjoint() * tmp; mat.noalias() -= V * tmp; } diff --git a/eigen/Eigen/src/Householder/CMakeLists.txt b/eigen/Eigen/src/Householder/CMakeLists.txt deleted file mode 100644 index ce4937d..0000000 --- a/eigen/Eigen/src/Householder/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_Householder_SRCS "*.h") - -INSTALL(FILES - ${Eigen_Householder_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Householder COMPONENT Devel - ) diff --git a/eigen/Eigen/src/Householder/Householder.h b/eigen/Eigen/src/Householder/Householder.h index 4c1f499..80de2c3 100644 --- a/eigen/Eigen/src/Householder/Householder.h +++ b/eigen/Eigen/src/Householder/Householder.h @@ -119,7 +119,7 @@ void MatrixBase::applyHouseholderOnTheLeft( { *this *= Scalar(1)-tau; } - else + else if(tau!=Scalar(0)) { Map::type> tmp(workspace,cols()); Block bottom(derived(), 1, 0, rows()-1, cols()); @@ -156,7 +156,7 @@ void MatrixBase::applyHouseholderOnTheRight( { *this *= Scalar(1)-tau; } - else + else if(tau!=Scalar(0)) { Map::type> tmp(workspace,rows()); Block right(derived(), 0, 1, rows(), cols()-1); diff --git a/eigen/Eigen/src/Householder/HouseholderSequence.h b/eigen/Eigen/src/Householder/HouseholderSequence.h index aea2439..3ce0a69 100644 --- a/eigen/Eigen/src/Householder/HouseholderSequence.h +++ b/eigen/Eigen/src/Householder/HouseholderSequence.h @@ -60,7 +60,7 @@ template struct traits > { typedef typename VectorsType::Scalar Scalar; - typedef typename VectorsType::Index Index; + typedef typename VectorsType::StorageIndex StorageIndex; typedef typename VectorsType::StorageKind StorageKind; enum { RowsAtCompileTime = Side==OnTheLeft ? traits::RowsAtCompileTime @@ -73,12 +73,20 @@ struct traits > }; }; +struct HouseholderSequenceShape {}; + +template +struct evaluator_traits > + : public evaluator_traits_base > +{ + typedef HouseholderSequenceShape Shape; +}; + template struct hseq_side_dependent_impl { typedef Block EssentialVectorType; typedef HouseholderSequence HouseholderSequenceType; - typedef typename VectorsType::Index Index; static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k) { Index start = k+1+h.m_shift; @@ -91,7 +99,6 @@ struct hseq_side_dependent_impl { typedef Transpose > EssentialVectorType; typedef HouseholderSequence HouseholderSequenceType; - typedef typename VectorsType::Index Index; static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k) { Index start = k+1+h.m_shift; @@ -101,7 +108,7 @@ struct hseq_side_dependent_impl template struct matrix_type_times_scalar_type { - typedef typename scalar_product_traits::ReturnType + typedef typename ScalarBinaryOpTraits::ReturnType ResultScalar; typedef Matrix Type; @@ -122,7 +129,6 @@ template class HouseholderS MaxColsAtCompileTime = internal::traits::MaxColsAtCompileTime }; typedef typename internal::traits::Scalar Scalar; - typedef typename VectorsType::Index Index; typedef HouseholderSequence< typename internal::conditional::IsComplex, @@ -237,9 +243,7 @@ template class HouseholderS { workspace.resize(rows()); Index vecs = m_length; - const typename Dest::Scalar *dst_data = internal::extract_data(dst); - if( internal::is_same::type,Dest>::value - && dst_data!=0 && dst_data == internal::extract_data(m_vectors)) + if(internal::is_same_dense(dst,m_vectors)) { // in-place dst.diagonal().setOnes(); @@ -300,7 +304,7 @@ template class HouseholderS /** \internal */ template inline void applyThisOnTheLeft(Dest& dst) const { - Matrix workspace(dst.cols()); + Matrix workspace; applyThisOnTheLeft(dst, workspace); } @@ -308,12 +312,36 @@ template class HouseholderS template inline void applyThisOnTheLeft(Dest& dst, Workspace& workspace) const { - workspace.resize(dst.cols()); - for(Index k = 0; k < m_length; ++k) + const Index BlockSize = 48; + // if the entries are large enough, then apply the reflectors by block + if(m_length>=BlockSize && dst.cols()>1) { - Index actual_k = m_trans ? k : m_length-k-1; - dst.bottomRows(rows()-m_shift-actual_k) - .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data()); + for(Index i = 0; i < m_length; i+=BlockSize) + { + Index end = m_trans ? (std::min)(m_length,i+BlockSize) : m_length-i; + Index k = m_trans ? i : (std::max)(Index(0),end-BlockSize); + Index bs = end-k; + Index start = k + m_shift; + + typedef Block::type,Dynamic,Dynamic> SubVectorsType; + SubVectorsType sub_vecs1(m_vectors.const_cast_derived(), Side==OnTheRight ? k : start, + Side==OnTheRight ? start : k, + Side==OnTheRight ? bs : m_vectors.rows()-start, + Side==OnTheRight ? m_vectors.cols()-start : bs); + typename internal::conditional, SubVectorsType&>::type sub_vecs(sub_vecs1); + Block sub_dst(dst,dst.rows()-rows()+m_shift+k,0, rows()-m_shift-k,dst.cols()); + apply_block_householder_on_the_left(sub_dst, sub_vecs, m_coeffs.segment(k, bs), !m_trans); + } + } + else + { + workspace.resize(dst.cols()); + for(Index k = 0; k < m_length; ++k) + { + Index actual_k = m_trans ? k : m_length-k-1; + dst.bottomRows(rows()-m_shift-actual_k) + .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data()); + } } } -- cgit v1.2.3