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