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 
17 namespace Eigen {
18 
34 template <typename MatrixType,
35  typename AtomicType,
36  int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
38 {
39  public:
40 
49  MatrixFunction(const MatrixType& A, AtomicType& atomic);
50 
59  template <typename ResultType>
60  void compute(ResultType &result);
61 };
62 
63 
67 template <typename MatrixType, typename AtomicType>
68 class 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;
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;
115  MatrixFunction& operator=(const MatrixFunction&);
116 };
117 
118 
122 template <typename MatrixType, typename AtomicType>
123 class 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;
178  static const RealScalar separation() { return static_cast<RealScalar>(0.1); }
179 
180  MatrixFunction& operator=(const MatrixFunction&);
181 };
182 
188 template <typename MatrixType, typename AtomicType>
189 MatrixFunction<MatrixType,AtomicType,1>::MatrixFunction(const MatrixType& A, AtomicType& atomic)
190  : m_A(A), m_atomic(atomic)
191 {
192  /* empty body */
193 }
194 
200 template <typename MatrixType, typename AtomicType>
201 template <typename ResultType>
202 void 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 
216 template <typename MatrixType, typename AtomicType>
217 void 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 
235 template <typename MatrixType, typename AtomicType>
236 void 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 
272 template <typename MatrixType, typename AtomicType>
273 typename 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 
285 template <typename MatrixType, typename AtomicType>
286 void 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 
307 template <typename MatrixType, typename AtomicType>
308 void 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 
318 template <typename MatrixType, typename AtomicType>
319 void 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 
331 template <typename MatrixType, typename AtomicType>
332 void 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 
349 template <typename MatrixType, typename AtomicType>
350 void 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 
365 template <typename MatrixType, typename AtomicType>
366 void 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 
376 template <typename MatrixType, typename AtomicType>
377 Block<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 
389 template <typename MatrixType, typename AtomicType>
390 void 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 
431 template <typename MatrixType, typename AtomicType>
432 typename 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 
487 template<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;
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 
535 };
536 
537 namespace internal {
538 template<typename Derived>
539 struct traits<MatrixFunctionReturnValue<Derived> >
540 {
541  typedef typename Derived::PlainObject ReturnType;
542 };
543 }
544 
545 
546 /********** MatrixBase methods **********/
547 
548 
549 template <typename Derived>
550 const 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 
556 template <typename Derived>
557 const 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 
564 template <typename Derived>
565 const 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 
572 template <typename Derived>
573 const 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 
580 template <typename Derived>
581 const 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