diff options
author | Gael Guennebaud <g.gael@free.fr> | 2009-07-16 00:03:17 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2009-07-16 00:03:17 +0200 |
commit | 34490f1493f8111c131e471a3dc7f6fbe5687404 (patch) | |
tree | ed6a86263f90ffc965139ce2b91814d29c0c3311 /Eigen | |
parent | 97c9445c60f307f6a59a212872173c256adf2acd (diff) |
* bugfixes in Product, and test/product_selfadjoint
* speed up in the extraction of the matrix Q in Tridiagonalization
Diffstat (limited to 'Eigen')
-rw-r--r-- | Eigen/src/Core/BandMatrix.h | 2 | ||||
-rw-r--r-- | Eigen/src/Core/Product.h | 2 | ||||
-rw-r--r-- | Eigen/src/QR/SelfAdjointEigenSolver.h | 1 | ||||
-rw-r--r-- | Eigen/src/QR/Tridiagonalization.h | 30 |
4 files changed, 24 insertions, 11 deletions
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<BandMatrix<_Scalar,Rows,Cols,Supers,Sub * \sa class BandMatrix */ template<typename Scalar, int Size, int Options> -class TridiagonalMatrix : public BandMatrix<Scalar,Size,Size,1,Options&SelfAdjoint?0:1,Options|RowMajor> +class TridiagonalMatrix : public BandMatrix<Scalar,Size,Size,Options&SelfAdjoint?0:1,1,Options|RowMajor> { typedef BandMatrix<Scalar,Size,Size,1,Options&SelfAdjoint?0:1,Options|RowMajor> 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<Derived>::lazyAssign(const Product<Lhs,Rhs,CacheFrien } else { - lazyAssign(static_cast<const MatrixBase<Product<Lhs,Rhs,CacheFriendlyProduct> > &>(product)); + lazyAssign(Product<Lhs,Rhs,NormalProduct>(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<typename _MatrixType> class SelfAdjointEigenSolver typedef Matrix<RealScalar, MatrixType::ColsAtCompileTime, 1> RealVectorType; typedef Matrix<RealScalar, Dynamic, 1> RealVectorTypeX; typedef Tridiagonalization<MatrixType> 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<typename _MatrixType> class Tridiagonalization */ inline const MatrixType& packedMatrix(void) const { return m_matrix; } - MatrixType matrixQ(void) const; - MatrixType matrixT(void) const; + MatrixType matrixQ() const; + template<typename QDerived> void matrixQInPlace(MatrixBase<QDerived>* q) const; + MatrixType matrixT() const; const DiagonalReturnType diagonal(void) const; const SubDiagonalReturnType subDiagonal(void) const; @@ -270,21 +271,32 @@ template<typename MatrixType> typename Tridiagonalization<MatrixType>::MatrixType Tridiagonalization<MatrixType>::matrixQ(void) const { + MatrixType matQ; + matrixQInPlace(&matQ); + return matQ; +} + +template<typename MatrixType> +template<typename QDerived> +void Tridiagonalization<MatrixType>::matrixQInPlace(MatrixBase<QDerived>* q) const +{ + QDerived& matQ = q->derived(); int n = m_matrix.rows(); - MatrixType matQ = MatrixType::Identity(n,n); + matQ = MatrixType::Identity(n,n); + 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; - // 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<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; } - return matQ; } /** Performs a full decomposition in place */ @@ -303,7 +315,7 @@ void Tridiagonalization<MatrixType>::decomposeInPlace(MatrixType& mat, DiagonalT diag = tridiag.diagonal(); subdiag = tridiag.subDiagonal(); if (extractQ) - mat = tridiag.matrixQ(); + tridiag.matrixQInPlace(&mat); } } |