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

Fix segfault and bug with equal eivals in matrix power (bug #614).

Browse Source
This commit is contained in:
Chen-Pang He 2013-06-24 13:58:51 +01:00
parent 1330ca611b
commit ee8a28fb85
2 changed files with 21 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
View File

@ -294,6 +294,7 @@ MatrixPowerTriangularAtomic<MatrixType>::MatrixPowerTriangularAtomic(const Matri
template<typename MatrixType> template<typename MatrixType>
void MatrixPowerTriangularAtomic<MatrixType>::compute(MatrixType& res, RealScalar p) const void MatrixPowerTriangularAtomic<MatrixType>::compute(MatrixType& res, RealScalar p) const
{ {
res.resizeLike(m_A);
switch (m_A.rows()) { switch (m_A.rows()) {
case 0: case 0:
break; break;
@ -320,6 +321,7 @@ void MatrixPowerTriangularAtomic<MatrixType>::computePade(int degree, const Matr
res += MatrixType::Identity(IminusT.rows(), IminusT.cols()); res += MatrixType::Identity(IminusT.rows(), IminusT.cols());
} }
// this function assumes that res has the correct size (see bug 614)
template<typename MatrixType> template<typename MatrixType>
void MatrixPowerTriangularAtomic<MatrixType>::compute2x2(MatrixType& res, RealScalar p) const void MatrixPowerTriangularAtomic<MatrixType>::compute2x2(MatrixType& res, RealScalar p) const
{ {
@ -332,17 +334,18 @@ void MatrixPowerTriangularAtomic<MatrixType>::compute2x2(MatrixType& res, RealSc
for (Index i=1; i < m_A.cols(); ++i) { for (Index i=1; i < m_A.cols(); ++i) {
res.coeffRef(i,i) = pow(m_A.coeff(i,i), p); res.coeffRef(i,i) = pow(m_A.coeff(i,i), p);
if (m_A.coeff(i-1,i-1) == m_A.coeff(i,i)) { if (m_A.coeff(i-1,i-1) == m_A.coeff(i,i)) {
res.coeffRef(i-1,i) = p * pow(m_A.coeff(i-1,i), p-1); res.coeffRef(i-1,i) = p * pow(m_A.coeff(i,i), p-1);
} }
else if (2*abs(m_A.coeff(i-1,i-1)) < abs(m_A.coeff(i,i)) || 2*abs(m_A.coeff(i,i)) < abs(m_A.coeff(i-1,i-1))) { else if (2*abs(m_A.coeff(i-1,i-1)) < abs(m_A.coeff(i,i)) || 2*abs(m_A.coeff(i,i)) < abs(m_A.coeff(i-1,i-1))) {
res.coeffRef(i-1,i) = m_A.coeff(i-1,i) * (res.coeff(i,i)-res.coeff(i-1,i-1)) / (m_A.coeff(i,i)-m_A.coeff(i-1,i-1)); res.coeffRef(i-1,i) = (res.coeff(i,i)-res.coeff(i-1,i-1)) / (m_A.coeff(i,i)-m_A.coeff(i-1,i-1));
} }
else { else {
int unwindingNumber = std::ceil((numext::imag(logTdiag[i]-logTdiag[i-1]) - M_PI) / (2*M_PI)); int unwindingNumber = std::ceil((numext::imag(logTdiag[i]-logTdiag[i-1]) - M_PI) / (2*M_PI));
Scalar w = internal::matrix_power_unwinder<Scalar>::run(m_A.coeff(i,i), m_A.coeff(i-1,i-1), unwindingNumber); Scalar w = internal::matrix_power_unwinder<Scalar>::run(m_A.coeff(i,i), m_A.coeff(i-1,i-1), unwindingNumber);
res.coeffRef(i-1,i) = m_A.coeff(i-1,i) * RealScalar(2) * std::exp(RealScalar(0.5)*p*(logTdiag[i]+logTdiag[i-1])) * res.coeffRef(i-1,i) = RealScalar(2) * std::exp(RealScalar(0.5)*p*(logTdiag[i]+logTdiag[i-1])) * std::sinh(p * w)
std::sinh(p * w) / (m_A.coeff(i,i) - m_A.coeff(i-1,i-1)); / (m_A.coeff(i,i) - m_A.coeff(i-1,i-1));
} }
res.coeffRef(i-1,i) *= m_A.coeff(i-1,i);
} }
} }
@ -360,6 +363,19 @@ void MatrixPowerTriangularAtomic<MatrixType>::computeBig(MatrixType& res, RealSc
int degree, degree2, numberOfSquareRoots = 0; int degree, degree2, numberOfSquareRoots = 0;
bool hasExtraSquareRoot = false; bool hasExtraSquareRoot = false;
/* FIXME
* For singular T, norm(I - T) >= 1 but maxNormForPade < 1, leads to infinite
* loop. We should move 0 eigenvalues to bottom right corner. We need not
* worry about tiny values (e.g. 1e-300) because they will reach 1 if
* repetitively sqrt'ed.
*
* [ T A ]^p [ T^p (T^-1 T^p A) ]
* [ ] = [ ]
* [ 0 0 ] [ 0 0 ]
*/
for (Index i=0; i < m_A.cols(); ++i)
eigen_assert(m_A(i,i) != RealScalar(0));
while (true) { while (true) {
IminusT = MatrixType::Identity(m_A.rows(), m_A.cols()) - T; IminusT = MatrixType::Identity(m_A.rows(), m_A.cols()) - T;
normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff(); normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff();

2
unsupported/test/matrix_power.cpp
View File

@ -181,6 +181,6 @@ void test_matrix_power()
CALL_SUBTEST_1(testMatrixVector(Matrix2f(), Vector2f(), 1e-4)); CALL_SUBTEST_1(testMatrixVector(Matrix2f(), Vector2f(), 1e-4));
CALL_SUBTEST_5(testMatrixVector(Matrix3cf(), Vector3cf(), 1e-4)); CALL_SUBTEST_5(testMatrixVector(Matrix3cf(), Vector3cf(), 1e-4));
CALL_SUBTEST_8(testMatrixVector(Matrix4f(), Vector4f(), 1e-4)); CALL_SUBTEST_8(testMatrixVector(Matrix4f(), Vector4f(), 1e-4));
CALL_SUBTEST_6(testMatrixVector(MatrixXf(8,8), VectorXf(8), 1e-3)); CALL_SUBTEST_6(testMatrixVector(MatrixXf(2,2), VectorXf(2), 1e-3)); // see bug 614
CALL_SUBTEST_9(testMatrixVector(MatrixXe(7,7), VectorXe(7), 1e-13)); CALL_SUBTEST_9(testMatrixVector(MatrixXe(7,7), VectorXe(7), 1e-13));
} }