aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/src/Eigenvalues/HessenbergDecomposition.h2
-rw-r--r--Eigen/src/Householder/Householder.h8
-rw-r--r--Eigen/src/Householder/HouseholderSequence.h36
-rw-r--r--Eigen/src/QR/ColPivHouseholderQR.h6
-rw-r--r--Eigen/src/QR/CompleteOrthogonalDecomposition.h7
-rw-r--r--Eigen/src/QR/HouseholderQR.h6
-rw-r--r--Eigen/src/SVD/UpperBidiagonalization.h2
-rw-r--r--test/householder.cpp29
-rw-r--r--test/qr.cpp2
9 files changed, 60 insertions, 38 deletions
diff --git a/Eigen/src/Eigenvalues/HessenbergDecomposition.h b/Eigen/src/Eigenvalues/HessenbergDecomposition.h
index f647f69b0..d947dac4e 100644
--- a/Eigen/src/Eigenvalues/HessenbergDecomposition.h
+++ b/Eigen/src/Eigenvalues/HessenbergDecomposition.h
@@ -315,7 +315,7 @@ void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVector
// A = A H'
matA.rightCols(remainingSize)
- .applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1).conjugate(), numext::conj(h), &temp.coeffRef(0));
+ .applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1), numext::conj(h), &temp.coeffRef(0));
}
}
diff --git a/Eigen/src/Householder/Householder.h b/Eigen/src/Householder/Householder.h
index 80de2c305..a5f336d18 100644
--- a/Eigen/src/Householder/Householder.h
+++ b/Eigen/src/Householder/Householder.h
@@ -103,7 +103,7 @@ void MatrixBase<Derived>::makeHouseholder(
* \param essential the essential part of the vector \c v
* \param tau the scaling factor of the Householder transformation
* \param workspace a pointer to working space with at least
- * this->cols() * essential.size() entries
+ * this->cols() entries
*
* \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(),
* MatrixBase::applyHouseholderOnTheRight()
@@ -140,7 +140,7 @@ void MatrixBase<Derived>::applyHouseholderOnTheLeft(
* \param essential the essential part of the vector \c v
* \param tau the scaling factor of the Householder transformation
* \param workspace a pointer to working space with at least
- * this->cols() * essential.size() entries
+ * this->rows() entries
*
* \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(),
* MatrixBase::applyHouseholderOnTheLeft()
@@ -160,10 +160,10 @@ void MatrixBase<Derived>::applyHouseholderOnTheRight(
{
Map<typename internal::plain_col_type<PlainObject>::type> tmp(workspace,rows());
Block<Derived, Derived::RowsAtCompileTime, EssentialPart::SizeAtCompileTime> right(derived(), 0, 1, rows(), cols()-1);
- tmp.noalias() = right * essential.conjugate();
+ tmp.noalias() = right * essential;
tmp += this->col(0);
this->col(0) -= tau * tmp;
- right.noalias() -= tau * tmp * essential.transpose();
+ right.noalias() -= tau * tmp * essential.adjoint();
}
}
diff --git a/Eigen/src/Householder/HouseholderSequence.h b/Eigen/src/Householder/HouseholderSequence.h
index 3ce0a693d..13030c664 100644
--- a/Eigen/src/Householder/HouseholderSequence.h
+++ b/Eigen/src/Householder/HouseholderSequence.h
@@ -140,6 +140,22 @@ template<typename VectorsType, typename CoeffsType, int Side> class HouseholderS
Side
> ConjugateReturnType;
+ typedef HouseholderSequence<
+ VectorsType,
+ typename internal::conditional<NumTraits<Scalar>::IsComplex,
+ typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type,
+ CoeffsType>::type,
+ Side
+ > AdjointReturnType;
+
+ typedef HouseholderSequence<
+ typename internal::conditional<NumTraits<Scalar>::IsComplex,
+ typename internal::remove_all<typename VectorsType::ConjugateReturnType>::type,
+ VectorsType>::type,
+ CoeffsType,
+ Side
+ > TransposeReturnType;
+
/** \brief Constructor.
* \param[in] v %Matrix containing the essential parts of the Householder vectors
* \param[in] h Vector containing the Householder coefficients
@@ -206,9 +222,12 @@ template<typename VectorsType, typename CoeffsType, int Side> class HouseholderS
}
/** \brief %Transpose of the Householder sequence. */
- HouseholderSequence transpose() const
+ TransposeReturnType transpose() const
{
- return HouseholderSequence(*this).setTrans(!m_trans);
+ return TransposeReturnType(m_vectors.conjugate(), m_coeffs)
+ .setTrans(!m_trans)
+ .setLength(m_length)
+ .setShift(m_shift);
}
/** \brief Complex conjugate of the Householder sequence. */
@@ -221,13 +240,16 @@ template<typename VectorsType, typename CoeffsType, int Side> class HouseholderS
}
/** \brief Adjoint (conjugate transpose) of the Householder sequence. */
- ConjugateReturnType adjoint() const
+ AdjointReturnType adjoint() const
{
- return conjugate().setTrans(!m_trans);
+ return AdjointReturnType(m_vectors, m_coeffs.conjugate())
+ .setTrans(!m_trans)
+ .setLength(m_length)
+ .setShift(m_shift);
}
/** \brief Inverse of the Householder sequence (equals the adjoint). */
- ConjugateReturnType inverse() const { return adjoint(); }
+ AdjointReturnType inverse() const { return adjoint(); }
/** \internal */
template<typename DestType> inline void evalTo(DestType& dst) const
@@ -273,10 +295,10 @@ template<typename VectorsType, typename CoeffsType, int Side> class HouseholderS
Index cornerSize = rows() - k - m_shift;
if(m_trans)
dst.bottomRightCorner(cornerSize, cornerSize)
- .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
+ .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data());
else
dst.bottomRightCorner(cornerSize, cornerSize)
- .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
+ .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data());
}
}
}
diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h
index ed47b05e3..1faa3442e 100644
--- a/Eigen/src/QR/ColPivHouseholderQR.h
+++ b/Eigen/src/QR/ColPivHouseholderQR.h
@@ -595,11 +595,7 @@ void ColPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &
typename RhsType::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(m_qr, m_hCoeffs)
- .setLength(nonzero_pivots)
- .transpose()
- );
+ c.applyOnTheLeft(householderQ().setLength(nonzero_pivots).adjoint() );
m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots)
.template triangularView<Upper>()
diff --git a/Eigen/src/QR/CompleteOrthogonalDecomposition.h b/Eigen/src/QR/CompleteOrthogonalDecomposition.h
index 880becb25..03017a375 100644
--- a/Eigen/src/QR/CompleteOrthogonalDecomposition.h
+++ b/Eigen/src/QR/CompleteOrthogonalDecomposition.h
@@ -452,7 +452,7 @@ void CompleteOrthogonalDecomposition<MatrixType>::computeInPlace()
// Apply Z(k) to the first k rows of X_k
m_cpqr.m_qr.topRightCorner(k, cols - rank + 1)
.applyHouseholderOnTheRight(
- m_cpqr.m_qr.row(k).tail(cols - rank).transpose(), m_zCoeffs(k),
+ m_cpqr.m_qr.row(k).tail(cols - rank).adjoint(), m_zCoeffs(k),
&m_temp(0));
}
if (k != rank - 1) {
@@ -500,11 +500,8 @@ void CompleteOrthogonalDecomposition<_MatrixType>::_solve_impl(
}
// Compute c = Q^* * rhs
- // Note that the matrix Q = H_0^* H_1^*... so its inverse is
- // Q^* = (H_0 H_1 ...)^T
typename RhsType::PlainObject c(rhs);
- c.applyOnTheLeft(
- householderSequence(matrixQTZ(), hCoeffs()).setLength(rank).transpose());
+ c.applyOnTheLeft(matrixQ().setLength(rank).adjoint());
// Solve T z = c(1:rank, :)
dst.topRows(rank) = matrixT()
diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h
index 762b21c36..13c02685d 100644
--- a/Eigen/src/QR/HouseholderQR.h
+++ b/Eigen/src/QR/HouseholderQR.h
@@ -353,11 +353,7 @@ void HouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) c
typename RhsType::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(
- m_qr.leftCols(rank),
- m_hCoeffs.head(rank)).transpose()
- );
+ c.applyOnTheLeft(householderQ().setLength(rank).adjoint() );
m_qr.topLeftCorner(rank, rank)
.template triangularView<Upper>()
diff --git a/Eigen/src/SVD/UpperBidiagonalization.h b/Eigen/src/SVD/UpperBidiagonalization.h
index 0526ac931..997defc47 100644
--- a/Eigen/src/SVD/UpperBidiagonalization.h
+++ b/Eigen/src/SVD/UpperBidiagonalization.h
@@ -127,7 +127,7 @@ void upperbidiagonalization_inplace_unblocked(MatrixType& mat,
.makeHouseholderInPlace(mat.coeffRef(k,k+1), upper_diagonal[k]);
// apply householder transform to remaining part of mat on the left
mat.bottomRightCorner(remainingRows-1, remainingCols)
- .applyHouseholderOnTheRight(mat.row(k).tail(remainingCols-1).transpose(), mat.coeff(k,k+1), tempData);
+ .applyHouseholderOnTheRight(mat.row(k).tail(remainingCols-1).adjoint(), mat.coeff(k,k+1), tempData);
}
}
diff --git a/test/householder.cpp b/test/householder.cpp
index c5f6b5e4f..78d44d48a 100644
--- a/test/householder.cpp
+++ b/test/householder.cpp
@@ -49,6 +49,17 @@ template<typename MatrixType> void householder(const MatrixType& m)
v1.applyHouseholderOnTheLeft(essential,beta,tmp);
VERIFY_IS_APPROX(v1.norm(), v2.norm());
+ // reconstruct householder matrix:
+ SquareMatrixType id, H1, H2;
+ id.setIdentity(rows, rows);
+ H1 = H2 = id;
+ VectorType vv(rows);
+ vv << Scalar(1), essential;
+ H1.applyHouseholderOnTheLeft(essential, beta, tmp);
+ H2.applyHouseholderOnTheRight(essential, beta, tmp);
+ VERIFY_IS_APPROX(H1, H2);
+ VERIFY_IS_APPROX(H1, id - beta * vv*vv.adjoint());
+
MatrixType m1(rows, cols),
m2(rows, cols);
@@ -69,7 +80,7 @@ template<typename MatrixType> void householder(const MatrixType& m)
m3.rowwise() = v1.transpose();
m4 = m3;
m3.row(0).makeHouseholder(essential, beta, alpha);
- m3.applyHouseholderOnTheRight(essential,beta,tmp);
+ m3.applyHouseholderOnTheRight(essential.conjugate(),beta,tmp);
VERIFY_IS_APPROX(m3.norm(), m4.norm());
if(rows>=2) VERIFY_IS_MUCH_SMALLER_THAN(m3.block(0,1,rows,rows-1).norm(), m3.norm());
VERIFY_IS_MUCH_SMALLER_THAN(numext::imag(m3(0,0)), numext::real(m3(0,0)));
@@ -104,14 +115,14 @@ template<typename MatrixType> void householder(const MatrixType& m)
VERIFY_IS_APPROX(hseq_mat.adjoint(), hseq_mat_adj);
VERIFY_IS_APPROX(hseq_mat.conjugate(), hseq_mat_conj);
VERIFY_IS_APPROX(hseq_mat.transpose(), hseq_mat_trans);
- VERIFY_IS_APPROX(hseq_mat * m6, hseq_mat * m6);
- VERIFY_IS_APPROX(hseq_mat.adjoint() * m6, hseq_mat_adj * m6);
- VERIFY_IS_APPROX(hseq_mat.conjugate() * m6, hseq_mat_conj * m6);
- VERIFY_IS_APPROX(hseq_mat.transpose() * m6, hseq_mat_trans * m6);
- VERIFY_IS_APPROX(m6 * hseq_mat, m6 * hseq_mat);
- VERIFY_IS_APPROX(m6 * hseq_mat.adjoint(), m6 * hseq_mat_adj);
- VERIFY_IS_APPROX(m6 * hseq_mat.conjugate(), m6 * hseq_mat_conj);
- VERIFY_IS_APPROX(m6 * hseq_mat.transpose(), m6 * hseq_mat_trans);
+ VERIFY_IS_APPROX(hseq * m6, hseq_mat * m6);
+ VERIFY_IS_APPROX(hseq.adjoint() * m6, hseq_mat_adj * m6);
+ VERIFY_IS_APPROX(hseq.conjugate() * m6, hseq_mat_conj * m6);
+ VERIFY_IS_APPROX(hseq.transpose() * m6, hseq_mat_trans * m6);
+ VERIFY_IS_APPROX(m6 * hseq, m6 * hseq_mat);
+ VERIFY_IS_APPROX(m6 * hseq.adjoint(), m6 * hseq_mat_adj);
+ VERIFY_IS_APPROX(m6 * hseq.conjugate(), m6 * hseq_mat_conj);
+ VERIFY_IS_APPROX(m6 * hseq.transpose(), m6 * hseq_mat_trans);
// test householder sequence on the right with a shift
diff --git a/test/qr.cpp b/test/qr.cpp
index dfcc1e8f9..00ef99ef9 100644
--- a/test/qr.cpp
+++ b/test/qr.cpp
@@ -85,7 +85,7 @@ template<typename MatrixType> void qr_invertible()
qr.compute(m1);
VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant());
// This test is tricky if the determinant becomes too small.
- // Since we generate random numbers with magnitude rrange [0,1], the average determinant is 0.5^size
+ // Since we generate random numbers with magnitude range [0,1], the average determinant is 0.5^size
VERIFY_IS_MUCH_SMALLER_THAN( abs(absdet-qr.absDeterminant()), numext::maxi(RealScalar(pow(0.5,size)),numext::maxi<RealScalar>(abs(absdet),abs(qr.absDeterminant()))) );
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-29 08:24:41 +0000