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