• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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 
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       */
MatrixFunctionAtomic(StemFunction f)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>
matrix_function_compute_mu(const MatrixType & A)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   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 
63 template <typename MatrixType>
compute(const MatrixType & A)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
105   * \returns Iterator to cluster containing \p key, or \c clusters.end() if no cluster in \p m_clusters
106   * contains \p key.
107   */
108 template <typename Index, typename ListOfClusters>
matrix_function_find_cluster(Index key,ListOfClusters & clusters)109 typename 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   */
131 template <typename EivalsType, typename Cluster>
matrix_function_partition_eigenvalues(const EivalsType & eivals,std::list<Cluster> & clusters)132 void 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 */
164 template <typename ListOfClusters, typename Index>
matrix_function_compute_cluster_size(const ListOfClusters & clusters,Matrix<Index,Dynamic,1> & clusterSize)165 void 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 */
177 template <typename VectorType>
matrix_function_compute_block_start(const VectorType & clusterSize,VectorType & blockStart)178 void 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 */
188 template <typename EivalsType, typename ListOfClusters, typename VectorType>
matrix_function_compute_map(const EivalsType & eivals,const ListOfClusters & clusters,VectorType & eivalToCluster)189 void 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 */
205 template <typename DynVectorType, typename VectorType>
matrix_function_compute_permutation(const DynVectorType & blockStart,const DynVectorType & eivalToCluster,VectorType & permutation)206 void 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 */
219 template <typename VectorType, typename MatrixType>
matrix_function_permute_schur(VectorType & permutation,MatrixType & U,MatrixType & T)220 void 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   */
246 template <typename MatrixType, typename AtomicType, typename VectorType>
matrix_function_compute_block_atomic(const MatrixType & T,AtomicType & atomic,const VectorType & blockStart,const VectorType & clusterSize,MatrixType & fT)247 void 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   */
278 template <typename MatrixType>
matrix_function_solve_triangular_sylvester(const MatrixType & A,const MatrixType & B,const MatrixType & C)279 MatrixType 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   */
328 template <typename MatrixType, typename VectorType>
matrix_function_compute_above_diagonal(const MatrixType & T,const VectorType & blockStart,const VectorType & clusterSize,MatrixType & fT)329 void 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   */
375 template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
376 struct 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   */
398 template <typename MatrixType>
399 struct 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   */
422 template <typename MatrixType>
423 struct 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   */
478 template<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 
527 namespace internal {
528 template<typename Derived>
529 struct traits<MatrixFunctionReturnValue<Derived> >
530 {
531   typedef typename Derived::PlainObject ReturnType;
532 };
533 }
534 
535 
536 /********** MatrixBase methods **********/
537 
538 
539 template <typename Derived>
540 const 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 
546 template <typename Derived>
547 const 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 
554 template <typename Derived>
555 const 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 
562 template <typename Derived>
563 const 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 
570 template <typename Derived>
571 const 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