mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-18 14:34:17 +08:00
Code Issues Packages Projects Releases Wiki Activity

fix HouseholderSequence API, bug #50:

Browse Source 
* remove ctors taking more than 2 ints
 * rename actualVectors to length
 * add length/shift/trans accessors/mutators
This commit is contained in:
Benoit Jacob 2010-12-30 04:18:40 -05:00
parent e112ad8124
commit dbd9c5fd50
6 changed files with 78 additions and 49 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
Eigen/src/Eigenvalues/HessenbergDecomposition.h
View File

@ -245,7 +245,9 @@ template<typename _MatrixType> class HessenbergDecomposition
HouseholderSequenceType matrixQ() const
{
eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate(), false, m_matrix.rows() - 1, 1);
return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
.setLength(m_matrix.rows() - 1)
.setShift(1);
}
/** \brief Constructs the Hessenberg matrix H in the decomposition

8
Eigen/src/Eigenvalues/Tridiagonalization.h
View File

@ -251,7 +251,9 @@ template<typename _MatrixType> class Tridiagonalization
HouseholderSequenceType matrixQ() const
{
eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate(), false, m_matrix.rows() - 1, 1);
return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
.setLength(m_matrix.rows() - 1)
.setShift(1);
}
/** \brief Returns an expression of the tridiagonal matrix T in the decomposition
@ -459,7 +461,9 @@ struct tridiagonalization_inplace_selector
diag = mat.diagonal().real();
subdiag = mat.template diagonal<-1>().real();
if(extractQ)
mat = HouseholderSequenceType(mat, hCoeffs.conjugate(), false, mat.rows() - 1, 1);
mat = HouseholderSequenceType(mat, hCoeffs.conjugate())
.setLength(mat.rows() - 1)
.setShift(1);
}
};

87
Eigen/src/Householder/HouseholderSequence.h
View File

@ -129,14 +129,18 @@ template<typename VectorsType, typename CoeffsType, int Side> class HouseholderS
Side
> ConjugateReturnType;
HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans = false)
: m_vectors(v), m_coeffs(h), m_trans(trans), m_actualVectors(v.diagonalSize()),
HouseholderSequence(const VectorsType& v, const CoeffsType& h)
: m_vectors(v), m_coeffs(h), m_trans(false), m_length(v.diagonalSize()),
m_shift(0)
{
}
HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans, Index actualVectors, Index shift)
: m_vectors(v), m_coeffs(h), m_trans(trans), m_actualVectors(actualVectors), m_shift(shift)
HouseholderSequence(const HouseholderSequence& other)
: m_vectors(other.m_vectors),
m_coeffs(other.m_coeffs),
m_trans(other.m_trans),
m_length(other.m_length),
m_shift(other.m_shift)
{
}
@ -145,25 +149,34 @@ template<typename VectorsType, typename CoeffsType, int Side> class HouseholderS
const EssentialVectorType essentialVector(Index k) const
{
eigen_assert(k >= 0 && k < m_actualVectors);
eigen_assert(k >= 0 && k < m_length);
return internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*this, k);
}
HouseholderSequence transpose() const
{ return HouseholderSequence(m_vectors, m_coeffs, !m_trans, m_actualVectors, m_shift); }
{
return HouseholderSequence(*this).setTrans(!m_trans);
}
ConjugateReturnType conjugate() const
{ return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans, m_actualVectors, m_shift); }
{
return ConjugateReturnType(m_vectors, m_coeffs.conjugate())
.setTrans(m_trans)
.setLength(m_length)
.setShift(m_shift);
}
ConjugateReturnType adjoint() const
{ return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), !m_trans, m_actualVectors, m_shift); }
{
return conjugate().setTrans(!m_trans);
}
ConjugateReturnType inverse() const { return adjoint(); }
/** \internal */
template<typename DestType> void evalTo(DestType& dst) const
{
Index vecs = m_actualVectors;
Index vecs = m_length;
// FIXME find a way to pass this temporary if the user wants to
Matrix<Scalar, DestType::RowsAtCompileTime, 1,
AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> temp(rows());
@ -210,9 +223,9 @@ template<typename VectorsType, typename CoeffsType, int Side> class HouseholderS
template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
{
Matrix<Scalar,1,Dest::RowsAtCompileTime> temp(dst.rows());
for(Index k = 0; k < m_actualVectors; ++k)
for(Index k = 0; k < m_length; ++k)
{
Index actual_k = m_trans ? m_actualVectors-k-1 : k;
Index actual_k = m_trans ? m_length-k-1 : k;
dst.rightCols(rows()-m_shift-actual_k)
.applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
}
@ -222,9 +235,9 @@ template<typename VectorsType, typename CoeffsType, int Side> class HouseholderS
template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
{
Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.cols());
for(Index k = 0; k < m_actualVectors; ++k)
for(Index k = 0; k < m_length; ++k)
{
Index actual_k = m_trans ? k : m_actualVectors-k-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), &temp.coeffRef(0));
}
@ -250,40 +263,46 @@ template<typename VectorsType, typename CoeffsType, int Side> class HouseholderS
template<typename _VectorsType, typename _CoeffsType, int _Side> friend struct internal::hseq_side_dependent_impl;
HouseholderSequence& setTrans(bool trans)
{
m_trans = trans;
return *this;
}
HouseholderSequence& setLength(Index length)
{
m_length = length;
return *this;
}
HouseholderSequence& setShift(Index shift)
{
m_shift = shift;
return *this;
}
bool trans() const { return m_trans; }
Index length() const { return m_length; }
Index shift() const { return m_shift; }
protected:
typename VectorsType::Nested m_vectors;
typename CoeffsType::Nested m_coeffs;
bool m_trans;
Index m_actualVectors;
Index m_length;
Index m_shift;
};
template<typename VectorsType, typename CoeffsType>
HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h, bool trans=false)
HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h)
{
return HouseholderSequence<VectorsType,CoeffsType,OnTheLeft>(v, h, trans);
return HouseholderSequence<VectorsType,CoeffsType,OnTheLeft>(v, h);
}
template<typename VectorsType, typename CoeffsType>
HouseholderSequence<VectorsType,CoeffsType> householderSequence
(const VectorsType& v, const CoeffsType& h,
bool trans, typename VectorsType::Index actualVectors, typename VectorsType::Index shift)
HouseholderSequence<VectorsType,CoeffsType,OnTheRight> rightHouseholderSequence(const VectorsType& v, const CoeffsType& h)
{
return HouseholderSequence<VectorsType,CoeffsType,OnTheLeft>(v, h, trans, actualVectors, shift);
}
template<typename VectorsType, typename CoeffsType>
HouseholderSequence<VectorsType,CoeffsType,OnTheRight> rightHouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans=false)
{
return HouseholderSequence<VectorsType,CoeffsType,OnTheRight>(v, h, trans);
}
template<typename VectorsType, typename CoeffsType>
HouseholderSequence<VectorsType,CoeffsType,OnTheRight> rightHouseholderSequence
(const VectorsType& v, const CoeffsType& h, bool trans,
typename VectorsType::Index actualVectors, typename VectorsType::Index shift)
{
return HouseholderSequence<VectorsType,CoeffsType,OnTheRight>(v, h, trans, actualVectors, shift);
return HouseholderSequence<VectorsType,CoeffsType,OnTheRight>(v, h);
}
#endif // EIGEN_HOUSEHOLDER_SEQUENCE_H

13
Eigen/src/QR/ColPivHouseholderQR.h
View File

@ -483,13 +483,10 @@ struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
typename Rhs::PlainObject c(rhs());
// Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
c.applyOnTheLeft(householderSequence(
dec().matrixQR(),
dec().hCoeffs(),
true,
dec().nonzeroPivots(),
0
));
c.applyOnTheLeft(householderSequence(dec().matrixQR(), dec().hCoeffs())
.setTrans(true)
.setLength(dec().nonzeroPivots())
);
dec().matrixQR()
.topLeftCorner(nonzero_pivots, nonzero_pivots)
@ -517,7 +514,7 @@ typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHousehol
::householderQ() const
{
eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate(), false, m_nonzero_pivots, 0);
return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()).setLength(m_nonzero_pivots);
}
/** \return the column-pivoting Householder QR decomposition of \c *this.

5
Eigen/src/SVD/UpperBidiagonalization.h
View File

@ -87,8 +87,9 @@ template<typename _MatrixType> class UpperBidiagonalization
const HouseholderVSequenceType householderV() // const here gives nasty errors and i'm lazy
{
eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized.");
return HouseholderVSequenceType(m_householder, m_householder.const_derived().template diagonal<1>(),
false, m_householder.cols()-1, 1);
return HouseholderVSequenceType(m_householder, m_householder.const_derived().template diagonal<1>())
.setLength(m_householder.cols()-1)
.setShift(1);
}
protected:

10
test/householder.cpp
View File

@ -102,7 +102,12 @@ template<typename MatrixType> void householder(const MatrixType& m)
m2 = m1;
m2.block(shift,0,brows,cols) = qr.matrixQR();
HCoeffsVectorType hc = qr.hCoeffs().conjugate();
HouseholderSequence<MatrixType, HCoeffsVectorType> hseq(m2, hc, false, hc.size(), shift);
HouseholderSequence<MatrixType, HCoeffsVectorType> hseq(m2, hc);
hseq.setLength(hc.size()).setShift(shift);
VERIFY(hseq.trans() == false);
VERIFY(hseq.length() == hc.size());
VERIFY(hseq.shift() == shift);
MatrixType m5 = m2;
m5.block(shift,0,brows,cols).template triangularView<StrictlyLower>().setZero();
VERIFY_IS_APPROX(hseq * m5, m1); // test applying hseq directly
@ -112,7 +117,8 @@ template<typename MatrixType> void householder(const MatrixType& m)
// test householder sequence on the right with a shift
TMatrixType tm2 = m2.transpose();
HouseholderSequence<TMatrixType, HCoeffsVectorType, OnTheRight> rhseq(tm2, hc, false, hc.size(), shift);
HouseholderSequence<TMatrixType, HCoeffsVectorType, OnTheRight> rhseq(tm2, hc);
rhseq.setLength(hc.size()).setShift(shift);
VERIFY_IS_APPROX(rhseq * m5, m1); // test applying rhseq directly
m3 = rhseq;
VERIFY_IS_APPROX(m3 * m5, m1); // test evaluating rhseq to a dense matrix, then applying