summaryrefslogtreecommitdiffhomepage
path: root/eigen/Eigen/src/Householder
diff options
context:
space:
mode:
Diffstat (limited to 'eigen/Eigen/src/Householder')
-rw-r--r--eigen/Eigen/src/Householder/BlockHouseholder.h77
-rw-r--r--eigen/Eigen/src/Householder/CMakeLists.txt6
-rw-r--r--eigen/Eigen/src/Householder/Householder.h4
-rw-r--r--eigen/Eigen/src/Householder/HouseholderSequence.h56
4 files changed, 100 insertions, 43 deletions
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<typename TriangularFactorType,typename VectorsType,typename CoeffsType>
+// 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<Scalar*>(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<Upper>()
+// * triFactor.col(i).head(i);
+// triFactor(i,i) = hCoeffs(i);
+// }
+// }
/** \internal */
+// This variant avoid modifications in vectors
template<typename TriangularFactorType,typename VectorsType,typename CoeffsType>
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<Upper>()
- * 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<UnitLower>();
+
+ // 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<Upper>();
+
+ }
triFactor(i,i) = hCoeffs(i);
}
}
-/** \internal */
+/** \internal
+ * if forward then perform mat = H0 * H1 * H2 * mat
+ * otherwise perform mat = H2 * H1 * H0 * mat
+ */
template<typename MatrixType,typename VectorsType,typename CoeffsType>
-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<typename MatrixType::Scalar, TFactorSize, TFactorSize, ColMajor> T(nbVecs,nbVecs);
- make_block_householder_triangular_factor(T, vectors, hCoeffs);
-
- const TriangularView<const VectorsType, UnitLower>& V(vectors);
+ Matrix<typename MatrixType::Scalar, TFactorSize, TFactorSize, RowMajor> 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<const VectorsType, UnitLower> V(vectors);
// A -= V T V^* A
- Matrix<typename MatrixType::Scalar,VectorsType::ColsAtCompileTime,MatrixType::ColsAtCompileTime,0,
+ Matrix<typename MatrixType::Scalar,VectorsType::ColsAtCompileTime,MatrixType::ColsAtCompileTime,
+ (VectorsType::MaxColsAtCompileTime==1 && MatrixType::MaxColsAtCompileTime!=1)?RowMajor:ColMajor,
VectorsType::MaxColsAtCompileTime,MatrixType::MaxColsAtCompileTime> tmp = V.adjoint() * mat;
// FIXME add .noalias() once the triangular product can work inplace
- tmp = T.template triangularView<Upper>().adjoint() * tmp;
+ if(forward) tmp = T.template triangularView<Upper>() * tmp;
+ else tmp = T.template triangularView<Upper>().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<Derived>::applyHouseholderOnTheLeft(
{
*this *= Scalar(1)-tau;
}
- else
+ else if(tau!=Scalar(0))
{
Map<typename internal::plain_row_type<PlainObject>::type> tmp(workspace,cols());
Block<Derived, EssentialPart::SizeAtCompileTime, Derived::ColsAtCompileTime> bottom(derived(), 1, 0, rows()-1, cols());
@@ -156,7 +156,7 @@ void MatrixBase<Derived>::applyHouseholderOnTheRight(
{
*this *= Scalar(1)-tau;
}
- else
+ else if(tau!=Scalar(0))
{
Map<typename internal::plain_col_type<PlainObject>::type> tmp(workspace,rows());
Block<Derived, Derived::RowsAtCompileTime, EssentialPart::SizeAtCompileTime> 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<typename VectorsType, typename CoeffsType, int Side>
struct traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
{
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<VectorsType>::RowsAtCompileTime
@@ -73,12 +73,20 @@ struct traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
};
};
+struct HouseholderSequenceShape {};
+
+template<typename VectorsType, typename CoeffsType, int Side>
+struct evaluator_traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
+ : public evaluator_traits_base<HouseholderSequence<VectorsType,CoeffsType,Side> >
+{
+ typedef HouseholderSequenceShape Shape;
+};
+
template<typename VectorsType, typename CoeffsType, int Side>
struct hseq_side_dependent_impl
{
typedef Block<const VectorsType, Dynamic, 1> EssentialVectorType;
typedef HouseholderSequence<VectorsType, CoeffsType, OnTheLeft> 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<VectorsType, CoeffsType, OnTheRight>
{
typedef Transpose<Block<const VectorsType, 1, Dynamic> > EssentialVectorType;
typedef HouseholderSequence<VectorsType, CoeffsType, OnTheRight> 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<VectorsType, CoeffsType, OnTheRight>
template<typename OtherScalarType, typename MatrixType> struct matrix_type_times_scalar_type
{
- typedef typename scalar_product_traits<OtherScalarType, typename MatrixType::Scalar>::ReturnType
+ typedef typename ScalarBinaryOpTraits<OtherScalarType, typename MatrixType::Scalar>::ReturnType
ResultScalar;
typedef Matrix<ResultScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime,
0, MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime> Type;
@@ -122,7 +129,6 @@ template<typename VectorsType, typename CoeffsType, int Side> class HouseholderS
MaxColsAtCompileTime = internal::traits<HouseholderSequence>::MaxColsAtCompileTime
};
typedef typename internal::traits<HouseholderSequence>::Scalar Scalar;
- typedef typename VectorsType::Index Index;
typedef HouseholderSequence<
typename internal::conditional<NumTraits<Scalar>::IsComplex,
@@ -237,9 +243,7 @@ template<typename VectorsType, typename CoeffsType, int Side> class HouseholderS
{
workspace.resize(rows());
Index vecs = m_length;
- const typename Dest::Scalar *dst_data = internal::extract_data(dst);
- if( internal::is_same<typename internal::remove_all<VectorsType>::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<typename VectorsType, typename CoeffsType, int Side> class HouseholderS
/** \internal */
template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
{
- Matrix<Scalar,1,Dest::ColsAtCompileTime,RowMajor,1,Dest::MaxColsAtCompileTime> workspace(dst.cols());
+ Matrix<Scalar,1,Dest::ColsAtCompileTime,RowMajor,1,Dest::MaxColsAtCompileTime> workspace;
applyThisOnTheLeft(dst, workspace);
}
@@ -308,12 +312,36 @@ template<typename VectorsType, typename CoeffsType, int Side> class HouseholderS
template<typename Dest, typename Workspace>
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<typename internal::remove_all<VectorsType>::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<Side==OnTheRight, Transpose<SubVectorsType>, SubVectorsType&>::type sub_vecs(sub_vecs1);
+ Block<Dest,Dynamic,Dynamic> 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());
+ }
}
}
generated by cgit v1.2.3 (git 2.39.1) at 2025-02-23 17:30:05 +0000