MatrixFunction.h (22671B)
1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2009-2011, 2013 Jitse Niesen <jitse@maths.leeds.ac.uk> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #ifndef EIGEN_MATRIX_FUNCTION_H 11 #define EIGEN_MATRIX_FUNCTION_H 12 13 #include "StemFunction.h" 14 15 16 namespace Eigen { 17 18 namespace internal { 19 20 /** \brief Maximum distance allowed between eigenvalues to be considered "close". */ 21 static const float matrix_function_separation = 0.1f; 22 23 /** \ingroup MatrixFunctions_Module 24 * \class MatrixFunctionAtomic 25 * \brief Helper class for computing matrix functions of atomic matrices. 26 * 27 * Here, an atomic matrix is a triangular matrix whose diagonal entries are close to each other. 28 */ 29 template <typename MatrixType> 30 class MatrixFunctionAtomic 31 { 32 public: 33 34 typedef typename MatrixType::Scalar Scalar; 35 typedef typename stem_function<Scalar>::type StemFunction; 36 37 /** \brief Constructor 38 * \param[in] f matrix function to compute. 39 */ 40 MatrixFunctionAtomic(StemFunction f) : m_f(f) { } 41 42 /** \brief Compute matrix function of atomic matrix 43 * \param[in] A argument of matrix function, should be upper triangular and atomic 44 * \returns f(A), the matrix function evaluated at the given matrix 45 */ 46 MatrixType compute(const MatrixType& A); 47 48 private: 49 StemFunction* m_f; 50 }; 51 52 template <typename MatrixType> 53 typename NumTraits<typename MatrixType::Scalar>::Real matrix_function_compute_mu(const MatrixType& A) 54 { 55 typedef typename plain_col_type<MatrixType>::type VectorType; 56 Index rows = A.rows(); 57 const MatrixType N = MatrixType::Identity(rows, rows) - A; 58 VectorType e = VectorType::Ones(rows); 59 N.template triangularView<Upper>().solveInPlace(e); 60 return e.cwiseAbs().maxCoeff(); 61 } 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 Index rows = A.rows(); 69 Scalar avgEival = A.trace() / Scalar(RealScalar(rows)); 70 MatrixType Ashifted = A - avgEival * MatrixType::Identity(rows, rows); 71 RealScalar mu = matrix_function_compute_mu(Ashifted); 72 MatrixType F = m_f(avgEival, 0) * MatrixType::Identity(rows, rows); 73 MatrixType P = Ashifted; 74 MatrixType Fincr; 75 for (Index s = 1; double(s) < 1.1 * double(rows) + 10.0; s++) { // upper limit is fairly arbitrary 76 Fincr = m_f(avgEival, static_cast<int>(s)) * P; 77 F += Fincr; 78 P = Scalar(RealScalar(1)/RealScalar(s + 1)) * P * Ashifted; 79 80 // test whether Taylor series converged 81 const RealScalar F_norm = F.cwiseAbs().rowwise().sum().maxCoeff(); 82 const RealScalar Fincr_norm = Fincr.cwiseAbs().rowwise().sum().maxCoeff(); 83 if (Fincr_norm < NumTraits<Scalar>::epsilon() * F_norm) { 84 RealScalar delta = 0; 85 RealScalar rfactorial = 1; 86 for (Index r = 0; r < rows; r++) { 87 RealScalar mx = 0; 88 for (Index i = 0; i < rows; i++) 89 mx = (std::max)(mx, std::abs(m_f(Ashifted(i, i) + avgEival, static_cast<int>(s+r)))); 90 if (r != 0) 91 rfactorial *= RealScalar(r); 92 delta = (std::max)(delta, mx / rfactorial); 93 } 94 const RealScalar P_norm = P.cwiseAbs().rowwise().sum().maxCoeff(); 95 if (mu * delta * P_norm < NumTraits<Scalar>::epsilon() * F_norm) // series converged 96 break; 97 } 98 } 99 return F; 100 } 101 102 /** \brief Find cluster in \p clusters containing some value 103 * \param[in] key Value to find 104 * \returns Iterator to cluster containing \p key, or \c clusters.end() if no cluster in \p m_clusters 105 * contains \p key. 106 */ 107 template <typename Index, typename ListOfClusters> 108 typename ListOfClusters::iterator matrix_function_find_cluster(Index key, ListOfClusters& clusters) 109 { 110 typename std::list<Index>::iterator j; 111 for (typename ListOfClusters::iterator i = clusters.begin(); i != clusters.end(); ++i) { 112 j = std::find(i->begin(), i->end(), key); 113 if (j != i->end()) 114 return i; 115 } 116 return clusters.end(); 117 } 118 119 /** \brief Partition eigenvalues in clusters of ei'vals close to each other 120 * 121 * \param[in] eivals Eigenvalues 122 * \param[out] clusters Resulting partition of eigenvalues 123 * 124 * The partition satisfies the following two properties: 125 * # Any eigenvalue in a certain cluster is at most matrix_function_separation() away from another eigenvalue 126 * in the same cluster. 127 * # The distance between two eigenvalues in different clusters is more than matrix_function_separation(). 128 * The implementation follows Algorithm 4.1 in the paper of Davies and Higham. 129 */ 130 template <typename EivalsType, typename Cluster> 131 void matrix_function_partition_eigenvalues(const EivalsType& eivals, std::list<Cluster>& clusters) 132 { 133 typedef typename EivalsType::RealScalar RealScalar; 134 for (Index i=0; i<eivals.rows(); ++i) { 135 // Find cluster containing i-th ei'val, adding a new cluster if necessary 136 typename std::list<Cluster>::iterator qi = matrix_function_find_cluster(i, clusters); 137 if (qi == clusters.end()) { 138 Cluster l; 139 l.push_back(i); 140 clusters.push_back(l); 141 qi = clusters.end(); 142 --qi; 143 } 144 145 // Look for other element to add to the set 146 for (Index j=i+1; j<eivals.rows(); ++j) { 147 if (abs(eivals(j) - eivals(i)) <= RealScalar(matrix_function_separation) 148 && std::find(qi->begin(), qi->end(), j) == qi->end()) { 149 typename std::list<Cluster>::iterator qj = matrix_function_find_cluster(j, clusters); 150 if (qj == clusters.end()) { 151 qi->push_back(j); 152 } else { 153 qi->insert(qi->end(), qj->begin(), qj->end()); 154 clusters.erase(qj); 155 } 156 } 157 } 158 } 159 } 160 161 /** \brief Compute size of each cluster given a partitioning */ 162 template <typename ListOfClusters, typename Index> 163 void matrix_function_compute_cluster_size(const ListOfClusters& clusters, Matrix<Index, Dynamic, 1>& clusterSize) 164 { 165 const Index numClusters = static_cast<Index>(clusters.size()); 166 clusterSize.setZero(numClusters); 167 Index clusterIndex = 0; 168 for (typename ListOfClusters::const_iterator cluster = clusters.begin(); cluster != clusters.end(); ++cluster) { 169 clusterSize[clusterIndex] = cluster->size(); 170 ++clusterIndex; 171 } 172 } 173 174 /** \brief Compute start of each block using clusterSize */ 175 template <typename VectorType> 176 void matrix_function_compute_block_start(const VectorType& clusterSize, VectorType& blockStart) 177 { 178 blockStart.resize(clusterSize.rows()); 179 blockStart(0) = 0; 180 for (Index i = 1; i < clusterSize.rows(); i++) { 181 blockStart(i) = blockStart(i-1) + clusterSize(i-1); 182 } 183 } 184 185 /** \brief Compute mapping of eigenvalue indices to cluster indices */ 186 template <typename EivalsType, typename ListOfClusters, typename VectorType> 187 void matrix_function_compute_map(const EivalsType& eivals, const ListOfClusters& clusters, VectorType& eivalToCluster) 188 { 189 eivalToCluster.resize(eivals.rows()); 190 Index clusterIndex = 0; 191 for (typename ListOfClusters::const_iterator cluster = clusters.begin(); cluster != clusters.end(); ++cluster) { 192 for (Index i = 0; i < eivals.rows(); ++i) { 193 if (std::find(cluster->begin(), cluster->end(), i) != cluster->end()) { 194 eivalToCluster[i] = clusterIndex; 195 } 196 } 197 ++clusterIndex; 198 } 199 } 200 201 /** \brief Compute permutation which groups ei'vals in same cluster together */ 202 template <typename DynVectorType, typename VectorType> 203 void matrix_function_compute_permutation(const DynVectorType& blockStart, const DynVectorType& eivalToCluster, VectorType& permutation) 204 { 205 DynVectorType indexNextEntry = blockStart; 206 permutation.resize(eivalToCluster.rows()); 207 for (Index i = 0; i < eivalToCluster.rows(); i++) { 208 Index cluster = eivalToCluster[i]; 209 permutation[i] = indexNextEntry[cluster]; 210 ++indexNextEntry[cluster]; 211 } 212 } 213 214 /** \brief Permute Schur decomposition in U and T according to permutation */ 215 template <typename VectorType, typename MatrixType> 216 void matrix_function_permute_schur(VectorType& permutation, MatrixType& U, MatrixType& T) 217 { 218 for (Index i = 0; i < permutation.rows() - 1; i++) { 219 Index j; 220 for (j = i; j < permutation.rows(); j++) { 221 if (permutation(j) == i) break; 222 } 223 eigen_assert(permutation(j) == i); 224 for (Index k = j-1; k >= i; k--) { 225 JacobiRotation<typename MatrixType::Scalar> rotation; 226 rotation.makeGivens(T(k, k+1), T(k+1, k+1) - T(k, k)); 227 T.applyOnTheLeft(k, k+1, rotation.adjoint()); 228 T.applyOnTheRight(k, k+1, rotation); 229 U.applyOnTheRight(k, k+1, rotation); 230 std::swap(permutation.coeffRef(k), permutation.coeffRef(k+1)); 231 } 232 } 233 } 234 235 /** \brief Compute block diagonal part of matrix function. 236 * 237 * This routine computes the matrix function applied to the block diagonal part of \p T (which should be 238 * upper triangular), with the blocking given by \p blockStart and \p clusterSize. The matrix function of 239 * each diagonal block is computed by \p atomic. The off-diagonal parts of \p fT are set to zero. 240 */ 241 template <typename MatrixType, typename AtomicType, typename VectorType> 242 void matrix_function_compute_block_atomic(const MatrixType& T, AtomicType& atomic, const VectorType& blockStart, const VectorType& clusterSize, MatrixType& fT) 243 { 244 fT.setZero(T.rows(), T.cols()); 245 for (Index i = 0; i < clusterSize.rows(); ++i) { 246 fT.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i)) 247 = atomic.compute(T.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i))); 248 } 249 } 250 251 /** \brief Solve a triangular Sylvester equation AX + XB = C 252 * 253 * \param[in] A the matrix A; should be square and upper triangular 254 * \param[in] B the matrix B; should be square and upper triangular 255 * \param[in] C the matrix C; should have correct size. 256 * 257 * \returns the solution X. 258 * 259 * If A is m-by-m and B is n-by-n, then both C and X are m-by-n. The (i,j)-th component of the Sylvester 260 * equation is 261 * \f[ 262 * \sum_{k=i}^m A_{ik} X_{kj} + \sum_{k=1}^j X_{ik} B_{kj} = C_{ij}. 263 * \f] 264 * This can be re-arranged to yield: 265 * \f[ 266 * X_{ij} = \frac{1}{A_{ii} + B_{jj}} \Bigl( C_{ij} 267 * - \sum_{k=i+1}^m A_{ik} X_{kj} - \sum_{k=1}^{j-1} X_{ik} B_{kj} \Bigr). 268 * \f] 269 * It is assumed that A and B are such that the numerator is never zero (otherwise the Sylvester equation 270 * does not have a unique solution). In that case, these equations can be evaluated in the order 271 * \f$ i=m,\ldots,1 \f$ and \f$ j=1,\ldots,n \f$. 272 */ 273 template <typename MatrixType> 274 MatrixType matrix_function_solve_triangular_sylvester(const MatrixType& A, const MatrixType& B, const MatrixType& C) 275 { 276 eigen_assert(A.rows() == A.cols()); 277 eigen_assert(A.isUpperTriangular()); 278 eigen_assert(B.rows() == B.cols()); 279 eigen_assert(B.isUpperTriangular()); 280 eigen_assert(C.rows() == A.rows()); 281 eigen_assert(C.cols() == B.rows()); 282 283 typedef typename MatrixType::Scalar Scalar; 284 285 Index m = A.rows(); 286 Index n = B.rows(); 287 MatrixType X(m, n); 288 289 for (Index i = m - 1; i >= 0; --i) { 290 for (Index j = 0; j < n; ++j) { 291 292 // Compute AX = \sum_{k=i+1}^m A_{ik} X_{kj} 293 Scalar AX; 294 if (i == m - 1) { 295 AX = 0; 296 } else { 297 Matrix<Scalar,1,1> AXmatrix = A.row(i).tail(m-1-i) * X.col(j).tail(m-1-i); 298 AX = AXmatrix(0,0); 299 } 300 301 // Compute XB = \sum_{k=1}^{j-1} X_{ik} B_{kj} 302 Scalar XB; 303 if (j == 0) { 304 XB = 0; 305 } else { 306 Matrix<Scalar,1,1> XBmatrix = X.row(i).head(j) * B.col(j).head(j); 307 XB = XBmatrix(0,0); 308 } 309 310 X(i,j) = (C(i,j) - AX - XB) / (A(i,i) + B(j,j)); 311 } 312 } 313 return X; 314 } 315 316 /** \brief Compute part of matrix function above block diagonal. 317 * 318 * This routine completes the computation of \p fT, denoting a matrix function applied to the triangular 319 * matrix \p T. It assumes that the block diagonal part of \p fT has already been computed. The part below 320 * the diagonal is zero, because \p T is upper triangular. 321 */ 322 template <typename MatrixType, typename VectorType> 323 void matrix_function_compute_above_diagonal(const MatrixType& T, const VectorType& blockStart, const VectorType& clusterSize, MatrixType& fT) 324 { 325 typedef internal::traits<MatrixType> Traits; 326 typedef typename MatrixType::Scalar Scalar; 327 static const int Options = MatrixType::Options; 328 typedef Matrix<Scalar, Dynamic, Dynamic, Options, Traits::RowsAtCompileTime, Traits::ColsAtCompileTime> DynMatrixType; 329 330 for (Index k = 1; k < clusterSize.rows(); k++) { 331 for (Index i = 0; i < clusterSize.rows() - k; i++) { 332 // compute (i, i+k) block 333 DynMatrixType A = T.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i)); 334 DynMatrixType B = -T.block(blockStart(i+k), blockStart(i+k), clusterSize(i+k), clusterSize(i+k)); 335 DynMatrixType C = fT.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i)) 336 * T.block(blockStart(i), blockStart(i+k), clusterSize(i), clusterSize(i+k)); 337 C -= T.block(blockStart(i), blockStart(i+k), clusterSize(i), clusterSize(i+k)) 338 * fT.block(blockStart(i+k), blockStart(i+k), clusterSize(i+k), clusterSize(i+k)); 339 for (Index m = i + 1; m < i + k; m++) { 340 C += fT.block(blockStart(i), blockStart(m), clusterSize(i), clusterSize(m)) 341 * T.block(blockStart(m), blockStart(i+k), clusterSize(m), clusterSize(i+k)); 342 C -= T.block(blockStart(i), blockStart(m), clusterSize(i), clusterSize(m)) 343 * fT.block(blockStart(m), blockStart(i+k), clusterSize(m), clusterSize(i+k)); 344 } 345 fT.block(blockStart(i), blockStart(i+k), clusterSize(i), clusterSize(i+k)) 346 = matrix_function_solve_triangular_sylvester(A, B, C); 347 } 348 } 349 } 350 351 /** \ingroup MatrixFunctions_Module 352 * \brief Class for computing matrix functions. 353 * \tparam MatrixType type of the argument of the matrix function, 354 * expected to be an instantiation of the Matrix class template. 355 * \tparam AtomicType type for computing matrix function of atomic blocks. 356 * \tparam IsComplex used internally to select correct specialization. 357 * 358 * This class implements the Schur-Parlett algorithm for computing matrix functions. The spectrum of the 359 * matrix is divided in clustered of eigenvalues that lies close together. This class delegates the 360 * computation of the matrix function on every block corresponding to these clusters to an object of type 361 * \p AtomicType and uses these results to compute the matrix function of the whole matrix. The class 362 * \p AtomicType should have a \p compute() member function for computing the matrix function of a block. 363 * 364 * \sa class MatrixFunctionAtomic, class MatrixLogarithmAtomic 365 */ 366 template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex> 367 struct matrix_function_compute 368 { 369 /** \brief Compute the matrix function. 370 * 371 * \param[in] A argument of matrix function, should be a square matrix. 372 * \param[in] atomic class for computing matrix function of atomic blocks. 373 * \param[out] result the function \p f applied to \p A, as 374 * specified in the constructor. 375 * 376 * See MatrixBase::matrixFunction() for details on how this computation 377 * is implemented. 378 */ 379 template <typename AtomicType, typename ResultType> 380 static void run(const MatrixType& A, AtomicType& atomic, ResultType &result); 381 }; 382 383 /** \internal \ingroup MatrixFunctions_Module 384 * \brief Partial specialization of MatrixFunction for real matrices 385 * 386 * This converts the real matrix to a complex matrix, compute the matrix function of that matrix, and then 387 * converts the result back to a real matrix. 388 */ 389 template <typename MatrixType> 390 struct matrix_function_compute<MatrixType, 0> 391 { 392 template <typename MatA, typename AtomicType, typename ResultType> 393 static void run(const MatA& A, AtomicType& atomic, ResultType &result) 394 { 395 typedef internal::traits<MatrixType> Traits; 396 typedef typename Traits::Scalar Scalar; 397 static const int Rows = Traits::RowsAtCompileTime, Cols = Traits::ColsAtCompileTime; 398 static const int MaxRows = Traits::MaxRowsAtCompileTime, MaxCols = Traits::MaxColsAtCompileTime; 399 400 typedef std::complex<Scalar> ComplexScalar; 401 typedef Matrix<ComplexScalar, Rows, Cols, 0, MaxRows, MaxCols> ComplexMatrix; 402 403 ComplexMatrix CA = A.template cast<ComplexScalar>(); 404 ComplexMatrix Cresult; 405 matrix_function_compute<ComplexMatrix>::run(CA, atomic, Cresult); 406 result = Cresult.real(); 407 } 408 }; 409 410 /** \internal \ingroup MatrixFunctions_Module 411 * \brief Partial specialization of MatrixFunction for complex matrices 412 */ 413 template <typename MatrixType> 414 struct matrix_function_compute<MatrixType, 1> 415 { 416 template <typename MatA, typename AtomicType, typename ResultType> 417 static void run(const MatA& A, AtomicType& atomic, ResultType &result) 418 { 419 typedef internal::traits<MatrixType> Traits; 420 421 // compute Schur decomposition of A 422 const ComplexSchur<MatrixType> schurOfA(A); 423 eigen_assert(schurOfA.info()==Success); 424 MatrixType T = schurOfA.matrixT(); 425 MatrixType U = schurOfA.matrixU(); 426 427 // partition eigenvalues into clusters of ei'vals "close" to each other 428 std::list<std::list<Index> > clusters; 429 matrix_function_partition_eigenvalues(T.diagonal(), clusters); 430 431 // compute size of each cluster 432 Matrix<Index, Dynamic, 1> clusterSize; 433 matrix_function_compute_cluster_size(clusters, clusterSize); 434 435 // blockStart[i] is row index at which block corresponding to i-th cluster starts 436 Matrix<Index, Dynamic, 1> blockStart; 437 matrix_function_compute_block_start(clusterSize, blockStart); 438 439 // compute map so that eivalToCluster[i] = j means that i-th ei'val is in j-th cluster 440 Matrix<Index, Dynamic, 1> eivalToCluster; 441 matrix_function_compute_map(T.diagonal(), clusters, eivalToCluster); 442 443 // compute permutation which groups ei'vals in same cluster together 444 Matrix<Index, Traits::RowsAtCompileTime, 1> permutation; 445 matrix_function_compute_permutation(blockStart, eivalToCluster, permutation); 446 447 // permute Schur decomposition 448 matrix_function_permute_schur(permutation, U, T); 449 450 // compute result 451 MatrixType fT; // matrix function applied to T 452 matrix_function_compute_block_atomic(T, atomic, blockStart, clusterSize, fT); 453 matrix_function_compute_above_diagonal(T, blockStart, clusterSize, fT); 454 result = U * (fT.template triangularView<Upper>() * U.adjoint()); 455 } 456 }; 457 458 } // end of namespace internal 459 460 /** \ingroup MatrixFunctions_Module 461 * 462 * \brief Proxy for the matrix function of some matrix (expression). 463 * 464 * \tparam Derived Type of the argument to the matrix function. 465 * 466 * This class holds the argument to the matrix function until it is assigned or evaluated for some other 467 * reason (so the argument should not be changed in the meantime). It is the return type of 468 * matrixBase::matrixFunction() and related functions and most of the time this is the only way it is used. 469 */ 470 template<typename Derived> class MatrixFunctionReturnValue 471 : public ReturnByValue<MatrixFunctionReturnValue<Derived> > 472 { 473 public: 474 typedef typename Derived::Scalar Scalar; 475 typedef typename internal::stem_function<Scalar>::type StemFunction; 476 477 protected: 478 typedef typename internal::ref_selector<Derived>::type DerivedNested; 479 480 public: 481 482 /** \brief Constructor. 483 * 484 * \param[in] A %Matrix (expression) forming the argument of the matrix function. 485 * \param[in] f Stem function for matrix function under consideration. 486 */ 487 MatrixFunctionReturnValue(const Derived& A, StemFunction f) : m_A(A), m_f(f) { } 488 489 /** \brief Compute the matrix function. 490 * 491 * \param[out] result \p f applied to \p A, where \p f and \p A are as in the constructor. 492 */ 493 template <typename ResultType> 494 inline void evalTo(ResultType& result) const 495 { 496 typedef typename internal::nested_eval<Derived, 10>::type NestedEvalType; 497 typedef typename internal::remove_all<NestedEvalType>::type NestedEvalTypeClean; 498 typedef internal::traits<NestedEvalTypeClean> Traits; 499 typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar; 500 typedef Matrix<ComplexScalar, Dynamic, Dynamic, 0, Traits::RowsAtCompileTime, Traits::ColsAtCompileTime> DynMatrixType; 501 502 typedef internal::MatrixFunctionAtomic<DynMatrixType> AtomicType; 503 AtomicType atomic(m_f); 504 505 internal::matrix_function_compute<typename NestedEvalTypeClean::PlainObject>::run(m_A, atomic, result); 506 } 507 508 Index rows() const { return m_A.rows(); } 509 Index cols() const { return m_A.cols(); } 510 511 private: 512 const DerivedNested m_A; 513 StemFunction *m_f; 514 }; 515 516 namespace internal { 517 template<typename Derived> 518 struct traits<MatrixFunctionReturnValue<Derived> > 519 { 520 typedef typename Derived::PlainObject ReturnType; 521 }; 522 } 523 524 525 /********** MatrixBase methods **********/ 526 527 528 template <typename Derived> 529 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename internal::stem_function<typename internal::traits<Derived>::Scalar>::type f) const 530 { 531 eigen_assert(rows() == cols()); 532 return MatrixFunctionReturnValue<Derived>(derived(), f); 533 } 534 535 template <typename Derived> 536 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const 537 { 538 eigen_assert(rows() == cols()); 539 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar; 540 return MatrixFunctionReturnValue<Derived>(derived(), internal::stem_function_sin<ComplexScalar>); 541 } 542 543 template <typename Derived> 544 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const 545 { 546 eigen_assert(rows() == cols()); 547 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar; 548 return MatrixFunctionReturnValue<Derived>(derived(), internal::stem_function_cos<ComplexScalar>); 549 } 550 551 template <typename Derived> 552 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const 553 { 554 eigen_assert(rows() == cols()); 555 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar; 556 return MatrixFunctionReturnValue<Derived>(derived(), internal::stem_function_sinh<ComplexScalar>); 557 } 558 559 template <typename Derived> 560 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const 561 { 562 eigen_assert(rows() == cols()); 563 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar; 564 return MatrixFunctionReturnValue<Derived>(derived(), internal::stem_function_cosh<ComplexScalar>); 565 } 566 567 } // end namespace Eigen 568 569 #endif // EIGEN_MATRIX_FUNCTION_H