aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/src/Core/MatrixBase.h8
-rw-r--r--Eigen/src/Householder/Householder.h34
-rw-r--r--Eigen/src/QR/QR.h23
-rw-r--r--Eigen/src/QR/Tridiagonalization.h80
-rw-r--r--test/eigensolver_selfadjoint.cpp2
-rw-r--r--test/householder.cpp31
-rw-r--r--test/qr.cpp2
7 files changed, 75 insertions, 105 deletions
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index 9e92c043f..1f4c6bf7a 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -786,17 +786,17 @@ template<typename Derived> class MatrixBase
////////// Householder module ///////////
- void makeHouseholderInPlace(RealScalar *tau, Scalar *beta);
+ void makeHouseholderInPlace(Scalar *tau, RealScalar *beta);
template<typename EssentialPart>
void makeHouseholder(EssentialPart *essential,
- RealScalar *tau, Scalar *beta) const;
+ Scalar *tau, RealScalar *beta) const;
template<typename EssentialPart>
void applyHouseholderOnTheLeft(const EssentialPart& essential,
- const RealScalar& tau,
+ const Scalar& tau,
Scalar* workspace);
template<typename EssentialPart>
void applyHouseholderOnTheRight(const EssentialPart& essential,
- const RealScalar& tau,
+ const Scalar& tau,
Scalar* workspace);
///////// Jacobi module /////////
diff --git a/Eigen/src/Householder/Householder.h b/Eigen/src/Householder/Householder.h
index a7bbe96ce..769ba3d43 100644
--- a/Eigen/src/Householder/Householder.h
+++ b/Eigen/src/Householder/Householder.h
@@ -42,7 +42,7 @@ void makeTrivialHouseholder(
}
template<typename Derived>
-void MatrixBase<Derived>::makeHouseholderInPlace(RealScalar *tau, Scalar *beta)
+void MatrixBase<Derived>::makeHouseholderInPlace(Scalar *tau, RealScalar *beta)
{
VectorBlock<Derived, ei_decrement_size<SizeAtCompileTime>::ret> essentialPart(derived(), 1, size()-1);
makeHouseholder(&essentialPart, tau, beta);
@@ -67,33 +67,35 @@ template<typename Derived>
template<typename EssentialPart>
void MatrixBase<Derived>::makeHouseholder(
EssentialPart *essential,
- RealScalar *tau,
- Scalar *beta) const
+ Scalar *tau,
+ RealScalar *beta) const
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart)
- RealScalar _squaredNorm = squaredNorm();
- Scalar c0;
- if(ei_abs2(coeff(0)) <= ei_abs2(precision<Scalar>()) * _squaredNorm)
+ VectorBlock<Derived, EssentialPart::SizeAtCompileTime> tail(derived(), 1, size()-1);
+
+ RealScalar tailSqNorm;
+ Scalar c0 = coeff(0);
+
+ if( (size()==1 || (tailSqNorm=tail.squaredNorm()) == RealScalar(0)) && ei_imag(c0)==RealScalar(0))
{
- c0 = ei_sqrt(_squaredNorm);
+ *tau = 0;
+ *beta = ei_real(c0);
}
else
{
- Scalar sign = coeff(0) / ei_abs(coeff(0));
- c0 = coeff(0) + sign * ei_sqrt(_squaredNorm);
+ *beta = ei_sqrt(ei_abs2(c0) + tailSqNorm);
+ if (ei_real(c0)>=0.)
+ *beta = -*beta;
+ *essential = tail / (c0 - *beta);
+ *tau = ei_conj((*beta - c0) / *beta);
}
- VectorBlock<Derived, EssentialPart::SizeAtCompileTime> tail(derived(), 1, size()-1);
- *essential = tail / c0;
- const RealScalar c0abs2 = ei_abs2(c0);
- *tau = RealScalar(2) * c0abs2 / (c0abs2 + _squaredNorm - ei_abs2(coeff(0)));
- *beta = coeff(0) - c0;
}
template<typename Derived>
template<typename EssentialPart>
void MatrixBase<Derived>::applyHouseholderOnTheLeft(
const EssentialPart& essential,
- const RealScalar& tau,
+ const Scalar& tau,
Scalar* workspace)
{
Map<Matrix<Scalar, 1, ColsAtCompileTime, PlainMatrixType::Options, 1, MaxColsAtCompileTime> > tmp(workspace,cols());
@@ -108,7 +110,7 @@ template<typename Derived>
template<typename EssentialPart>
void MatrixBase<Derived>::applyHouseholderOnTheRight(
const EssentialPart& essential,
- const RealScalar& tau,
+ const Scalar& tau,
Scalar* workspace)
{
Map<Matrix<Scalar, RowsAtCompileTime, 1, PlainMatrixType::Options, MaxRowsAtCompileTime, 1> > tmp(workspace,rows());
diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h
index 0c6142129..90e6f8132 100644
--- a/Eigen/src/QR/QR.h
+++ b/Eigen/src/QR/QR.h
@@ -53,7 +53,7 @@ template<typename MatrixType> class HouseholderQR
typedef typename MatrixType::RealScalar RealScalar;
typedef Block<MatrixType, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> MatrixRBlockType;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> MatrixTypeR;
- typedef Matrix<RealScalar, MinSizeAtCompileTime, 1> HCoeffsType;
+ typedef Matrix<Scalar, MinSizeAtCompileTime, 1> HCoeffsType;
/**
* \brief Default Constructor.
@@ -132,11 +132,13 @@ HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType&
int remainingRows = rows - k;
int remainingCols = cols - k -1;
- m_qr.col(k).end(remainingRows).makeHouseholderInPlace(&m_hCoeffs.coeffRef(k), &m_qr.coeffRef(k,k));
-
- if (remainingRows>1 && remainingCols>0)
- m_qr.corner(BottomRight, remainingRows, remainingCols)
- .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingRows-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1));
+ RealScalar beta;
+ m_qr.col(k).end(remainingRows).makeHouseholderInPlace(&m_hCoeffs.coeffRef(k), &beta);
+ m_qr.coeffRef(k,k) = beta;
+
+ // apply H to remaining part of m_qr from the left
+ m_qr.corner(BottomRight, remainingRows, remainingCols)
+ .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingRows-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1));
}
m_isInitialized = true;
return *this;
@@ -163,7 +165,7 @@ void HouseholderQR<MatrixType>::solve(
int remainingSize = rows-k;
result->corner(BottomRight, remainingSize, cols)
- .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), ei_real(m_hCoeffs.coeff(k)), &temp.coeffRef(0));
+ .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), m_hCoeffs.coeff(k), &temp.coeffRef(0));
}
const int rank = std::min(result->rows(), result->cols());
@@ -177,8 +179,8 @@ template<typename MatrixType>
MatrixType HouseholderQR<MatrixType>::matrixQ() const
{
ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
- // compute the product Q_0 Q_1 ... Q_n-1,
- // where Q_k is the k-th Householder transformation I - h_k v_k v_k'
+ // compute the product H'_0 H'_1 ... H'_n-1,
+ // where H_k is the k-th Householder transformation I - h_k v_k v_k'
// and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
int rows = m_qr.rows();
int cols = m_qr.cols();
@@ -187,9 +189,8 @@ MatrixType HouseholderQR<MatrixType>::matrixQ() const
for (int k = cols-1; k >= 0; k--)
{
int remainingSize = rows-k;
-
res.corner(BottomRight, remainingSize, cols-k)
- .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), ei_real(m_hCoeffs.coeff(k)), &temp.coeffRef(k));
+ .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), ei_conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k));
}
return res;
}
diff --git a/Eigen/src/QR/Tridiagonalization.h b/Eigen/src/QR/Tridiagonalization.h
index d98322fac..4cc7433c1 100644
--- a/Eigen/src/QR/Tridiagonalization.h
+++ b/Eigen/src/QR/Tridiagonalization.h
@@ -198,65 +198,29 @@ void Tridiagonalization<MatrixType>::_compute(MatrixType& matA, CoeffVectorType&
{
assert(matA.rows()==matA.cols());
int n = matA.rows();
- for (int i = 0; i<n-2; ++i)
+ Matrix<Scalar,1,Dynamic> aux(n);
+ for (int i = 0; i<n-1; ++i)
{
- // let's consider the vector v = i-th column starting at position i+1
-
- // start of the householder transformation
- // squared norm of the vector v skipping the first element
- RealScalar v1norm2 = matA.col(i).end(n-(i+2)).squaredNorm();
-
- // FIXME comparing against 1
- if (ei_isMuchSmallerThan(v1norm2,static_cast<Scalar>(1)))
- {
- hCoeffs.coeffRef(i) = 0.;
- }
- else
- {
- Scalar v0 = matA.col(i).coeff(i+1);
- RealScalar beta = ei_sqrt(ei_abs2(v0)+v1norm2);
- if (ei_real(v0)>=0.)
- beta = -beta;
- matA.col(i).end(n-(i+2)) *= (Scalar(1)/(v0-beta));
- matA.col(i).coeffRef(i+1) = beta;
- Scalar h = (beta - v0) / beta;
- // end of the householder transformation
+
+ int remainingSize = n-i-1;
+ RealScalar beta;
+ Scalar h;
+ matA.col(i).end(remainingSize).makeHouseholderInPlace(&h, &beta);
- // Apply similarity transformation to remaining columns,
- // i.e., A = H' A H where H = I - h v v' and v = matA.col(i).end(n-i-1)
- matA.col(i).coeffRef(i+1) = 1;
+ // Apply similarity transformation to remaining columns,
+ // i.e., A = H A H' where H = I - h v v' and v = matA.col(i).end(n-i-1)
+ matA.col(i).coeffRef(i+1) = 1;
- hCoeffs.end(n-i-1) = (matA.corner(BottomRight,n-i-1,n-i-1).template selfadjointView<LowerTriangular>()
- * (h * matA.col(i).end(n-i-1)));
+ hCoeffs.end(n-i-1) = (matA.corner(BottomRight,remainingSize,remainingSize).template selfadjointView<LowerTriangular>()
+ * (ei_conj(h) * matA.col(i).end(remainingSize)));
- hCoeffs.end(n-i-1) += (h*Scalar(-0.5)*(hCoeffs.end(n-i-1).dot(matA.col(i).end(n-i-1)))) * matA.col(i).end(n-i-1);
+ hCoeffs.end(n-i-1) += (ei_conj(h)*Scalar(-0.5)*(hCoeffs.end(remainingSize).dot(matA.col(i).end(remainingSize)))) * matA.col(i).end(n-i-1);
- matA.corner(BottomRight, n-i-1, n-i-1).template selfadjointView<LowerTriangular>()
- .rankUpdate(matA.col(i).end(n-i-1), hCoeffs.end(n-i-1), -1);
+ matA.corner(BottomRight, remainingSize, remainingSize).template selfadjointView<LowerTriangular>()
+ .rankUpdate(matA.col(i).end(remainingSize), hCoeffs.end(remainingSize), -1);
- // note: at that point matA(i+1,i+1) is the (i+1)-th element of the final diagonal
- // note: the sequence of the beta values leads to the subdiagonal entries
- matA.col(i).coeffRef(i+1) = beta;
-
- hCoeffs.coeffRef(i) = h;
- }
- }
- if (NumTraits<Scalar>::IsComplex)
- {
- // Householder transformation on the remaining single scalar
- int i = n-2;
- Scalar v0 = matA.col(i).coeff(i+1);
- RealScalar beta = ei_abs(v0);
- if (ei_real(v0)>=RealScalar(0))
- beta = -beta;
matA.col(i).coeffRef(i+1) = beta;
- // FIXME comparing against 1
- if(ei_isMuchSmallerThan(beta, Scalar(1))) hCoeffs.coeffRef(i) = Scalar(0);
- else hCoeffs.coeffRef(i) = (beta - v0) / beta;
- }
- else
- {
- hCoeffs.coeffRef(n-2) = 0;
+ hCoeffs.coeffRef(i) = h;
}
}
@@ -280,16 +244,8 @@ void Tridiagonalization<MatrixType>::matrixQInPlace(MatrixBase<QDerived>* q) con
Matrix<Scalar,1,Dynamic> aux(n);
for (int i = n-2; i>=0; i--)
{
- Scalar tmp = m_matrix.coeff(i+1,i);
- m_matrix.const_cast_derived().coeffRef(i+1,i) = 1;
-
- aux.end(n-i-1) = (m_hCoeffs.coeff(i) * m_matrix.col(i).end(n-i-1).adjoint() * matQ.corner(BottomRight,n-i-1,n-i-1)).lazy();
- // rank one update, TODO ! make it works efficiently as expected
- for (int j=i+1;j<n;++j)
- matQ.col(j).end(n-i-1) -= aux.coeff(j) * m_matrix.col(i).end(n-i-1);
-// matQ.corner(BottomRight,n-i-1,n-i-1) -= (m_matrix.col(i).end(n-i-1) * aux.end(n-i-1)).lazy();
-
- m_matrix.const_cast_derived().coeffRef(i+1,i) = tmp;
+ matQ.corner(BottomRight,n-i-1,n-i-1)
+ .applyHouseholderOnTheLeft(m_matrix.col(i).end(n-i-2), ei_conj(m_hCoeffs.coeff(i)), &aux.coeffRef(0,0));
}
}
diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp
index 1f10b131b..2571dbf43 100644
--- a/test/eigensolver_selfadjoint.cpp
+++ b/test/eigensolver_selfadjoint.cpp
@@ -120,8 +120,8 @@ void test_eigensolver_selfadjoint()
CALL_SUBTEST( selfadjointeigensolver(Matrix3f()) );
CALL_SUBTEST( selfadjointeigensolver(Matrix4d()) );
CALL_SUBTEST( selfadjointeigensolver(MatrixXf(10,10)) );
- CALL_SUBTEST( selfadjointeigensolver(MatrixXcd(17,17)) );
CALL_SUBTEST( selfadjointeigensolver(MatrixXd(19,19)) );
+ CALL_SUBTEST( selfadjointeigensolver(MatrixXcd(17,17)) );
// some trivial but implementation-wise tricky cases
CALL_SUBTEST( selfadjointeigensolver(MatrixXd(1,1)) );
diff --git a/test/householder.cpp b/test/householder.cpp
index d8b11c9f1..7d300899f 100644
--- a/test/householder.cpp
+++ b/test/householder.cpp
@@ -27,6 +27,8 @@
template<typename MatrixType> void householder(const MatrixType& m)
{
+ static bool even = true;
+ even = !even;
/* this test covers the following files:
Householder.h
*/
@@ -38,46 +40,55 @@ template<typename MatrixType> void householder(const MatrixType& m)
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
typedef Matrix<Scalar, ei_decrement_size<MatrixType::RowsAtCompileTime>::ret, 1> EssentialVectorType;
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;
+
+ Matrix<Scalar, EIGEN_ENUM_MAX(MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime), 1> _tmp(std::max(rows,cols));
+ Scalar* tmp = &_tmp.coeffRef(0,0);
- RealScalar beta;
+ Scalar beta;
+ RealScalar alpha;
EssentialVectorType essential;
VectorType v1 = VectorType::Random(rows), v2;
v2 = v1;
- v1.makeHouseholder(&essential, &beta);
- v1.applyHouseholderOnTheLeft(essential,beta);
-
+ v1.makeHouseholder(&essential, &beta, &alpha);
+ v1.applyHouseholderOnTheLeft(essential,beta,tmp);
VERIFY_IS_APPROX(v1.norm(), v2.norm());
VERIFY_IS_MUCH_SMALLER_THAN(v1.end(rows-1).norm(), v1.norm());
v1 = VectorType::Random(rows);
v2 = v1;
- v1.applyHouseholderOnTheLeft(essential,beta);
+ v1.applyHouseholderOnTheLeft(essential,beta,tmp);
VERIFY_IS_APPROX(v1.norm(), v2.norm());
MatrixType m1(rows, cols),
m2(rows, cols);
v1 = VectorType::Random(rows);
+ if(even) v1.end(rows-1).setZero();
m1.colwise() = v1;
m2 = m1;
- m1.col(0).makeHouseholder(&essential, &beta);
- m1.applyHouseholderOnTheLeft(essential,beta);
+ m1.col(0).makeHouseholder(&essential, &beta, &alpha);
+ m1.applyHouseholderOnTheLeft(essential,beta,tmp);
VERIFY_IS_APPROX(m1.norm(), m2.norm());
VERIFY_IS_MUCH_SMALLER_THAN(m1.block(1,0,rows-1,cols).norm(), m1.norm());
+ VERIFY_IS_MUCH_SMALLER_THAN(ei_imag(m1(0,0)), ei_real(m1(0,0)));
+ VERIFY_IS_APPROX(ei_real(m1(0,0)), alpha);
v1 = VectorType::Random(rows);
+ if(even) v1.end(rows-1).setZero();
SquareMatrixType m3(rows,rows), m4(rows,rows);
m3.rowwise() = v1.transpose();
m4 = m3;
- m3.row(0).makeHouseholder(&essential, &beta);
- m3.applyHouseholderOnTheRight(essential,beta);
+ m3.row(0).makeHouseholder(&essential, &beta, &alpha);
+ m3.applyHouseholderOnTheRight(essential,beta,tmp);
VERIFY_IS_APPROX(m3.norm(), m4.norm());
VERIFY_IS_MUCH_SMALLER_THAN(m3.block(0,1,rows,rows-1).norm(), m3.norm());
+ VERIFY_IS_MUCH_SMALLER_THAN(ei_imag(m3(0,0)), ei_real(m3(0,0)));
+ VERIFY_IS_APPROX(ei_real(m3(0,0)), alpha);
}
void test_householder()
{
- for(int i = 0; i < g_repeat; i++) {
+ for(int i = 0; i < 2*g_repeat; i++) {
CALL_SUBTEST( householder(Matrix<double,2,2>()) );
CALL_SUBTEST( householder(Matrix<float,2,3>()) );
CALL_SUBTEST( householder(Matrix<double,3,5>()) );
diff --git a/test/qr.cpp b/test/qr.cpp
index 99df1d89b..f004a36ca 100644
--- a/test/qr.cpp
+++ b/test/qr.cpp
@@ -93,8 +93,8 @@ void test_qr()
// FIXME : very weird bug here
// CALL_SUBTEST( qr(Matrix2f()) );
CALL_SUBTEST( qr(Matrix4d()) );
- CALL_SUBTEST( qr(MatrixXcd(17,7)) );
CALL_SUBTEST( qr(MatrixXf(47,40)) );
+ CALL_SUBTEST( qr(MatrixXcd(17,7)) );
}
for(int i = 0; i < g_repeat; i++) {
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-04 14:46:39 +0000