diff --git a/Eigen/src/Core/BandMatrix.h b/Eigen/src/Core/BandMatrix.h index 8136724dd..7a37b3624 100644 --- a/Eigen/src/Core/BandMatrix.h +++ b/Eigen/src/Core/BandMatrix.h @@ -198,7 +198,7 @@ class BandMatrix : public MultiplierBase -class TridiagonalMatrix : public BandMatrix +class TridiagonalMatrix : public BandMatrix { typedef BandMatrix Base; public: diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index 6feede458..93e318035 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -880,7 +880,7 @@ inline Derived& MatrixBase::lazyAssign(const Product > &>(product)); + lazyAssign(Product(product.lhs(),product.rhs())); } return derived(); } diff --git a/Eigen/src/QR/SelfAdjointEigenSolver.h b/Eigen/src/QR/SelfAdjointEigenSolver.h index ff6b98aa1..b003fb4d7 100644 --- a/Eigen/src/QR/SelfAdjointEigenSolver.h +++ b/Eigen/src/QR/SelfAdjointEigenSolver.h @@ -50,6 +50,7 @@ template class SelfAdjointEigenSolver typedef Matrix RealVectorType; typedef Matrix RealVectorTypeX; typedef Tridiagonalization TridiagonalizationType; +// typedef typename TridiagonalizationType::TridiagonalMatrixType TridiagonalMatrixType; SelfAdjointEigenSolver() : m_eivec(int(Size), int(Size)), diff --git a/Eigen/src/QR/Tridiagonalization.h b/Eigen/src/QR/Tridiagonalization.h index 8051f7a55..cf737bd30 100644 --- a/Eigen/src/QR/Tridiagonalization.h +++ b/Eigen/src/QR/Tridiagonalization.h @@ -121,8 +121,9 @@ template class Tridiagonalization */ inline const MatrixType& packedMatrix(void) const { return m_matrix; } - MatrixType matrixQ(void) const; - MatrixType matrixT(void) const; + MatrixType matrixQ() const; + template void matrixQInPlace(MatrixBase* q) const; + MatrixType matrixT() const; const DiagonalReturnType diagonal(void) const; const SubDiagonalReturnType subDiagonal(void) const; @@ -270,21 +271,32 @@ template typename Tridiagonalization::MatrixType Tridiagonalization::matrixQ(void) const { + MatrixType matQ; + matrixQInPlace(&matQ); + return matQ; +} + +template +template +void Tridiagonalization::matrixQInPlace(MatrixBase* q) const +{ + QDerived& matQ = q->derived(); int n = m_matrix.rows(); - MatrixType matQ = MatrixType::Identity(n,n); + matQ = MatrixType::Identity(n,n); + Matrix 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; - // TODO this product could be optimized by processing the submatrix per panel of at least 4 columns - matQ.corner(BottomRight,n-i-1,n-i-1) -= - ((m_hCoeffs.coeff(i) * m_matrix.col(i).end(n-i-1)) * - (m_matrix.col(i).end(n-i-1).adjoint() * matQ.corner(BottomRight,n-i-1,n-i-1)).lazy()).lazy(); + 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::decomposeInPlace(MatrixType& mat, DiagonalT diag = tridiag.diagonal(); subdiag = tridiag.subDiagonal(); if (extractQ) - mat = tridiag.matrixQ(); + tridiag.matrixQInPlace(&mat); } } diff --git a/test/product_selfadjoint.cpp b/test/product_selfadjoint.cpp index 297bab1a9..29fbf11bf 100644 --- a/test/product_selfadjoint.cpp +++ b/test/product_selfadjoint.cpp @@ -77,12 +77,14 @@ template void product_selfadjoint(const MatrixType& m) m2.template selfadjointView().rank2update(-r1.adjoint(),r2.adjoint()*s3,s1); VERIFY_IS_APPROX(m2, (m1 + (-s3*s1) * (r1.adjoint() * r2 + r2.adjoint() * r1)).template triangularView().toDense()); - m2 = m1.template triangularView(); - m2.block(1,1,rows-1,cols-1).template selfadjointView().rank2update(v1.end(rows-1),v2.start(cols-1)); - m3 = m1; - m3.block(1,1,rows-1,cols-1) += v1.end(rows-1) * v2.start(cols-1).adjoint()+ v2.start(cols-1) * v1.end(rows-1).adjoint(); - VERIFY_IS_APPROX(m2, m3.template triangularView().toDense()); - + if (rows>1) + { + m2 = m1.template triangularView(); + m2.block(1,1,rows-1,cols-1).template selfadjointView().rank2update(v1.end(rows-1),v2.start(cols-1)); + m3 = m1; + m3.block(1,1,rows-1,cols-1) += v1.end(rows-1) * v2.start(cols-1).adjoint()+ v2.start(cols-1) * v1.end(rows-1).adjoint(); + VERIFY_IS_APPROX(m2, m3.template triangularView().toDense()); + } } void test_product_selfadjoint() @@ -93,7 +95,7 @@ void test_product_selfadjoint() CALL_SUBTEST( product_selfadjoint(Matrix3d()) ); CALL_SUBTEST( product_selfadjoint(MatrixXcf(4, 4)) ); CALL_SUBTEST( product_selfadjoint(MatrixXcd(21,21)) ); - CALL_SUBTEST( product_selfadjoint(MatrixXd(4,4)) ); + CALL_SUBTEST( product_selfadjoint(MatrixXd(14,14)) ); CALL_SUBTEST( product_selfadjoint(Matrix(17,17)) ); CALL_SUBTEST( product_selfadjoint(Matrix,Dynamic,Dynamic,RowMajor>(19, 19)) ); }