25 #ifndef EIGEN_MATRIX_SQUARE_ROOT
26 #define EIGEN_MATRIX_SQUARE_ROOT
41 template <
typename MatrixType>
57 eigen_assert(A.rows() == A.cols());
68 template <
typename ResultType>
void compute(ResultType &result);
71 typedef typename MatrixType::Index Index;
72 typedef typename MatrixType::Scalar Scalar;
74 void computeDiagonalPartOfSqrt(MatrixType& sqrtT,
const MatrixType& T);
75 void computeOffDiagonalPartOfSqrt(MatrixType& sqrtT,
const MatrixType& T);
76 void compute2x2diagonalBlock(MatrixType& sqrtT,
const MatrixType& T,
typename MatrixType::Index i);
77 void compute1x1offDiagonalBlock(MatrixType& sqrtT,
const MatrixType& T,
78 typename MatrixType::Index i,
typename MatrixType::Index j);
79 void compute1x2offDiagonalBlock(MatrixType& sqrtT,
const MatrixType& T,
80 typename MatrixType::Index i,
typename MatrixType::Index j);
81 void compute2x1offDiagonalBlock(MatrixType& sqrtT,
const MatrixType& T,
82 typename MatrixType::Index i,
typename MatrixType::Index j);
83 void compute2x2offDiagonalBlock(MatrixType& sqrtT,
const MatrixType& T,
84 typename MatrixType::Index i,
typename MatrixType::Index j);
86 template <
typename SmallMatrixType>
87 static void solveAuxiliaryEquation(SmallMatrixType& X,
const SmallMatrixType& A,
88 const SmallMatrixType& B,
const SmallMatrixType& C);
90 const MatrixType& m_A;
93 template <
typename MatrixType>
94 template <
typename ResultType>
98 const RealSchur<MatrixType> schurOfA(m_A);
99 const MatrixType& T = schurOfA.matrixT();
100 const MatrixType& U = schurOfA.matrixU();
103 MatrixType sqrtT = MatrixType::Zero(m_A.rows(), m_A.rows());
104 computeDiagonalPartOfSqrt(sqrtT, T);
105 computeOffDiagonalPartOfSqrt(sqrtT, T);
108 result = U * sqrtT * U.adjoint();
113 template <
typename MatrixType>
117 const Index size = m_A.rows();
118 for (Index i = 0; i < size; i++) {
119 if (i == size - 1 || T.coeff(i+1, i) == 0) {
120 eigen_assert(T(i,i) > 0);
121 sqrtT.coeffRef(i,i) = internal::sqrt(T.coeff(i,i));
124 compute2x2diagonalBlock(sqrtT, T, i);
132 template <
typename MatrixType>
133 void MatrixSquareRootQuasiTriangular<MatrixType>::computeOffDiagonalPartOfSqrt(MatrixType& sqrtT,
136 const Index size = m_A.rows();
137 for (Index j = 1; j < size; j++) {
138 if (T.coeff(j, j-1) != 0)
140 for (Index i = j-1; i >= 0; i--) {
141 if (i > 0 && T.coeff(i, i-1) != 0)
143 bool iBlockIs2x2 = (i < size - 1) && (T.coeff(i+1, i) != 0);
144 bool jBlockIs2x2 = (j < size - 1) && (T.coeff(j+1, j) != 0);
145 if (iBlockIs2x2 && jBlockIs2x2)
146 compute2x2offDiagonalBlock(sqrtT, T, i, j);
147 else if (iBlockIs2x2 && !jBlockIs2x2)
148 compute2x1offDiagonalBlock(sqrtT, T, i, j);
149 else if (!iBlockIs2x2 && jBlockIs2x2)
150 compute1x2offDiagonalBlock(sqrtT, T, i, j);
151 else if (!iBlockIs2x2 && !jBlockIs2x2)
152 compute1x1offDiagonalBlock(sqrtT, T, i, j);
159 template <
typename MatrixType>
160 void MatrixSquareRootQuasiTriangular<MatrixType>
161 ::compute2x2diagonalBlock(MatrixType& sqrtT,
const MatrixType& T,
typename MatrixType::Index i)
165 Matrix<Scalar,2,2>
block = T.template block<2,2>(i,i);
166 EigenSolver<Matrix<Scalar,2,2> > es(block);
167 sqrtT.template block<2,2>(i,i)
168 = (es.eigenvectors() * es.eigenvalues().cwiseSqrt().asDiagonal() * es.eigenvectors().inverse()).
real();
174 template <
typename MatrixType>
175 void MatrixSquareRootQuasiTriangular<MatrixType>
176 ::compute1x1offDiagonalBlock(MatrixType& sqrtT,
const MatrixType& T,
177 typename MatrixType::Index i,
typename MatrixType::Index j)
179 Scalar tmp = (sqrtT.row(i).segment(i+1,j-i-1) * sqrtT.col(j).segment(i+1,j-i-1)).value();
180 sqrtT.coeffRef(i,j) = (T.coeff(i,j) - tmp) / (sqrtT.coeff(i,i) + sqrtT.coeff(j,j));
184 template <
typename MatrixType>
185 void MatrixSquareRootQuasiTriangular<MatrixType>
186 ::compute1x2offDiagonalBlock(MatrixType& sqrtT,
const MatrixType& T,
187 typename MatrixType::Index i,
typename MatrixType::Index j)
189 Matrix<Scalar,1,2> rhs = T.template block<1,2>(i,j);
191 rhs -= sqrtT.block(i, i+1, 1, j-i-1) * sqrtT.block(i+1, j, j-i-1, 2);
192 Matrix<Scalar,2,2> A = sqrtT.coeff(i,i) * Matrix<Scalar,2,2>::Identity();
193 A += sqrtT.template block<2,2>(j,j).transpose();
194 sqrtT.template block<1,2>(i,j).transpose() = A.fullPivLu().solve(rhs.transpose());
198 template <
typename MatrixType>
199 void MatrixSquareRootQuasiTriangular<MatrixType>
200 ::compute2x1offDiagonalBlock(MatrixType& sqrtT,
const MatrixType& T,
201 typename MatrixType::Index i,
typename MatrixType::Index j)
203 Matrix<Scalar,2,1> rhs = T.template block<2,1>(i,j);
205 rhs -= sqrtT.block(i, i+2, 2, j-i-2) * sqrtT.block(i+2, j, j-i-2, 1);
206 Matrix<Scalar,2,2> A = sqrtT.coeff(j,j) * Matrix<Scalar,2,2>::Identity();
207 A += sqrtT.template block<2,2>(i,i);
208 sqrtT.template block<2,1>(i,j) = A.fullPivLu().solve(rhs);
212 template <
typename MatrixType>
213 void MatrixSquareRootQuasiTriangular<MatrixType>
214 ::compute2x2offDiagonalBlock(MatrixType& sqrtT,
const MatrixType& T,
215 typename MatrixType::Index i,
typename MatrixType::Index j)
217 Matrix<Scalar,2,2> A = sqrtT.template block<2,2>(i,i);
218 Matrix<Scalar,2,2> B = sqrtT.template block<2,2>(j,j);
219 Matrix<Scalar,2,2> C = T.template block<2,2>(i,j);
221 C -= sqrtT.block(i, i+2, 2, j-i-2) * sqrtT.block(i+2, j, j-i-2, 2);
222 Matrix<Scalar,2,2> X;
223 solveAuxiliaryEquation(X, A, B, C);
224 sqrtT.template block<2,2>(i,j) = X;
228 template <
typename MatrixType>
229 template <
typename SmallMatrixType>
230 void MatrixSquareRootQuasiTriangular<MatrixType>
231 ::solveAuxiliaryEquation(SmallMatrixType& X,
const SmallMatrixType& A,
232 const SmallMatrixType& B,
const SmallMatrixType& C)
234 EIGEN_STATIC_ASSERT((internal::is_same<SmallMatrixType, Matrix<Scalar,2,2> >::value),
235 EIGEN_INTERNAL_ERROR_PLEASE_FILE_A_BUG_REPORT);
237 Matrix<Scalar,4,4> coeffMatrix = Matrix<Scalar,4,4>::Zero();
238 coeffMatrix.coeffRef(0,0) = A.coeff(0,0) + B.coeff(0,0);
239 coeffMatrix.coeffRef(1,1) = A.coeff(0,0) + B.coeff(1,1);
240 coeffMatrix.coeffRef(2,2) = A.coeff(1,1) + B.coeff(0,0);
241 coeffMatrix.coeffRef(3,3) = A.coeff(1,1) + B.coeff(1,1);
242 coeffMatrix.coeffRef(0,1) = B.coeff(1,0);
243 coeffMatrix.coeffRef(0,2) = A.coeff(0,1);
244 coeffMatrix.coeffRef(1,0) = B.coeff(0,1);
245 coeffMatrix.coeffRef(1,3) = A.coeff(0,1);
246 coeffMatrix.coeffRef(2,0) = A.coeff(1,0);
247 coeffMatrix.coeffRef(2,3) = B.coeff(1,0);
248 coeffMatrix.coeffRef(3,1) = A.coeff(1,0);
249 coeffMatrix.coeffRef(3,2) = B.coeff(0,1);
251 Matrix<Scalar,4,1> rhs;
252 rhs.coeffRef(0) = C.coeff(0,0);
253 rhs.coeffRef(1) = C.coeff(0,1);
254 rhs.coeffRef(2) = C.coeff(1,0);
255 rhs.coeffRef(3) = C.coeff(1,1);
257 Matrix<Scalar,4,1> result;
258 result = coeffMatrix.fullPivLu().solve(rhs);
260 X.coeffRef(0,0) = result.coeff(0);
261 X.coeffRef(0,1) = result.coeff(1);
262 X.coeffRef(1,0) = result.coeff(2);
263 X.coeffRef(1,1) = result.coeff(3);
278 template <
typename MatrixType>
285 eigen_assert(A.rows() == A.cols());
297 template <
typename ResultType>
void compute(ResultType &result);
300 const MatrixType& m_A;
303 template <
typename MatrixType>
304 template <
typename ResultType>
308 const ComplexSchur<MatrixType> schurOfA(m_A);
309 const MatrixType& T = schurOfA.matrixT();
310 const MatrixType& U = schurOfA.matrixU();
314 result.resize(m_A.rows(), m_A.cols());
315 typedef typename MatrixType::Index Index;
316 for (Index i = 0; i < m_A.rows(); i++) {
317 result.coeffRef(i,i) = internal::sqrt(T.coeff(i,i));
319 for (Index j = 1; j < m_A.cols(); j++) {
320 for (Index i = j-1; i >= 0; i--) {
321 typedef typename MatrixType::Scalar Scalar;
323 Scalar tmp = (result.row(i).segment(i+1,j-i-1) * result.col(j).segment(i+1,j-i-1)).value();
325 result.coeffRef(i,j) = (T.coeff(i,j) - tmp) / (result.coeff(i,i) + result.coeff(j,j));
331 tmp.noalias() = U * result.template triangularView<Upper>();
332 result.noalias() = tmp * U.adjoint();
343 template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
364 template <
typename ResultType>
void compute(ResultType &result);
370 template <
typename MatrixType>
378 eigen_assert(A.rows() == A.cols());
381 template <
typename ResultType>
void compute(ResultType &result)
384 const RealSchur<MatrixType> schurOfA(m_A);
385 const MatrixType& T = schurOfA.matrixT();
386 const MatrixType& U = schurOfA.matrixU();
389 MatrixSquareRootQuasiTriangular<MatrixType> tmp(T);
390 MatrixType sqrtT = MatrixType::Zero(m_A.rows(), m_A.rows());
394 result = U * sqrtT * U.adjoint();
398 const MatrixType& m_A;
404 template <
typename MatrixType>
405 class MatrixSquareRoot<MatrixType, 1>
412 eigen_assert(A.rows() == A.cols());
415 template <
typename ResultType>
void compute(ResultType &result)
418 const ComplexSchur<MatrixType> schurOfA(m_A);
419 const MatrixType& T = schurOfA.matrixT();
420 const MatrixType& U = schurOfA.matrixU();
423 MatrixSquareRootTriangular<MatrixType> tmp(T);
424 MatrixType sqrtT = MatrixType::Zero(m_A.rows(), m_A.rows());
428 result = U * sqrtT * U.adjoint();
432 const MatrixType& m_A;
449 :
public ReturnByValue<MatrixSquareRootReturnValue<Derived> >
451 typedef typename Derived::Index Index;
465 template <
typename ResultType>
466 inline void evalTo(ResultType& result)
const
468 const typename Derived::PlainObject srcEvaluated = m_src.eval();
473 Index rows()
const {
return m_src.rows(); }
474 Index cols()
const {
return m_src.cols(); }
477 const Derived& m_src;
483 template<
typename Derived>
484 struct traits<MatrixSquareRootReturnValue<Derived> >
486 typedef typename Derived::PlainObject ReturnType;
490 template <
typename Derived>
491 const MatrixSquareRootReturnValue<Derived> MatrixBase<Derived>::sqrt()
const
493 eigen_assert(rows() == cols());
494 return MatrixSquareRootReturnValue<Derived>(derived());
499 #endif // EIGEN_MATRIX_FUNCTION