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

fix bug #107: SelfAdjointEigenSolver and RowMajor (and add unit test)

Browse Source
This commit is contained in:
Gael Guennebaud 2010-11-04 09:33:05 +01:00
parent 20fcef9656
commit 5a4f77716d
2 changed files with 7 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
View File

@ -336,7 +336,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
* "implicit symmetric QR step with Wilkinson shift"
*/
namespace internal {
template<typename RealScalar, typename Scalar, typename Index>
template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n);
}
@ -398,7 +398,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
while (start>0 && m_subdiag[start-1]!=0)
start--;
internal::tridiagonal_qr_step(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n);
internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n);
}
if (iter <= m_maxIterations)
@ -430,7 +430,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
}
namespace internal {
template<typename RealScalar, typename Scalar, typename Index>
template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n)
{
RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
@ -466,7 +466,8 @@ static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index sta
// apply the givens rotation to the unit matrix Q = Q * G
if (matrixQ)
{
Map<Matrix<Scalar,Dynamic,Dynamic> > q(matrixQ,n,n);
// FIXME if StorageOrder == RowMajor this operation is not very efficient
Map<Matrix<Scalar,Dynamic,Dynamic,StorageOrder> > q(matrixQ,n,n);
q.applyOnTheRight(k,k+1,rot);
}
}

2
test/eigensolver_selfadjoint.cpp
View File

@ -173,6 +173,8 @@ void test_eigensolver_selfadjoint()
CALL_SUBTEST_3( selfadjointeigensolver(MatrixXf(10,10)) );
CALL_SUBTEST_4( selfadjointeigensolver(MatrixXd(19,19)) );
CALL_SUBTEST_5( selfadjointeigensolver(MatrixXcd(17,17)) );
CALL_SUBTEST_9( selfadjointeigensolver(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(17,17)) );
// some trivial but implementation-wise tricky cases
CALL_SUBTEST_4( selfadjointeigensolver(MatrixXd(1,1)) );