diff options
Diffstat (limited to 'unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h')
| -rw-r--r-- | unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h | 33 |
1 files changed, 16 insertions, 17 deletions
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h b/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h index 9f08c6162..afd88ec4d 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h @@ -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 |