mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-03-07 18:27:40 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix MatrixFunctions module.

Browse Source
This commit is contained in:
Gael Guennebaud 2016-06-03 09:21:35 +02:00
parent 82293f38d6
commit e8b922ca63
1 changed files with 16 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h
View File

@ -65,21 +65,6 @@ void matrix_sqrt_quasi_triangular_2x1_off_diagonal_block(const MatrixType& T, ty
sqrtT.template block<2,1>(i,j) = A.fullPivLu().solve(rhs);
}
// similar to compute1x1offDiagonalBlock()
template <typename MatrixType, typename ResultType>
void matrix_sqrt_quasi_triangular_2x2_off_diagonal_block(const MatrixType& T, typename MatrixType::Index i, typename MatrixType::Index j, ResultType& sqrtT)
{
typedef typename traits<MatrixType>::Scalar Scalar;
Matrix<Scalar,2,2> A = sqrtT.template block<2,2>(i,i);
Matrix<Scalar,2,2> B = sqrtT.template block<2,2>(j,j);
Matrix<Scalar,2,2> C = T.template block<2,2>(i,j);
if (j-i > 2)
C -= sqrtT.block(i, i+2, 2, j-i-2) * sqrtT.block(i+2, j, j-i-2, 2);
Matrix<Scalar,2,2> X;
matrix_sqrt_quasi_triangular_solve_auxiliary_equation(X, A, B, C);
sqrtT.template block<2,2>(i,j) = X;
}
// solves the equation A X + X B = C where all matrices are 2-by-2
template <typename MatrixType>
void matrix_sqrt_quasi_triangular_solve_auxiliary_equation(MatrixType& X, const MatrixType& A, const MatrixType& B, const MatrixType& C)
@ -98,13 +83,13 @@ void matrix_sqrt_quasi_triangular_solve_auxiliary_equation(MatrixType& X, const
coeffMatrix.coeffRef(2,3) = B.coeff(1,0);
coeffMatrix.coeffRef(3,1) = A.coeff(1,0);
coeffMatrix.coeffRef(3,2) = B.coeff(0,1);
Matrix<Scalar,4,1> rhs;
rhs.coeffRef(0) = C.coeff(0,0);
rhs.coeffRef(1) = C.coeff(0,1);
rhs.coeffRef(2) = C.coeff(1,0);
rhs.coeffRef(3) = C.coeff(1,1);
Matrix<Scalar,4,1> result;
result = coeffMatrix.fullPivLu().solve(rhs);
@ -114,6 +99,20 @@ void matrix_sqrt_quasi_triangular_solve_auxiliary_equation(MatrixType& X, const
X.coeffRef(1,1) = result.coeff(3);
}
// similar to compute1x1offDiagonalBlock()
template <typename MatrixType, typename ResultType>
void matrix_sqrt_quasi_triangular_2x2_off_diagonal_block(const MatrixType& T, typename MatrixType::Index i, typename MatrixType::Index j, ResultType& sqrtT)
{
typedef typename traits<MatrixType>::Scalar Scalar;
Matrix<Scalar,2,2> A = sqrtT.template block<2,2>(i,i);
Matrix<Scalar,2,2> B = sqrtT.template block<2,2>(j,j);
Matrix<Scalar,2,2> C = T.template block<2,2>(i,j);
if (j-i > 2)
C -= sqrtT.block(i, i+2, 2, j-i-2) * sqrtT.block(i+2, j, j-i-2, 2);
Matrix<Scalar,2,2> X;
matrix_sqrt_quasi_triangular_solve_auxiliary_equation(X, A, B, C);
sqrtT.template block<2,2>(i,j) = X;
}
// pre: T is quasi-upper-triangular and sqrtT is a zero matrix of the same size
// post: the diagonal blocks of sqrtT are the square roots of the diagonal blocks of T