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 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_MATRIX_FUNCTION
26 #define EIGEN_MATRIX_FUNCTION
27 
28 #include "StemFunction.h"
29 #include "MatrixFunctionAtomic.h"
30 
31 
32 namespace Eigen {
33 
49 template <typename MatrixType,
50  typename AtomicType,
51  int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
53 {
54  public:
55 
64  MatrixFunction(const MatrixType& A, AtomicType& atomic);
65 
74  template <typename ResultType>
75  void compute(ResultType &result);
76 };
77 
78 
82 template <typename MatrixType, typename AtomicType>
83 class MatrixFunction<MatrixType, AtomicType, 0>
84 {
85  private:
86 
87  typedef internal::traits<MatrixType> Traits;
88  typedef typename Traits::Scalar Scalar;
89  static const int Rows = Traits::RowsAtCompileTime;
90  static const int Cols = Traits::ColsAtCompileTime;
91  static const int Options = MatrixType::Options;
92  static const int MaxRows = Traits::MaxRowsAtCompileTime;
93  static const int MaxCols = Traits::MaxColsAtCompileTime;
94 
95  typedef std::complex<Scalar> ComplexScalar;
96  typedef Matrix<ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols> ComplexMatrix;
97 
98  public:
99 
105  MatrixFunction(const MatrixType& A, AtomicType& atomic) : m_A(A), m_atomic(atomic) { }
106 
116  template <typename ResultType>
117  void compute(ResultType& result)
118  {
119  ComplexMatrix CA = m_A.template cast<ComplexScalar>();
120  ComplexMatrix Cresult;
121  MatrixFunction<ComplexMatrix, AtomicType> mf(CA, m_atomic);
122  mf.compute(Cresult);
123  result = Cresult.real();
124  }
125 
126  private:
127  typename internal::nested<MatrixType>::type m_A;
128  AtomicType& m_atomic;
130  MatrixFunction& operator=(const MatrixFunction&);
131 };
132 
133 
137 template <typename MatrixType, typename AtomicType>
138 class MatrixFunction<MatrixType, AtomicType, 1>
139 {
140  private:
141 
142  typedef internal::traits<MatrixType> Traits;
143  typedef typename MatrixType::Scalar Scalar;
144  typedef typename MatrixType::Index Index;
145  static const int RowsAtCompileTime = Traits::RowsAtCompileTime;
146  static const int ColsAtCompileTime = Traits::ColsAtCompileTime;
147  static const int Options = MatrixType::Options;
148  typedef typename NumTraits<Scalar>::Real RealScalar;
149  typedef Matrix<Scalar, Traits::RowsAtCompileTime, 1> VectorType;
150  typedef Matrix<Index, Traits::RowsAtCompileTime, 1> IntVectorType;
151  typedef Matrix<Index, Dynamic, 1> DynamicIntVectorType;
152  typedef std::list<Scalar> Cluster;
153  typedef std::list<Cluster> ListOfClusters;
154  typedef Matrix<Scalar, Dynamic, Dynamic, Options, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType;
155 
156  public:
157 
158  MatrixFunction(const MatrixType& A, AtomicType& atomic);
159  template <typename ResultType> void compute(ResultType& result);
160 
161  private:
162 
163  void computeSchurDecomposition();
164  void partitionEigenvalues();
165  typename ListOfClusters::iterator findCluster(Scalar key);
166  void computeClusterSize();
167  void computeBlockStart();
168  void constructPermutation();
169  void permuteSchur();
170  void swapEntriesInSchur(Index index);
171  void computeBlockAtomic();
172  Block<MatrixType> block(MatrixType& A, Index i, Index j);
173  void computeOffDiagonal();
174  DynMatrixType solveTriangularSylvester(const DynMatrixType& A, const DynMatrixType& B, const DynMatrixType& C);
175 
176  typename internal::nested<MatrixType>::type m_A;
177  AtomicType& m_atomic;
178  MatrixType m_T;
179  MatrixType m_U;
180  MatrixType m_fT;
181  ListOfClusters m_clusters;
182  DynamicIntVectorType m_eivalToCluster;
183  DynamicIntVectorType m_clusterSize;
184  DynamicIntVectorType m_blockStart;
185  IntVectorType m_permutation;
193  static const RealScalar separation() { return static_cast<RealScalar>(0.1); }
194 
195  MatrixFunction& operator=(const MatrixFunction&);
196 };
197 
203 template <typename MatrixType, typename AtomicType>
204 MatrixFunction<MatrixType,AtomicType,1>::MatrixFunction(const MatrixType& A, AtomicType& atomic)
205  : m_A(A), m_atomic(atomic)
206 {
207  /* empty body */
208 }
209 
215 template <typename MatrixType, typename AtomicType>
216 template <typename ResultType>
217 void MatrixFunction<MatrixType,AtomicType,1>::compute(ResultType& result)
218 {
219  computeSchurDecomposition();
220  partitionEigenvalues();
221  computeClusterSize();
222  computeBlockStart();
223  constructPermutation();
224  permuteSchur();
225  computeBlockAtomic();
226  computeOffDiagonal();
227  result = m_U * m_fT * m_U.adjoint();
228 }
229 
231 template <typename MatrixType, typename AtomicType>
232 void MatrixFunction<MatrixType,AtomicType,1>::computeSchurDecomposition()
233 {
234  const ComplexSchur<MatrixType> schurOfA(m_A);
235  m_T = schurOfA.matrixT();
236  m_U = schurOfA.matrixU();
237 }
238 
250 template <typename MatrixType, typename AtomicType>
251 void MatrixFunction<MatrixType,AtomicType,1>::partitionEigenvalues()
252 {
253  const Index rows = m_T.rows();
254  VectorType diag = m_T.diagonal(); // contains eigenvalues of A
255 
256  for (Index i=0; i<rows; ++i) {
257  // Find set containing diag(i), adding a new set if necessary
258  typename ListOfClusters::iterator qi = findCluster(diag(i));
259  if (qi == m_clusters.end()) {
260  Cluster l;
261  l.push_back(diag(i));
262  m_clusters.push_back(l);
263  qi = m_clusters.end();
264  --qi;
265  }
266 
267  // Look for other element to add to the set
268  for (Index j=i+1; j<rows; ++j) {
269  if (internal::abs(diag(j) - diag(i)) <= separation() && std::find(qi->begin(), qi->end(), diag(j)) == qi->end()) {
270  typename ListOfClusters::iterator qj = findCluster(diag(j));
271  if (qj == m_clusters.end()) {
272  qi->push_back(diag(j));
273  } else {
274  qi->insert(qi->end(), qj->begin(), qj->end());
275  m_clusters.erase(qj);
276  }
277  }
278  }
279  }
280 }
281 
287 template <typename MatrixType, typename AtomicType>
288 typename MatrixFunction<MatrixType,AtomicType,1>::ListOfClusters::iterator MatrixFunction<MatrixType,AtomicType,1>::findCluster(Scalar key)
289 {
290  typename Cluster::iterator j;
291  for (typename ListOfClusters::iterator i = m_clusters.begin(); i != m_clusters.end(); ++i) {
292  j = std::find(i->begin(), i->end(), key);
293  if (j != i->end())
294  return i;
295  }
296  return m_clusters.end();
297 }
298 
300 template <typename MatrixType, typename AtomicType>
301 void MatrixFunction<MatrixType,AtomicType,1>::computeClusterSize()
302 {
303  const Index rows = m_T.rows();
304  VectorType diag = m_T.diagonal();
305  const Index numClusters = static_cast<Index>(m_clusters.size());
306 
307  m_clusterSize.setZero(numClusters);
308  m_eivalToCluster.resize(rows);
309  Index clusterIndex = 0;
310  for (typename ListOfClusters::const_iterator cluster = m_clusters.begin(); cluster != m_clusters.end(); ++cluster) {
311  for (Index i = 0; i < diag.rows(); ++i) {
312  if (std::find(cluster->begin(), cluster->end(), diag(i)) != cluster->end()) {
313  ++m_clusterSize[clusterIndex];
314  m_eivalToCluster[i] = clusterIndex;
315  }
316  }
317  ++clusterIndex;
318  }
319 }
320 
322 template <typename MatrixType, typename AtomicType>
323 void MatrixFunction<MatrixType,AtomicType,1>::computeBlockStart()
324 {
325  m_blockStart.resize(m_clusterSize.rows());
326  m_blockStart(0) = 0;
327  for (Index i = 1; i < m_clusterSize.rows(); i++) {
328  m_blockStart(i) = m_blockStart(i-1) + m_clusterSize(i-1);
329  }
330 }
331 
333 template <typename MatrixType, typename AtomicType>
334 void MatrixFunction<MatrixType,AtomicType,1>::constructPermutation()
335 {
336  DynamicIntVectorType indexNextEntry = m_blockStart;
337  m_permutation.resize(m_T.rows());
338  for (Index i = 0; i < m_T.rows(); i++) {
339  Index cluster = m_eivalToCluster[i];
340  m_permutation[i] = indexNextEntry[cluster];
341  ++indexNextEntry[cluster];
342  }
343 }
344 
346 template <typename MatrixType, typename AtomicType>
347 void MatrixFunction<MatrixType,AtomicType,1>::permuteSchur()
348 {
349  IntVectorType p = m_permutation;
350  for (Index i = 0; i < p.rows() - 1; i++) {
351  Index j;
352  for (j = i; j < p.rows(); j++) {
353  if (p(j) == i) break;
354  }
355  eigen_assert(p(j) == i);
356  for (Index k = j-1; k >= i; k--) {
357  swapEntriesInSchur(k);
358  std::swap(p.coeffRef(k), p.coeffRef(k+1));
359  }
360  }
361 }
362 
364 template <typename MatrixType, typename AtomicType>
365 void MatrixFunction<MatrixType,AtomicType,1>::swapEntriesInSchur(Index index)
366 {
367  JacobiRotation<Scalar> rotation;
368  rotation.makeGivens(m_T(index, index+1), m_T(index+1, index+1) - m_T(index, index));
369  m_T.applyOnTheLeft(index, index+1, rotation.adjoint());
370  m_T.applyOnTheRight(index, index+1, rotation);
371  m_U.applyOnTheRight(index, index+1, rotation);
372 }
373 
380 template <typename MatrixType, typename AtomicType>
381 void MatrixFunction<MatrixType,AtomicType,1>::computeBlockAtomic()
382 {
383  m_fT.resize(m_T.rows(), m_T.cols());
384  m_fT.setZero();
385  for (Index i = 0; i < m_clusterSize.rows(); ++i) {
386  block(m_fT, i, i) = m_atomic.compute(block(m_T, i, i));
387  }
388 }
389 
391 template <typename MatrixType, typename AtomicType>
392 Block<MatrixType> MatrixFunction<MatrixType,AtomicType,1>::block(MatrixType& A, Index i, Index j)
393 {
394  return A.block(m_blockStart(i), m_blockStart(j), m_clusterSize(i), m_clusterSize(j));
395 }
396 
404 template <typename MatrixType, typename AtomicType>
405 void MatrixFunction<MatrixType,AtomicType,1>::computeOffDiagonal()
406 {
407  for (Index diagIndex = 1; diagIndex < m_clusterSize.rows(); diagIndex++) {
408  for (Index blockIndex = 0; blockIndex < m_clusterSize.rows() - diagIndex; blockIndex++) {
409  // compute (blockIndex, blockIndex+diagIndex) block
410  DynMatrixType A = block(m_T, blockIndex, blockIndex);
411  DynMatrixType B = -block(m_T, blockIndex+diagIndex, blockIndex+diagIndex);
412  DynMatrixType C = block(m_fT, blockIndex, blockIndex) * block(m_T, blockIndex, blockIndex+diagIndex);
413  C -= block(m_T, blockIndex, blockIndex+diagIndex) * block(m_fT, blockIndex+diagIndex, blockIndex+diagIndex);
414  for (Index k = blockIndex + 1; k < blockIndex + diagIndex; k++) {
415  C += block(m_fT, blockIndex, k) * block(m_T, k, blockIndex+diagIndex);
416  C -= block(m_T, blockIndex, k) * block(m_fT, k, blockIndex+diagIndex);
417  }
418  block(m_fT, blockIndex, blockIndex+diagIndex) = solveTriangularSylvester(A, B, C);
419  }
420  }
421 }
422 
446 template <typename MatrixType, typename AtomicType>
447 typename MatrixFunction<MatrixType,AtomicType,1>::DynMatrixType MatrixFunction<MatrixType,AtomicType,1>::solveTriangularSylvester(
448  const DynMatrixType& A,
449  const DynMatrixType& B,
450  const DynMatrixType& C)
451 {
452  eigen_assert(A.rows() == A.cols());
453  eigen_assert(A.isUpperTriangular());
454  eigen_assert(B.rows() == B.cols());
455  eigen_assert(B.isUpperTriangular());
456  eigen_assert(C.rows() == A.rows());
457  eigen_assert(C.cols() == B.rows());
458 
459  Index m = A.rows();
460  Index n = B.rows();
461  DynMatrixType X(m, n);
462 
463  for (Index i = m - 1; i >= 0; --i) {
464  for (Index j = 0; j < n; ++j) {
465 
466  // Compute AX = \sum_{k=i+1}^m A_{ik} X_{kj}
467  Scalar AX;
468  if (i == m - 1) {
469  AX = 0;
470  } else {
471  Matrix<Scalar,1,1> AXmatrix = A.row(i).tail(m-1-i) * X.col(j).tail(m-1-i);
472  AX = AXmatrix(0,0);
473  }
474 
475  // Compute XB = \sum_{k=1}^{j-1} X_{ik} B_{kj}
476  Scalar XB;
477  if (j == 0) {
478  XB = 0;
479  } else {
480  Matrix<Scalar,1,1> XBmatrix = X.row(i).head(j) * B.col(j).head(j);
481  XB = XBmatrix(0,0);
482  }
483 
484  X(i,j) = (C(i,j) - AX - XB) / (A(i,i) + B(j,j));
485  }
486  }
487  return X;
488 }
489 
502 template<typename Derived> class MatrixFunctionReturnValue
503 : public ReturnByValue<MatrixFunctionReturnValue<Derived> >
504 {
505  public:
506 
507  typedef typename Derived::Scalar Scalar;
508  typedef typename Derived::Index Index;
509  typedef typename internal::stem_function<Scalar>::type StemFunction;
510 
517  MatrixFunctionReturnValue(const Derived& A, StemFunction f) : m_A(A), m_f(f) { }
518 
524  template <typename ResultType>
525  inline void evalTo(ResultType& result) const
526  {
527  typedef typename Derived::PlainObject PlainObject;
528  typedef internal::traits<PlainObject> Traits;
529  static const int RowsAtCompileTime = Traits::RowsAtCompileTime;
530  static const int ColsAtCompileTime = Traits::ColsAtCompileTime;
531  static const int Options = PlainObject::Options;
532  typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
533  typedef Matrix<ComplexScalar, Dynamic, Dynamic, Options, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType;
534  typedef MatrixFunctionAtomic<DynMatrixType> AtomicType;
535  AtomicType atomic(m_f);
536 
537  const PlainObject Aevaluated = m_A.eval();
538  MatrixFunction<PlainObject, AtomicType> mf(Aevaluated, atomic);
539  mf.compute(result);
540  }
541 
542  Index rows() const { return m_A.rows(); }
543  Index cols() const { return m_A.cols(); }
544 
545  private:
546  typename internal::nested<Derived>::type m_A;
547  StemFunction *m_f;
548 
550 };
551 
552 namespace internal {
553 template<typename Derived>
554 struct traits<MatrixFunctionReturnValue<Derived> >
555 {
556  typedef typename Derived::PlainObject ReturnType;
557 };
558 }
559 
560 
561 /********** MatrixBase methods **********/
562 
563 
564 template <typename Derived>
565 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename internal::stem_function<typename internal::traits<Derived>::Scalar>::type f) const
566 {
567  eigen_assert(rows() == cols());
568  return MatrixFunctionReturnValue<Derived>(derived(), f);
569 }
570 
571 template <typename Derived>
572 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const
573 {
574  eigen_assert(rows() == cols());
575  typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
576  return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::sin);
577 }
578 
579 template <typename Derived>
580 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const
581 {
582  eigen_assert(rows() == cols());
583  typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
584  return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::cos);
585 }
586 
587 template <typename Derived>
588 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const
589 {
590  eigen_assert(rows() == cols());
591  typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
592  return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::sinh);
593 }
594 
595 template <typename Derived>
596 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const
597 {
598  eigen_assert(rows() == cols());
599  typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
600  return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::cosh);
601 }
602 
603 } // end namespace Eigen
604 
605 #endif // EIGEN_MATRIX_FUNCTION