aboutsummaryrefslogtreecommitdiffhomepage
path: root/Eigen
diff options
context:
space:
mode:
authorGravatar Benoit Jacob <jacob.benoit.1@gmail.com>2009-12-02 11:11:09 -0500
committerGravatar Benoit Jacob <jacob.benoit.1@gmail.com>2009-12-02 11:11:09 -0500
commit3e73f6036c4f28ae8d11ae43641c213e608529e6 (patch)
treef9287eac5106ed6618f7b58b4dc8ab311b186bb2 /Eigen
parent3279e3934013d28b3870dd861eb64aec241a38b7 (diff)
* HouseholderSequence:
* be aware of number of actual householder vectors (optimization in non-full-rank case, no behavior change) * fix applyThisOnTheRight, it was using k instead of actual_k * QR: rename matrixQ() to householderQ() where applicable
Diffstat (limited to 'Eigen')
-rw-r--r--Eigen/src/Householder/HouseholderSequence.h37
-rw-r--r--Eigen/src/QR/ColPivHouseholderQR.h10
-rw-r--r--Eigen/src/QR/HouseholderQR.h13
3 files changed, 32 insertions, 28 deletions
diff --git a/Eigen/src/Householder/HouseholderSequence.h b/Eigen/src/Householder/HouseholderSequence.h
index afeb758b3..1159f03ed 100644
--- a/Eigen/src/Householder/HouseholderSequence.h
+++ b/Eigen/src/Householder/HouseholderSequence.h
@@ -69,31 +69,35 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence
typedef HouseholderSequence<VectorsType,
typename ei_meta_if<NumTraits<Scalar>::IsComplex,
- typename ei_cleantype<typename CoeffsType::ConjugateReturnType>::type,
+ NestByValue<typename ei_cleantype<typename CoeffsType::ConjugateReturnType>::type >,
CoeffsType>::ret> ConjugateReturnType;
HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans = false)
- : m_vectors(v), m_coeffs(h), m_trans(trans)
+ : m_vectors(v), m_coeffs(h), m_trans(trans), m_actualVectors(v.diagonalSize())
+ {}
+
+ HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans, int actualVectors)
+ : m_vectors(v), m_coeffs(h), m_trans(trans), m_actualVectors(actualVectors)
{}
int rows() const { return m_vectors.rows(); }
int cols() const { return m_vectors.rows(); }
HouseholderSequence transpose() const
- { return HouseholderSequence(m_vectors, m_coeffs, !m_trans); }
+ { return HouseholderSequence(m_vectors, m_coeffs, !m_trans, m_actualVectors); }
ConjugateReturnType conjugate() const
- { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans); }
+ { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans, m_actualVectors); }
ConjugateReturnType adjoint() const
- { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), !m_trans); }
+ { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), !m_trans, m_actualVectors); }
ConjugateReturnType inverse() const { return adjoint(); }
/** \internal */
template<typename DestType> void evalTo(DestType& dst) const
{
- int vecs = std::min(m_vectors.cols(),m_vectors.rows());
+ int vecs = m_actualVectors;
int length = m_vectors.rows();
dst.setIdentity();
Matrix<Scalar,1,DestType::RowsAtCompileTime> temp(dst.rows());
@@ -111,22 +115,22 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence
/** \internal */
template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
{
- int vecs = std::min(m_vectors.cols(),m_vectors.rows()); // number of householder vectors
- int length = m_vectors.rows(); // size of the largest householder vector
- Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.rows());
+ int vecs = m_actualVectors; // number of householder vectors
+ int length = m_vectors.rows(); // size of the largest householder vector
+ Matrix<Scalar,1,Dest::RowsAtCompileTime> temp(dst.rows());
for(int k = 0; k < vecs; ++k)
{
int actual_k = m_trans ? vecs-k-1 : k;
- dst.corner(BottomRight, dst.rows(), length-k)
- .applyHouseholderOnTheRight(m_vectors.col(k).end(length-k-1), m_coeffs.coeff(k), &temp.coeffRef(0));
+ dst.corner(BottomRight, dst.rows(), length-actual_k)
+ .applyHouseholderOnTheRight(m_vectors.col(actual_k).end(length-actual_k-1), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
}
}
/** \internal */
template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
{
- int vecs = std::min(m_vectors.cols(),m_vectors.rows()); // number of householder vectors
- int length = m_vectors.rows(); // size of the largest householder vector
+ int vecs = m_actualVectors; // number of householder vectors
+ int length = m_vectors.rows(); // size of the largest householder vector
Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.cols());
for(int k = 0; k < vecs; ++k)
{
@@ -156,6 +160,7 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence
typename VectorsType::Nested m_vectors;
typename CoeffsType::Nested m_coeffs;
bool m_trans;
+ int m_actualVectors;
private:
HouseholderSequence& operator=(const HouseholderSequence&);
@@ -167,4 +172,10 @@ HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsTyp
return HouseholderSequence<VectorsType,CoeffsType>(v, h, trans);
}
+template<typename VectorsType, typename CoeffsType>
+HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h, bool trans, int actualVectors)
+{
+ return HouseholderSequence<VectorsType,CoeffsType>(v, h, trans, actualVectors);
+}
+
#endif // EIGEN_HOUSEHOLDER_SEQUENCE_H
diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h
index 9c1e56dc6..1b2472a99 100644
--- a/Eigen/src/QR/ColPivHouseholderQR.h
+++ b/Eigen/src/QR/ColPivHouseholderQR.h
@@ -105,7 +105,7 @@ template<typename _MatrixType> class ColPivHouseholderQR
return ei_solve_retval<ColPivHouseholderQR, Rhs>(*this, b.derived());
}
- HouseholderSequenceType matrixQ(void) const;
+ HouseholderSequenceType householderQ(void) const;
/** \returns a reference to the matrix where the Householder QR decomposition is stored
*/
@@ -447,7 +447,8 @@ struct ei_solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
c.applyOnTheLeft(householderSequence(
dec().matrixQR(),
dec().hCoeffs(),
- true
+ true,
+ dec().nonzeroPivots()
));
dec().matrixQR()
@@ -470,10 +471,11 @@ struct ei_solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
/** \returns the matrix Q as a sequence of householder transformations */
template<typename MatrixType>
-typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>::matrixQ() const
+typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
+ ::householderQ() const
{
ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
- return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate(), false);
+ return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate(), false, m_nonzero_pivots);
}
#endif // EIGEN_HIDE_HEAVY_CODE
diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h
index 8d842d129..95496b943 100644
--- a/Eigen/src/QR/HouseholderQR.h
+++ b/Eigen/src/QR/HouseholderQR.h
@@ -104,11 +104,10 @@ template<typename _MatrixType> class HouseholderQR
ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
return ei_solve_retval<HouseholderQR, Rhs>(*this, b.derived());
}
-
- MatrixQType matrixQ() const;
- HouseholderSequenceType matrixQAsHouseholderSequence() const
+ HouseholderSequenceType householderQ() const
{
+ ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
}
@@ -240,14 +239,6 @@ struct ei_solve_retval<HouseholderQR<_MatrixType>, Rhs>
}
};
-/** \returns the matrix Q */
-template<typename MatrixType>
-typename HouseholderQR<MatrixType>::MatrixQType HouseholderQR<MatrixType>::matrixQ() const
-{
- ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
- return matrixQAsHouseholderSequence();
-}
-
#endif // EIGEN_HIDE_HEAVY_CODE
/** \return the Householder QR decomposition of \c *this.