| 62 | |
| 63 | template <typename MatrixType> |
| 64 | MatrixType MatrixFunctionAtomic<MatrixType>::compute(const MatrixType& A) |
| 65 | { |
| 66 | // TODO: Use that A is upper triangular |
| 67 | typedef typename NumTraits<Scalar>::Real RealScalar; |
| 68 | typedef typename MatrixType::Index Index; |
| 69 | Index rows = A.rows(); |
| 70 | Scalar avgEival = A.trace() / Scalar(RealScalar(rows)); |
| 71 | MatrixType Ashifted = A - avgEival * MatrixType::Identity(rows, rows); |
| 72 | RealScalar mu = matrix_function_compute_mu(Ashifted); |
| 73 | MatrixType F = m_f(avgEival, 0) * MatrixType::Identity(rows, rows); |
| 74 | MatrixType P = Ashifted; |
| 75 | MatrixType Fincr; |
| 76 | for (Index s = 1; s < 1.1 * rows + 10; s++) { // upper limit is fairly arbitrary |
| 77 | Fincr = m_f(avgEival, static_cast<int>(s)) * P; |
| 78 | F += Fincr; |
| 79 | P = Scalar(RealScalar(1.0/(s + 1))) * P * Ashifted; |
| 80 | |
| 81 | // test whether Taylor series converged |
| 82 | const RealScalar F_norm = F.cwiseAbs().rowwise().sum().maxCoeff(); |
| 83 | const RealScalar Fincr_norm = Fincr.cwiseAbs().rowwise().sum().maxCoeff(); |
| 84 | if (Fincr_norm < NumTraits<Scalar>::epsilon() * F_norm) { |
| 85 | RealScalar delta = 0; |
| 86 | RealScalar rfactorial = 1; |
| 87 | for (Index r = 0; r < rows; r++) { |
| 88 | RealScalar mx = 0; |
| 89 | for (Index i = 0; i < rows; i++) |
| 90 | mx = (std::max)(mx, std::abs(m_f(Ashifted(i, i) + avgEival, static_cast<int>(s+r)))); |
| 91 | if (r != 0) |
| 92 | rfactorial *= RealScalar(r); |
| 93 | delta = (std::max)(delta, mx / rfactorial); |
| 94 | } |
| 95 | const RealScalar P_norm = P.cwiseAbs().rowwise().sum().maxCoeff(); |
| 96 | if (mu * delta * P_norm < NumTraits<Scalar>::epsilon() * F_norm) // series converged |
| 97 | break; |
| 98 | } |
| 99 | } |
| 100 | return F; |
| 101 | } |
| 102 | |
| 103 | /** \brief Find cluster in \p clusters containing some value |
| 104 | * \param[in] key Value to find |
no test coverage detected