GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-09-15 02:43:14 +08:00
Code Issues Packages Projects Releases Wiki Activity

* bugfixes in Product, and test/product_selfadjoint

Browse Source 
* speed up in the extraction of the matrix Q in Tridiagonalization
This commit is contained in:
Gael Guennebaud 2009-07-16 00:03:17 +02:00
parent 97c9445c60
commit 34490f1493
5 changed files with 33 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Eigen/src/Core/BandMatrix.h
View File

@ -198,7 +198,7 @@ class BandMatrix : public MultiplierBase<BandMatrix<_Scalar,Rows,Cols,Supers,Sub
* \sa class BandMatrix * \sa class BandMatrix
*/ */
template<typename Scalar, int Size, int Options> 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; typedef BandMatrix<Scalar,Size,Size,1,Options&SelfAdjoint?0:1,Options|RowMajor> Base;
public: public:

2
Eigen/src/Core/Product.h
View File

@ -880,7 +880,7 @@ inline Derived& MatrixBase<Derived>::lazyAssign(const Product<Lhs,Rhs,CacheFrien
} }
else else
{ {
lazyAssign(static_cast<const MatrixBase<Product<Lhs,Rhs,CacheFriendlyProduct> > &>(product)); lazyAssign(Product<Lhs,Rhs,NormalProduct>(product.lhs(),product.rhs()));
} }
return derived(); return derived();
} }

1
Eigen/src/QR/SelfAdjointEigenSolver.h
View File

@ -50,6 +50,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
typedef Matrix<RealScalar, MatrixType::ColsAtCompileTime, 1> RealVectorType; typedef Matrix<RealScalar, MatrixType::ColsAtCompileTime, 1> RealVectorType;
typedef Matrix<RealScalar, Dynamic, 1> RealVectorTypeX; typedef Matrix<RealScalar, Dynamic, 1> RealVectorTypeX;
typedef Tridiagonalization<MatrixType> TridiagonalizationType; typedef Tridiagonalization<MatrixType> TridiagonalizationType;
// typedef typename TridiagonalizationType::TridiagonalMatrixType TridiagonalMatrixType;
SelfAdjointEigenSolver() SelfAdjointEigenSolver()
: m_eivec(int(Size), int(Size)), : m_eivec(int(Size), int(Size)),

30
Eigen/src/QR/Tridiagonalization.h
View File

@ -121,8 +121,9 @@ template<typename _MatrixType> class Tridiagonalization
*/ */
inline const MatrixType& packedMatrix(void) const { return m_matrix; } inline const MatrixType& packedMatrix(void) const { return m_matrix; }
MatrixType matrixQ(void) const; MatrixType matrixQ() const;
MatrixType matrixT(void) const; template<typename QDerived> void matrixQInPlace(MatrixBase<QDerived>* q) const;
MatrixType matrixT() const;
const DiagonalReturnType diagonal(void) const; const DiagonalReturnType diagonal(void) const;
const SubDiagonalReturnType subDiagonal(void) const; const SubDiagonalReturnType subDiagonal(void) const;
@ -270,21 +271,32 @@ template<typename MatrixType>
typename Tridiagonalization<MatrixType>::MatrixType typename Tridiagonalization<MatrixType>::MatrixType
Tridiagonalization<MatrixType>::matrixQ(void) const 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(); 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--) for (int i = n-2; i>=0; i--)
{ {
Scalar tmp = m_matrix.coeff(i+1,i); Scalar tmp = m_matrix.coeff(i+1,i);
m_matrix.const_cast_derived().coeffRef(i+1,i) = 1; 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 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();
matQ.corner(BottomRight,n-i-1,n-i-1) -= // rank one update, TODO ! make it works efficiently as expected
((m_hCoeffs.coeff(i) * m_matrix.col(i).end(n-i-1)) * for (int j=i+1;j<n;++j)
(m_matrix.col(i).end(n-i-1).adjoint() * matQ.corner(BottomRight,n-i-1,n-i-1)).lazy()).lazy(); 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; m_matrix.const_cast_derived().coeffRef(i+1,i) = tmp;
} }
return matQ;
} }
/** Performs a full decomposition in place */ /** Performs a full decomposition in place */
@ -303,7 +315,7 @@ void Tridiagonalization<MatrixType>::decomposeInPlace(MatrixType& mat, DiagonalT
diag = tridiag.diagonal(); diag = tridiag.diagonal();
subdiag = tridiag.subDiagonal(); subdiag = tridiag.subDiagonal();
if (extractQ) if (extractQ)
mat = tridiag.matrixQ(); tridiag.matrixQInPlace(&mat);
} }
} }

6
test/product_selfadjoint.cpp
View File

@ -77,12 +77,14 @@ template<typename MatrixType> void product_selfadjoint(const MatrixType& m)
m2.template selfadjointView<UpperTriangular>().rank2update(-r1.adjoint(),r2.adjoint()*s3,s1); m2.template selfadjointView<UpperTriangular>().rank2update(-r1.adjoint(),r2.adjoint()*s3,s1);
VERIFY_IS_APPROX(m2, (m1 + (-s3*s1) * (r1.adjoint() * r2 + r2.adjoint() * r1)).template triangularView<UpperTriangular>().toDense()); VERIFY_IS_APPROX(m2, (m1 + (-s3*s1) * (r1.adjoint() * r2 + r2.adjoint() * r1)).template triangularView<UpperTriangular>().toDense());
if (rows>1)
{
m2 = m1.template triangularView<LowerTriangular>(); m2 = m1.template triangularView<LowerTriangular>();
m2.block(1,1,rows-1,cols-1).template selfadjointView<LowerTriangular>().rank2update(v1.end(rows-1),v2.start(cols-1)); m2.block(1,1,rows-1,cols-1).template selfadjointView<LowerTriangular>().rank2update(v1.end(rows-1),v2.start(cols-1));
m3 = m1; 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(); 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<LowerTriangular>().toDense()); VERIFY_IS_APPROX(m2, m3.template triangularView<LowerTriangular>().toDense());
}
} }
void test_product_selfadjoint() void test_product_selfadjoint()
@ -93,7 +95,7 @@ void test_product_selfadjoint()
CALL_SUBTEST( product_selfadjoint(Matrix3d()) ); CALL_SUBTEST( product_selfadjoint(Matrix3d()) );
CALL_SUBTEST( product_selfadjoint(MatrixXcf(4, 4)) ); CALL_SUBTEST( product_selfadjoint(MatrixXcf(4, 4)) );
CALL_SUBTEST( product_selfadjoint(MatrixXcd(21,21)) ); 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<float,Dynamic,Dynamic,RowMajor>(17,17)) ); CALL_SUBTEST( product_selfadjoint(Matrix<float,Dynamic,Dynamic,RowMajor>(17,17)) );
CALL_SUBTEST( product_selfadjoint(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(19, 19)) ); CALL_SUBTEST( product_selfadjoint(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(19, 19)) );
} }