10 #ifndef EIGEN_MATRIX_SQUARE_ROOT
11 #define EIGEN_MATRIX_SQUARE_ROOT
19 template <
typename MatrixType,
typename ResultType>
20 void matrix_sqrt_quasi_triangular_2x2_diagonal_block(
const MatrixType& T,
typename MatrixType::Index i, ResultType& sqrtT)
24 typedef typename traits<MatrixType>::Scalar Scalar;
25 Matrix<Scalar,2,2> block = T.template block<2,2>(i,i);
26 EigenSolver<Matrix<Scalar,2,2> > es(block);
27 sqrtT.template block<2,2>(i,i)
28 = (es.eigenvectors() * es.eigenvalues().cwiseSqrt().asDiagonal() * es.eigenvectors().inverse()).
real();
34 template <
typename MatrixType,
typename ResultType>
35 void matrix_sqrt_quasi_triangular_1x1_off_diagonal_block(
const MatrixType& T,
typename MatrixType::Index i,
typename MatrixType::Index j, ResultType& sqrtT)
37 typedef typename traits<MatrixType>::Scalar Scalar;
38 Scalar tmp = (sqrtT.row(i).segment(i+1,j-i-1) * sqrtT.col(j).segment(i+1,j-i-1)).value();
39 sqrtT.coeffRef(i,j) = (T.coeff(i,j) - tmp) / (sqrtT.coeff(i,i) + sqrtT.coeff(j,j));
43 template <
typename MatrixType,
typename ResultType>
44 void matrix_sqrt_quasi_triangular_1x2_off_diagonal_block(
const MatrixType& T,
typename MatrixType::Index i,
typename MatrixType::Index j, ResultType& sqrtT)
46 typedef typename traits<MatrixType>::Scalar Scalar;
47 Matrix<Scalar,1,2> rhs = T.template block<1,2>(i,j);
49 rhs -= sqrtT.block(i, i+1, 1, j-i-1) * sqrtT.block(i+1, j, j-i-1, 2);
51 A += sqrtT.template block<2,2>(j,j).transpose();
52 sqrtT.template block<1,2>(i,j).transpose() = A.fullPivLu().solve(rhs.transpose());
56 template <
typename MatrixType,
typename ResultType>
57 void matrix_sqrt_quasi_triangular_2x1_off_diagonal_block(
const MatrixType& T,
typename MatrixType::Index i,
typename MatrixType::Index j, ResultType& sqrtT)
59 typedef typename traits<MatrixType>::Scalar Scalar;
60 Matrix<Scalar,2,1> rhs = T.template block<2,1>(i,j);
62 rhs -= sqrtT.block(i, i+2, 2, j-i-2) * sqrtT.block(i+2, j, j-i-2, 1);
64 A += sqrtT.template block<2,2>(i,i);
65 sqrtT.template block<2,1>(i,j) = A.fullPivLu().solve(rhs);
69 template <
typename MatrixType>
70 void matrix_sqrt_quasi_triangular_solve_auxiliary_equation(MatrixType& X,
const MatrixType& A,
const MatrixType& B,
const MatrixType& C)
72 typedef typename traits<MatrixType>::Scalar Scalar;
73 Matrix<Scalar,4,4> coeffMatrix = Matrix<Scalar,4,4>::Zero();
74 coeffMatrix.coeffRef(0,0) = A.coeff(0,0) + B.coeff(0,0);
75 coeffMatrix.coeffRef(1,1) = A.coeff(0,0) + B.coeff(1,1);
76 coeffMatrix.coeffRef(2,2) = A.coeff(1,1) + B.coeff(0,0);
77 coeffMatrix.coeffRef(3,3) = A.coeff(1,1) + B.coeff(1,1);
78 coeffMatrix.coeffRef(0,1) = B.coeff(1,0);
79 coeffMatrix.coeffRef(0,2) = A.coeff(0,1);
80 coeffMatrix.coeffRef(1,0) = B.coeff(0,1);
81 coeffMatrix.coeffRef(1,3) = A.coeff(0,1);
82 coeffMatrix.coeffRef(2,0) = A.coeff(1,0);
83 coeffMatrix.coeffRef(2,3) = B.coeff(1,0);
84 coeffMatrix.coeffRef(3,1) = A.coeff(1,0);
85 coeffMatrix.coeffRef(3,2) = B.coeff(0,1);
87 Matrix<Scalar,4,1> rhs;
88 rhs.coeffRef(0) = C.coeff(0,0);
89 rhs.coeffRef(1) = C.coeff(0,1);
90 rhs.coeffRef(2) = C.coeff(1,0);
91 rhs.coeffRef(3) = C.coeff(1,1);
93 Matrix<Scalar,4,1> result;
94 result = coeffMatrix.fullPivLu().solve(rhs);
96 X.coeffRef(0,0) = result.coeff(0);
97 X.coeffRef(0,1) = result.coeff(1);
98 X.coeffRef(1,0) = result.coeff(2);
99 X.coeffRef(1,1) = result.coeff(3);
103 template <
typename MatrixType,
typename ResultType>
104 void matrix_sqrt_quasi_triangular_2x2_off_diagonal_block(
const MatrixType& T,
typename MatrixType::Index i,
typename MatrixType::Index j, ResultType& sqrtT)
106 typedef typename traits<MatrixType>::Scalar Scalar;
107 Matrix<Scalar,2,2> A = sqrtT.template block<2,2>(i,i);
108 Matrix<Scalar,2,2> B = sqrtT.template block<2,2>(j,j);
109 Matrix<Scalar,2,2> C = T.template block<2,2>(i,j);
111 C -= sqrtT.block(i, i+2, 2, j-i-2) * sqrtT.block(i+2, j, j-i-2, 2);
112 Matrix<Scalar,2,2> X;
113 matrix_sqrt_quasi_triangular_solve_auxiliary_equation(X, A, B, C);
114 sqrtT.template block<2,2>(i,j) = X;
119 template <
typename MatrixType,
typename ResultType>
120 void matrix_sqrt_quasi_triangular_diagonal(
const MatrixType& T, ResultType& sqrtT)
123 const Index size = T.rows();
124 for (
Index i = 0; i < size; i++) {
125 if (i == size - 1 || T.coeff(i+1, i) == 0) {
126 eigen_assert(T(i,i) >= 0);
127 sqrtT.coeffRef(i,i) =
sqrt(T.coeff(i,i));
130 matrix_sqrt_quasi_triangular_2x2_diagonal_block(T, i, sqrtT);
138 template <
typename MatrixType,
typename ResultType>
139 void matrix_sqrt_quasi_triangular_off_diagonal(
const MatrixType& T, ResultType& sqrtT)
141 const Index size = T.rows();
142 for (
Index j = 1; j < size; j++) {
143 if (T.coeff(j, j-1) != 0)
145 for (
Index i = j-1; i >= 0; i--) {
146 if (i > 0 && T.coeff(i, i-1) != 0)
148 bool iBlockIs2x2 = (i < size - 1) && (T.coeff(i+1, i) != 0);
149 bool jBlockIs2x2 = (j < size - 1) && (T.coeff(j+1, j) != 0);
150 if (iBlockIs2x2 && jBlockIs2x2)
151 matrix_sqrt_quasi_triangular_2x2_off_diagonal_block(T, i, j, sqrtT);
152 else if (iBlockIs2x2 && !jBlockIs2x2)
153 matrix_sqrt_quasi_triangular_2x1_off_diagonal_block(T, i, j, sqrtT);
154 else if (!iBlockIs2x2 && jBlockIs2x2)
155 matrix_sqrt_quasi_triangular_1x2_off_diagonal_block(T, i, j, sqrtT);
156 else if (!iBlockIs2x2 && !jBlockIs2x2)
157 matrix_sqrt_quasi_triangular_1x1_off_diagonal_block(T, i, j, sqrtT);
179 template <
typename MatrixType,
typename ResultType>
182 eigen_assert(arg.rows() == arg.cols());
183 result.resize(arg.rows(), arg.cols());
184 internal::matrix_sqrt_quasi_triangular_diagonal(arg, result);
185 internal::matrix_sqrt_quasi_triangular_off_diagonal(arg, result);
203 template <
typename MatrixType,
typename ResultType>
207 typedef typename MatrixType::Scalar Scalar;
209 eigen_assert(arg.rows() == arg.cols());
213 result.resize(arg.rows(), arg.cols());
214 for (
Index i = 0; i < arg.rows(); i++) {
215 result.coeffRef(i,i) =
sqrt(arg.coeff(i,i));
217 for (
Index j = 1; j < arg.cols(); j++) {
218 for (
Index i = j-1; i >= 0; i--) {
220 Scalar tmp = (result.row(i).segment(i+1,j-i-1) * result.col(j).segment(i+1,j-i-1)).value();
222 result.coeffRef(i,j) = (arg.coeff(i,j) - tmp) / (result.coeff(i,i) + result.coeff(j,j));
237 template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
238 struct matrix_sqrt_compute
247 template <
typename ResultType>
static void run(
const MatrixType &
arg, ResultType &result);
253 template <
typename MatrixType>
254 struct matrix_sqrt_compute<MatrixType, 0>
256 template <
typename ResultType>
257 static void run(
const MatrixType &
arg, ResultType &result)
259 eigen_assert(arg.rows() == arg.cols());
262 const RealSchur<MatrixType> schurOfA(arg);
263 const MatrixType& T = schurOfA.matrixT();
264 const MatrixType& U = schurOfA.matrixU();
267 MatrixType sqrtT = MatrixType::Zero(arg.rows(), arg.cols());
271 result = U * sqrtT * U.adjoint();
278 template <
typename MatrixType>
279 struct matrix_sqrt_compute<MatrixType, 1>
281 template <
typename ResultType>
282 static void run(
const MatrixType &arg, ResultType &result)
284 eigen_assert(arg.rows() == arg.cols());
287 const ComplexSchur<MatrixType> schurOfA(arg);
288 const MatrixType& T = schurOfA.matrixT();
289 const MatrixType& U = schurOfA.matrixU();
296 result = U * (sqrtT.template triangularView<Upper>() * U.adjoint());
315 :
public ReturnByValue<MatrixSquareRootReturnValue<Derived> >
318 typedef typename internal::ref_selector<Derived>::type DerivedNested;
333 template <
typename ResultType>
334 inline void evalTo(ResultType& result)
const
336 typedef typename internal::nested_eval<Derived, 10>::type DerivedEvalType;
337 typedef typename internal::remove_all<DerivedEvalType>::type DerivedEvalTypeClean;
338 DerivedEvalType tmp(m_src);
339 internal::matrix_sqrt_compute<DerivedEvalTypeClean>::run(tmp, result);
342 Index rows()
const {
return m_src.rows(); }
343 Index cols()
const {
return m_src.cols(); }
346 const DerivedNested m_src;
350 template<
typename Derived>
351 struct traits<MatrixSquareRootReturnValue<Derived> >
353 typedef typename Derived::PlainObject ReturnType;
357 template <
typename Derived>
358 const MatrixSquareRootReturnValue<Derived> MatrixBase<Derived>::sqrt()
const
360 eigen_assert(rows() == cols());
361 return MatrixSquareRootReturnValue<Derived>(derived());
366 #endif // EIGEN_MATRIX_FUNCTION
Proxy for the matrix square root of some matrix (expression).
Definition: MatrixSquareRoot.h:314
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_arg_op< typename Derived::Scalar >, const Derived > arg(const Eigen::ArrayBase< Derived > &x)
void matrix_sqrt_quasi_triangular(const MatrixType &arg, ResultType &result)
Compute matrix square root of quasi-triangular matrix.
Definition: MatrixSquareRoot.h:180
void evalTo(ResultType &result) const
Compute the matrix square root.
Definition: MatrixSquareRoot.h:334
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
void matrix_sqrt_triangular(const MatrixType &arg, ResultType &result)
Compute matrix square root of triangular matrix.
Definition: MatrixSquareRoot.h:204
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_real_op< typename Derived::Scalar >, const Derived > real(const Eigen::ArrayBase< Derived > &x)
static const IdentityReturnType Identity()
MatrixSquareRootReturnValue(const Derived &src)
Constructor.
Definition: MatrixSquareRoot.h:326