MatrixFunction.h
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2009-2011 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#include "MatrixFunctionAtomic.h"
15
16
17namespace Eigen {
18
34template <typename MatrixType,
35 typename AtomicType,
36 int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
37class MatrixFunction
38{
39 public:
40
49 MatrixFunction(const MatrixType& A, AtomicType& atomic);
50
59 template <typename ResultType>
60 void compute(ResultType &result);
61};
62
63
67template <typename MatrixType, typename AtomicType>
68class MatrixFunction<MatrixType, AtomicType, 0>
69{
70 private:
71
72 typedef internal::traits<MatrixType> Traits;
73 typedef typename Traits::Scalar Scalar;
74 static const int Rows = Traits::RowsAtCompileTime;
75 static const int Cols = Traits::ColsAtCompileTime;
76 static const int Options = MatrixType::Options;
77 static const int MaxRows = Traits::MaxRowsAtCompileTime;
78 static const int MaxCols = Traits::MaxColsAtCompileTime;
79
80 typedef std::complex<Scalar> ComplexScalar;
81 typedef Matrix<ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols> ComplexMatrix;
82
83 public:
84
90 MatrixFunction(const MatrixType& A, AtomicType& atomic) : m_A(A), m_atomic(atomic) { }
91
101 template <typename ResultType>
102 void compute(ResultType& result)
103 {
104 ComplexMatrix CA = m_A.template cast<ComplexScalar>();
105 ComplexMatrix Cresult;
106 MatrixFunction<ComplexMatrix, AtomicType> mf(CA, m_atomic);
107 mf.compute(Cresult);
108 result = Cresult.real();
109 }
110
111 private:
112 typename internal::nested<MatrixType>::type m_A;
113 AtomicType& m_atomic;
114
115 MatrixFunction& operator=(const MatrixFunction&);
116};
117
118
122template <typename MatrixType, typename AtomicType>
123class MatrixFunction<MatrixType, AtomicType, 1>
124{
125 private:
126
127 typedef internal::traits<MatrixType> Traits;
128 typedef typename MatrixType::Scalar Scalar;
129 typedef typename MatrixType::Index Index;
130 static const int RowsAtCompileTime = Traits::RowsAtCompileTime;
131 static const int ColsAtCompileTime = Traits::ColsAtCompileTime;
132 static const int Options = MatrixType::Options;
133 typedef typename NumTraits<Scalar>::Real RealScalar;
134 typedef Matrix<Scalar, Traits::RowsAtCompileTime, 1> VectorType;
135 typedef Matrix<Index, Traits::RowsAtCompileTime, 1> IntVectorType;
136 typedef Matrix<Index, Dynamic, 1> DynamicIntVectorType;
137 typedef std::list<Scalar> Cluster;
138 typedef std::list<Cluster> ListOfClusters;
139 typedef Matrix<Scalar, Dynamic, Dynamic, Options, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType;
140
141 public:
142
143 MatrixFunction(const MatrixType& A, AtomicType& atomic);
144 template <typename ResultType> void compute(ResultType& result);
145
146 private:
147
148 void computeSchurDecomposition();
149 void partitionEigenvalues();
150 typename ListOfClusters::iterator findCluster(Scalar key);
151 void computeClusterSize();
152 void computeBlockStart();
153 void constructPermutation();
154 void permuteSchur();
155 void swapEntriesInSchur(Index index);
156 void computeBlockAtomic();
157 Block<MatrixType> block(MatrixType& A, Index i, Index j);
158 void computeOffDiagonal();
159 DynMatrixType solveTriangularSylvester(const DynMatrixType& A, const DynMatrixType& B, const DynMatrixType& C);
160
161 typename internal::nested<MatrixType>::type m_A;
162 AtomicType& m_atomic;
163 MatrixType m_T;
164 MatrixType m_U;
165 MatrixType m_fT;
166 ListOfClusters m_clusters;
167 DynamicIntVectorType m_eivalToCluster;
168 DynamicIntVectorType m_clusterSize;
169 DynamicIntVectorType m_blockStart;
170 IntVectorType m_permutation;
171
178 static const RealScalar separation() { return static_cast<RealScalar>(0.1); }
179
180 MatrixFunction& operator=(const MatrixFunction&);
181};
182
188template <typename MatrixType, typename AtomicType>
189MatrixFunction<MatrixType,AtomicType,1>::MatrixFunction(const MatrixType& A, AtomicType& atomic)
190 : m_A(A), m_atomic(atomic)
191{
192 /* empty body */
193}
194
200template <typename MatrixType, typename AtomicType>
201template <typename ResultType>
202void MatrixFunction<MatrixType,AtomicType,1>::compute(ResultType& result)
203{
204 computeSchurDecomposition();
205 partitionEigenvalues();
206 computeClusterSize();
207 computeBlockStart();
208 constructPermutation();
209 permuteSchur();
210 computeBlockAtomic();
211 computeOffDiagonal();
212 result = m_U * m_fT * m_U.adjoint();
213}
214
216template <typename MatrixType, typename AtomicType>
217void MatrixFunction<MatrixType,AtomicType,1>::computeSchurDecomposition()
218{
219 const ComplexSchur<MatrixType> schurOfA(m_A);
220 m_T = schurOfA.matrixT();
221 m_U = schurOfA.matrixU();
222}
223
235template <typename MatrixType, typename AtomicType>
236void MatrixFunction<MatrixType,AtomicType,1>::partitionEigenvalues()
237{
238 const Index rows = m_T.rows();
239 VectorType diag = m_T.diagonal(); // contains eigenvalues of A
240
241 for (Index i=0; i<rows; ++i) {
242 // Find set containing diag(i), adding a new set if necessary
243 typename ListOfClusters::iterator qi = findCluster(diag(i));
244 if (qi == m_clusters.end()) {
245 Cluster l;
246 l.push_back(diag(i));
247 m_clusters.push_back(l);
248 qi = m_clusters.end();
249 --qi;
250 }
251
252 // Look for other element to add to the set
253 for (Index j=i+1; j<rows; ++j) {
254 if (internal::abs(diag(j) - diag(i)) <= separation() && std::find(qi->begin(), qi->end(), diag(j)) == qi->end()) {
255 typename ListOfClusters::iterator qj = findCluster(diag(j));
256 if (qj == m_clusters.end()) {
257 qi->push_back(diag(j));
258 } else {
259 qi->insert(qi->end(), qj->begin(), qj->end());
260 m_clusters.erase(qj);
261 }
262 }
263 }
264 }
265}
266
272template <typename MatrixType, typename AtomicType>
273typename MatrixFunction<MatrixType,AtomicType,1>::ListOfClusters::iterator MatrixFunction<MatrixType,AtomicType,1>::findCluster(Scalar key)
274{
275 typename Cluster::iterator j;
276 for (typename ListOfClusters::iterator i = m_clusters.begin(); i != m_clusters.end(); ++i) {
277 j = std::find(i->begin(), i->end(), key);
278 if (j != i->end())
279 return i;
280 }
281 return m_clusters.end();
282}
283
285template <typename MatrixType, typename AtomicType>
286void MatrixFunction<MatrixType,AtomicType,1>::computeClusterSize()
287{
288 const Index rows = m_T.rows();
289 VectorType diag = m_T.diagonal();
290 const Index numClusters = static_cast<Index>(m_clusters.size());
291
292 m_clusterSize.setZero(numClusters);
293 m_eivalToCluster.resize(rows);
294 Index clusterIndex = 0;
295 for (typename ListOfClusters::const_iterator cluster = m_clusters.begin(); cluster != m_clusters.end(); ++cluster) {
296 for (Index i = 0; i < diag.rows(); ++i) {
297 if (std::find(cluster->begin(), cluster->end(), diag(i)) != cluster->end()) {
298 ++m_clusterSize[clusterIndex];
299 m_eivalToCluster[i] = clusterIndex;
300 }
301 }
302 ++clusterIndex;
303 }
304}
305
307template <typename MatrixType, typename AtomicType>
308void MatrixFunction<MatrixType,AtomicType,1>::computeBlockStart()
309{
310 m_blockStart.resize(m_clusterSize.rows());
311 m_blockStart(0) = 0;
312 for (Index i = 1; i < m_clusterSize.rows(); i++) {
313 m_blockStart(i) = m_blockStart(i-1) + m_clusterSize(i-1);
314 }
315}
316
318template <typename MatrixType, typename AtomicType>
319void MatrixFunction<MatrixType,AtomicType,1>::constructPermutation()
320{
321 DynamicIntVectorType indexNextEntry = m_blockStart;
322 m_permutation.resize(m_T.rows());
323 for (Index i = 0; i < m_T.rows(); i++) {
324 Index cluster = m_eivalToCluster[i];
325 m_permutation[i] = indexNextEntry[cluster];
326 ++indexNextEntry[cluster];
327 }
328}
329
331template <typename MatrixType, typename AtomicType>
332void MatrixFunction<MatrixType,AtomicType,1>::permuteSchur()
333{
334 IntVectorType p = m_permutation;
335 for (Index i = 0; i < p.rows() - 1; i++) {
336 Index j;
337 for (j = i; j < p.rows(); j++) {
338 if (p(j) == i) break;
339 }
340 eigen_assert(p(j) == i);
341 for (Index k = j-1; k >= i; k--) {
342 swapEntriesInSchur(k);
343 std::swap(p.coeffRef(k), p.coeffRef(k+1));
344 }
345 }
346}
347
349template <typename MatrixType, typename AtomicType>
350void MatrixFunction<MatrixType,AtomicType,1>::swapEntriesInSchur(Index index)
351{
352 JacobiRotation<Scalar> rotation;
353 rotation.makeGivens(m_T(index, index+1), m_T(index+1, index+1) - m_T(index, index));
354 m_T.applyOnTheLeft(index, index+1, rotation.adjoint());
355 m_T.applyOnTheRight(index, index+1, rotation);
356 m_U.applyOnTheRight(index, index+1, rotation);
357}
358
365template <typename MatrixType, typename AtomicType>
366void MatrixFunction<MatrixType,AtomicType,1>::computeBlockAtomic()
367{
368 m_fT.resize(m_T.rows(), m_T.cols());
369 m_fT.setZero();
370 for (Index i = 0; i < m_clusterSize.rows(); ++i) {
371 block(m_fT, i, i) = m_atomic.compute(block(m_T, i, i));
372 }
373}
374
376template <typename MatrixType, typename AtomicType>
377Block<MatrixType> MatrixFunction<MatrixType,AtomicType,1>::block(MatrixType& A, Index i, Index j)
378{
379 return A.block(m_blockStart(i), m_blockStart(j), m_clusterSize(i), m_clusterSize(j));
380}
381
389template <typename MatrixType, typename AtomicType>
390void MatrixFunction<MatrixType,AtomicType,1>::computeOffDiagonal()
391{
392 for (Index diagIndex = 1; diagIndex < m_clusterSize.rows(); diagIndex++) {
393 for (Index blockIndex = 0; blockIndex < m_clusterSize.rows() - diagIndex; blockIndex++) {
394 // compute (blockIndex, blockIndex+diagIndex) block
395 DynMatrixType A = block(m_T, blockIndex, blockIndex);
396 DynMatrixType B = -block(m_T, blockIndex+diagIndex, blockIndex+diagIndex);
397 DynMatrixType C = block(m_fT, blockIndex, blockIndex) * block(m_T, blockIndex, blockIndex+diagIndex);
398 C -= block(m_T, blockIndex, blockIndex+diagIndex) * block(m_fT, blockIndex+diagIndex, blockIndex+diagIndex);
399 for (Index k = blockIndex + 1; k < blockIndex + diagIndex; k++) {
400 C += block(m_fT, blockIndex, k) * block(m_T, k, blockIndex+diagIndex);
401 C -= block(m_T, blockIndex, k) * block(m_fT, k, blockIndex+diagIndex);
402 }
403 block(m_fT, blockIndex, blockIndex+diagIndex) = solveTriangularSylvester(A, B, C);
404 }
405 }
406}
407
431template <typename MatrixType, typename AtomicType>
432typename MatrixFunction<MatrixType,AtomicType,1>::DynMatrixType MatrixFunction<MatrixType,AtomicType,1>::solveTriangularSylvester(
433 const DynMatrixType& A,
434 const DynMatrixType& B,
435 const DynMatrixType& C)
436{
437 eigen_assert(A.rows() == A.cols());
438 eigen_assert(A.isUpperTriangular());
439 eigen_assert(B.rows() == B.cols());
440 eigen_assert(B.isUpperTriangular());
441 eigen_assert(C.rows() == A.rows());
442 eigen_assert(C.cols() == B.rows());
443
444 Index m = A.rows();
445 Index n = B.rows();
446 DynMatrixType X(m, n);
447
448 for (Index i = m - 1; i >= 0; --i) {
449 for (Index j = 0; j < n; ++j) {
450
451 // Compute AX = \sum_{k=i+1}^m A_{ik} X_{kj}
452 Scalar AX;
453 if (i == m - 1) {
454 AX = 0;
455 } else {
456 Matrix<Scalar,1,1> AXmatrix = A.row(i).tail(m-1-i) * X.col(j).tail(m-1-i);
457 AX = AXmatrix(0,0);
458 }
459
460 // Compute XB = \sum_{k=1}^{j-1} X_{ik} B_{kj}
461 Scalar XB;
462 if (j == 0) {
463 XB = 0;
464 } else {
465 Matrix<Scalar,1,1> XBmatrix = X.row(i).head(j) * B.col(j).head(j);
466 XB = XBmatrix(0,0);
467 }
468
469 X(i,j) = (C(i,j) - AX - XB) / (A(i,i) + B(j,j));
470 }
471 }
472 return X;
473}
474
487template<typename Derived> class MatrixFunctionReturnValue
488: public ReturnByValue<MatrixFunctionReturnValue<Derived> >
489{
490 public:
491
492 typedef typename Derived::Scalar Scalar;
493 typedef typename Derived::Index Index;
494 typedef typename internal::stem_function<Scalar>::type StemFunction;
495
502 MatrixFunctionReturnValue(const Derived& A, StemFunction f) : m_A(A), m_f(f) { }
503
509 template <typename ResultType>
510 inline void evalTo(ResultType& result) const
511 {
512 typedef typename Derived::PlainObject PlainObject;
513 typedef internal::traits<PlainObject> Traits;
514 static const int RowsAtCompileTime = Traits::RowsAtCompileTime;
515 static const int ColsAtCompileTime = Traits::ColsAtCompileTime;
516 static const int Options = PlainObject::Options;
517 typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
518 typedef Matrix<ComplexScalar, Dynamic, Dynamic, Options, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType;
519 typedef MatrixFunctionAtomic<DynMatrixType> AtomicType;
520 AtomicType atomic(m_f);
521
522 const PlainObject Aevaluated = m_A.eval();
523 MatrixFunction<PlainObject, AtomicType> mf(Aevaluated, atomic);
524 mf.compute(result);
525 }
526
527 Index rows() const { return m_A.rows(); }
528 Index cols() const { return m_A.cols(); }
529
530 private:
531 typename internal::nested<Derived>::type m_A;
532 StemFunction *m_f;
533
534 MatrixFunctionReturnValue& operator=(const MatrixFunctionReturnValue&);
535};
536
537namespace internal {
538template<typename Derived>
539struct traits<MatrixFunctionReturnValue<Derived> >
540{
541 typedef typename Derived::PlainObject ReturnType;
542};
543}
544
545
546/********** MatrixBase methods **********/
547
548
549template <typename Derived>
550const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename internal::stem_function<typename internal::traits<Derived>::Scalar>::type f) const
551{
552 eigen_assert(rows() == cols());
553 return MatrixFunctionReturnValue<Derived>(derived(), f);
554}
555
556template <typename Derived>
557const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const
558{
559 eigen_assert(rows() == cols());
560 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
561 return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::sin);
562}
563
564template <typename Derived>
565const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const
566{
567 eigen_assert(rows() == cols());
568 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
569 return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::cos);
570}
571
572template <typename Derived>
573const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const
574{
575 eigen_assert(rows() == cols());
576 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
577 return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::sinh);
578}
579
580template <typename Derived>
581const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const
582{
583 eigen_assert(rows() == cols());
584 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
585 return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::cosh);
586}
587
588} // end namespace Eigen
589
590#endif // EIGEN_MATRIX_FUNCTION
Eigen::Matrix::row
RowXpr row(Index i)
Eigen::JacobiRotation::makeGivens
void makeGivens(const Scalar &p, const Scalar &q, Scalar *z=0)
Eigen::MatrixFunctionAtomic
Helper class for computing matrix functions of atomic matrices.
Definition MatrixFunctionAtomic.h:25
Eigen::MatrixFunctionReturnValue
Proxy for the matrix function of some matrix (expression).
Definition MatrixFunction.h:489
Eigen::MatrixFunctionReturnValue::MatrixFunctionReturnValue
MatrixFunctionReturnValue(const Derived &A, StemFunction f)
Constructor.
Definition MatrixFunction.h:502
Eigen::MatrixFunctionReturnValue::evalTo
void evalTo(ResultType &result) const
Compute the matrix function.
Definition MatrixFunction.h:510
Eigen::MatrixFunction
Class for computing matrix functions.
Definition MatrixFunction.h:38
Eigen::MatrixFunction::compute
void compute(ResultType &result)
Compute the matrix function.
Eigen::MatrixFunction::MatrixFunction
MatrixFunction(const MatrixType &A, AtomicType &atomic)
Constructor.
Eigen::StdStemFunctions::cos
static Scalar cos(Scalar x, int n)
Cosine (and its derivatives).
Definition StemFunction.h:30
Eigen::StdStemFunctions::cosh
static Scalar cosh(Scalar x, int n)
Hyperbolic cosine (and its derivatives).
Definition StemFunction.h:72
Eigen::StdStemFunctions::sin
static Scalar sin(Scalar x, int n)
Sine (and its derivatives).
Definition StemFunction.h:51
Eigen::StdStemFunctions::sinh
static Scalar sinh(Scalar x, int n)
Hyperbolic sine (and its derivatives).
Definition StemFunction.h:87
Generated by doxygen 1.13.2