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);
50 Matrix<Scalar,2,2> A = sqrtT.coeff(i,i) * Matrix<Scalar,2,2>::Identity();
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);
63 Matrix<Scalar,2,2> A = sqrtT.coeff(j,j) * Matrix<Scalar,2,2>::Identity();
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,
typename ResultType>
70 void matrix_sqrt_quasi_triangular_2x2_off_diagonal_block(
const MatrixType& T,
typename MatrixType::Index i,
typename MatrixType::Index j, ResultType& sqrtT)
72 typedef typename traits<MatrixType>::Scalar Scalar;
73 Matrix<Scalar,2,2> A = sqrtT.template block<2,2>(i,i);
74 Matrix<Scalar,2,2> B = sqrtT.template block<2,2>(j,j);
75 Matrix<Scalar,2,2> C = T.template block<2,2>(i,j);
77 C -= sqrtT.block(i, i+2, 2, j-i-2) * sqrtT.block(i+2, j, j-i-2, 2);
79 matrix_sqrt_quasi_triangular_solve_auxiliary_equation(X, A, B, C);
80 sqrtT.template block<2,2>(i,j) = X;
84 template <
typename MatrixType>
85 void matrix_sqrt_quasi_triangular_solve_auxiliary_equation(MatrixType& X,
const MatrixType& A,
const MatrixType& B,
const MatrixType& C)
87 typedef typename traits<MatrixType>::Scalar Scalar;
88 Matrix<Scalar,4,4> coeffMatrix = Matrix<Scalar,4,4>::Zero();
89 coeffMatrix.coeffRef(0,0) = A.coeff(0,0) + B.coeff(0,0);
90 coeffMatrix.coeffRef(1,1) = A.coeff(0,0) + B.coeff(1,1);
91 coeffMatrix.coeffRef(2,2) = A.coeff(1,1) + B.coeff(0,0);
92 coeffMatrix.coeffRef(3,3) = A.coeff(1,1) + B.coeff(1,1);
93 coeffMatrix.coeffRef(0,1) = B.coeff(1,0);
94 coeffMatrix.coeffRef(0,2) = A.coeff(0,1);
95 coeffMatrix.coeffRef(1,0) = B.coeff(0,1);
96 coeffMatrix.coeffRef(1,3) = A.coeff(0,1);
97 coeffMatrix.coeffRef(2,0) = A.coeff(1,0);
98 coeffMatrix.coeffRef(2,3) = B.coeff(1,0);
99 coeffMatrix.coeffRef(3,1) = A.coeff(1,0);
100 coeffMatrix.coeffRef(3,2) = B.coeff(0,1);
102 Matrix<Scalar,4,1> rhs;
103 rhs.coeffRef(0) = C.coeff(0,0);
104 rhs.coeffRef(1) = C.coeff(0,1);
105 rhs.coeffRef(2) = C.coeff(1,0);
106 rhs.coeffRef(3) = C.coeff(1,1);
108 Matrix<Scalar,4,1> result;
109 result = coeffMatrix.fullPivLu().solve(rhs);
111 X.coeffRef(0,0) = result.coeff(0);
112 X.coeffRef(0,1) = result.coeff(1);
113 X.coeffRef(1,0) = result.coeff(2);
114 X.coeffRef(1,1) = result.coeff(3);
120 template <
typename MatrixType,
typename ResultType>
121 void matrix_sqrt_quasi_triangular_diagonal(
const MatrixType& T, ResultType& sqrtT)
124 typedef typename MatrixType::Index Index;
125 const Index size = T.rows();
126 for (Index i = 0; i < size; i++) {
127 if (i == size - 1 || T.coeff(i+1, i) == 0) {
128 eigen_assert(T(i,i) >= 0);
129 sqrtT.coeffRef(i,i) = sqrt(T.coeff(i,i));
132 matrix_sqrt_quasi_triangular_2x2_diagonal_block(T, i, sqrtT);
140 template <
typename MatrixType,
typename ResultType>
141 void matrix_sqrt_quasi_triangular_off_diagonal(
const MatrixType& T, ResultType& sqrtT)
143 typedef typename MatrixType::Index Index;
144 const Index size = T.rows();
145 for (Index j = 1; j < size; j++) {
146 if (T.coeff(j, j-1) != 0)
148 for (Index i = j-1; i >= 0; i--) {
149 if (i > 0 && T.coeff(i, i-1) != 0)
151 bool iBlockIs2x2 = (i < size - 1) && (T.coeff(i+1, i) != 0);
152 bool jBlockIs2x2 = (j < size - 1) && (T.coeff(j+1, j) != 0);
153 if (iBlockIs2x2 && jBlockIs2x2)
154 matrix_sqrt_quasi_triangular_2x2_off_diagonal_block(T, i, j, sqrtT);
155 else if (iBlockIs2x2 && !jBlockIs2x2)
156 matrix_sqrt_quasi_triangular_2x1_off_diagonal_block(T, i, j, sqrtT);
157 else if (!iBlockIs2x2 && jBlockIs2x2)
158 matrix_sqrt_quasi_triangular_1x2_off_diagonal_block(T, i, j, sqrtT);
159 else if (!iBlockIs2x2 && !jBlockIs2x2)
160 matrix_sqrt_quasi_triangular_1x1_off_diagonal_block(T, i, j, sqrtT);
182 template <
typename MatrixType,
typename ResultType>
185 eigen_assert(arg.rows() == arg.cols());
186 result.resize(arg.rows(), arg.cols());
187 internal::matrix_sqrt_quasi_triangular_diagonal(arg, result);
188 internal::matrix_sqrt_quasi_triangular_off_diagonal(arg, result);
206 template <
typename MatrixType,
typename ResultType>
210 typedef typename MatrixType::Index Index;
211 typedef typename MatrixType::Scalar Scalar;
213 eigen_assert(arg.rows() == arg.cols());
217 result.resize(arg.rows(), arg.cols());
218 for (Index i = 0; i < arg.rows(); i++) {
219 result.coeffRef(i,i) = sqrt(arg.coeff(i,i));
221 for (Index j = 1; j < arg.cols(); j++) {
222 for (Index i = j-1; i >= 0; i--) {
224 Scalar tmp = (result.row(i).segment(i+1,j-i-1) * result.col(j).segment(i+1,j-i-1)).value();
226 result.coeffRef(i,j) = (arg.coeff(i,j) - tmp) / (result.coeff(i,i) + result.coeff(j,j));
241 template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
242 struct matrix_sqrt_compute
251 template <
typename ResultType>
static void run(
const MatrixType &arg, ResultType &result);
257 template <
typename MatrixType>
258 struct matrix_sqrt_compute<MatrixType, 0>
260 template <
typename ResultType>
261 static void run(
const MatrixType &arg, ResultType &result)
263 eigen_assert(arg.rows() == arg.cols());
266 const RealSchur<MatrixType> schurOfA(arg);
267 const MatrixType& T = schurOfA.matrixT();
268 const MatrixType& U = schurOfA.matrixU();
271 MatrixType sqrtT = MatrixType::Zero(arg.rows(), arg.cols());
275 result = U * sqrtT * U.adjoint();
282 template <
typename MatrixType>
283 struct matrix_sqrt_compute<MatrixType, 1>
285 template <
typename ResultType>
286 static void run(
const MatrixType &arg, ResultType &result)
288 eigen_assert(arg.rows() == arg.cols());
291 const ComplexSchur<MatrixType> schurOfA(arg);
292 const MatrixType& T = schurOfA.matrixT();
293 const MatrixType& U = schurOfA.matrixU();
300 result = U * (sqrtT.template triangularView<Upper>() * U.adjoint());
319 :
public ReturnByValue<MatrixSquareRootReturnValue<Derived> >
322 typedef typename Derived::Index Index;
323 typedef typename internal::ref_selector<Derived>::type DerivedNested;
338 template <
typename ResultType>
339 inline void evalTo(ResultType& result)
const 341 typedef typename internal::nested_eval<Derived, 10>::type DerivedEvalType;
342 typedef typename internal::remove_all<DerivedEvalType>::type DerivedEvalTypeClean;
343 DerivedEvalType tmp(m_src);
344 internal::matrix_sqrt_compute<DerivedEvalTypeClean>::run(tmp, result);
347 Index rows()
const {
return m_src.rows(); }
348 Index cols()
const {
return m_src.cols(); }
351 const DerivedNested m_src;
355 template<
typename Derived>
356 struct traits<MatrixSquareRootReturnValue<Derived> >
358 typedef typename Derived::PlainObject ReturnType;
362 template <
typename Derived>
365 eigen_assert(rows() == cols());
371 #endif // EIGEN_MATRIX_FUNCTION Proxy for the matrix square root of some matrix (expression).
Definition: MatrixSquareRoot.h:318
Namespace containing all symbols from the Eigen library.
Definition: CXX11Meta.h:13
void evalTo(ResultType &result) const
Compute the matrix square root.
Definition: MatrixSquareRoot.h:339
MatrixSquareRootReturnValue(const Derived &src)
Constructor.
Definition: MatrixSquareRoot.h:331
void matrix_sqrt_quasi_triangular(const MatrixType &arg, ResultType &result)
Compute matrix square root of quasi-triangular matrix.
Definition: MatrixSquareRoot.h:183
void matrix_sqrt_triangular(const MatrixType &arg, ResultType &result)
Compute matrix square root of triangular matrix.
Definition: MatrixSquareRoot.h:207