Attachment 352 Details for Bug 614 – Proposed patch, including regression test

This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen

[patch] Proposed patch, including regression test

MatrixPower.patch (text/plain), 2.27 KB, created by Chen-Pang He on 2013-06-21 06:52:26 UTC

(hide)

Description:

Filename:

MIME Type:

Creator: Chen-Pang He

Created: 2013-06-21 06:52:26 UTC

Size: 2.27 KB

patch

obsolete

>diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
>--- a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
>+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
>@@ -294,6 +294,7 @@
> template<typename MatrixType>
> void MatrixPowerTriangularAtomic<MatrixType>::compute(MatrixType& res, RealScalar p) const
> {
>+  res.resizeLike(m_A);
>   switch (m_A.rows()) {
>     case 0:
>       break;
>@@ -320,6 +321,7 @@
>   res += MatrixType::Identity(IminusT.rows(), IminusT.cols());
> }
> 
>+// this function assumes that res has the correct size (see bug 614)
> template<typename MatrixType>
> void MatrixPowerTriangularAtomic<MatrixType>::compute2x2(MatrixType& res, RealScalar p) const
> {
>@@ -357,10 +359,11 @@
> 						   9.134603732914548552537150753385375e-2L; // quadruple precision
>   MatrixType IminusT, sqrtT, T = m_A.template triangularView<Upper>();
>   RealScalar normIminusT;
>-  int degree, degree2, numberOfSquareRoots = 0;
>+  int degree, degree2, numberOfSquareRoots;
>   bool hasExtraSquareRoot = false;
> 
>-  while (true) {
>+  // Sqrt at most 64 times to avoid infinite loop for ill conditioned (I - T)
>+  for (numberOfSquareRoots = 0; numberOfSquareRoots < 64; ++numberOfSquareRoots) {
>     IminusT = MatrixType::Identity(m_A.rows(), m_A.cols()) - T;
>     normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff();
>     if (normIminusT < maxNormForPade) {
>@@ -372,7 +375,6 @@
>     }
>     MatrixSquareRootTriangular<MatrixType>(T).compute(sqrtT);
>     T = sqrtT.template triangularView<Upper>();
>-    ++numberOfSquareRoots;
>   }
>   computePade(degree, IminusT, res, p);
> 
>diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp
>--- a/unsupported/test/matrix_power.cpp
>+++ b/unsupported/test/matrix_power.cpp
>@@ -181,6 +181,6 @@
>   CALL_SUBTEST_1(testMatrixVector(Matrix2f(),         Vector2f(),    1e-4));
>   CALL_SUBTEST_5(testMatrixVector(Matrix3cf(),        Vector3cf(),   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));
> }
>

Actions: View | Diff

Attachments on bug 614: 347 | 348 | 352 | 355