10#ifndef EIGEN_MATRIX_FUNCTION_H
11#define EIGEN_MATRIX_FUNCTION_H
13#include "StemFunction.h"
21static const float matrix_function_separation = 0.1f;
29template <
typename MatrixType>
30class MatrixFunctionAtomic
34 typedef typename MatrixType::Scalar Scalar;
35 typedef typename stem_function<Scalar>::type StemFunction;
40 MatrixFunctionAtomic(StemFunction f) : m_f(f) { }
46 MatrixType compute(
const MatrixType& A);
52template <
typename MatrixType>
53typename NumTraits<typename MatrixType::Scalar>::Real matrix_function_compute_mu(
const MatrixType& A)
55 typedef typename plain_col_type<MatrixType>::type VectorType;
56 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();
63template <
typename MatrixType>
64MatrixType MatrixFunctionAtomic<MatrixType>::compute(
const MatrixType& A)
67 typedef typename NumTraits<Scalar>::Real RealScalar;
68 Index rows = A.rows();
69 Scalar avgEival = A.trace() / Scalar(RealScalar(rows));
70 MatrixType Ashifted = A - avgEival * MatrixType::Identity(rows, rows);
71 RealScalar mu = matrix_function_compute_mu(Ashifted);
72 MatrixType F = m_f(avgEival, 0) * MatrixType::Identity(rows, rows);
73 MatrixType P = Ashifted;
75 for (
Index s = 1; double(s) < 1.1 * double(rows) + 10.0; s++) {
76 Fincr = m_f(avgEival,
static_cast<int>(s)) * P;
78 P = Scalar(RealScalar(1)/RealScalar(s + 1)) * P * Ashifted;
81 const RealScalar F_norm = F.cwiseAbs().rowwise().sum().maxCoeff();
82 const RealScalar Fincr_norm = Fincr.cwiseAbs().rowwise().sum().maxCoeff();
83 if (Fincr_norm < NumTraits<Scalar>::epsilon() * F_norm) {
85 RealScalar rfactorial = 1;
86 for (
Index r = 0; r < rows; r++) {
88 for (
Index i = 0; i < rows; i++)
89 mx = (std::max)(mx, std::abs(m_f(Ashifted(i, i) + avgEival,
static_cast<int>(s+r))));
91 rfactorial *= RealScalar(r);
92 delta = (std::max)(delta, mx / rfactorial);
94 const RealScalar P_norm = P.cwiseAbs().rowwise().sum().maxCoeff();
95 if (mu * delta * P_norm < NumTraits<Scalar>::epsilon() * F_norm)
107template <
typename Index,
typename ListOfClusters>
108typename ListOfClusters::iterator matrix_function_find_cluster(
Index key, ListOfClusters& clusters)
110 typename std::list<Index>::iterator j;
111 for (
typename ListOfClusters::iterator i = clusters.begin(); i != clusters.end(); ++i) {
112 j = std::find(i->begin(), i->end(), key);
116 return clusters.end();
130template <
typename EivalsType,
typename Cluster>
131void matrix_function_partition_eigenvalues(
const EivalsType& eivals, std::list<Cluster>& clusters)
133 typedef typename EivalsType::RealScalar RealScalar;
134 for (
Index i=0; i<eivals.rows(); ++i) {
136 typename std::list<Cluster>::iterator qi = matrix_function_find_cluster(i, clusters);
137 if (qi == clusters.end()) {
140 clusters.push_back(l);
146 for (
Index j=i+1; j<eivals.rows(); ++j) {
147 if (
abs(eivals(j) - eivals(i)) <= RealScalar(matrix_function_separation)
148 && std::find(qi->begin(), qi->end(), j) == qi->end()) {
149 typename std::list<Cluster>::iterator qj = matrix_function_find_cluster(j, clusters);
150 if (qj == clusters.end()) {
153 qi->insert(qi->end(), qj->begin(), qj->end());
162template <
typename ListOfClusters,
typename Index>
163void matrix_function_compute_cluster_size(
const ListOfClusters& clusters, Matrix<Index, Dynamic, 1>& clusterSize)
165 const Index numClusters =
static_cast<Index>(clusters.size());
166 clusterSize.setZero(numClusters);
167 Index clusterIndex = 0;
168 for (
typename ListOfClusters::const_iterator cluster = clusters.begin(); cluster != clusters.end(); ++cluster) {
169 clusterSize[clusterIndex] = cluster->size();
175template <
typename VectorType>
176void matrix_function_compute_block_start(
const VectorType& clusterSize, VectorType& blockStart)
178 blockStart.resize(clusterSize.rows());
180 for (
Index i = 1; i < clusterSize.rows(); i++) {
181 blockStart(i) = blockStart(i-1) + clusterSize(i-1);
186template <
typename EivalsType,
typename ListOfClusters,
typename VectorType>
187void matrix_function_compute_map(
const EivalsType& eivals,
const ListOfClusters& clusters, VectorType& eivalToCluster)
189 eivalToCluster.resize(eivals.rows());
190 Index clusterIndex = 0;
191 for (
typename ListOfClusters::const_iterator cluster = clusters.begin(); cluster != clusters.end(); ++cluster) {
192 for (
Index i = 0; i < eivals.rows(); ++i) {
193 if (std::find(cluster->begin(), cluster->end(), i) != cluster->end()) {
194 eivalToCluster[i] = clusterIndex;
202template <
typename DynVectorType,
typename VectorType>
203void matrix_function_compute_permutation(
const DynVectorType& blockStart,
const DynVectorType& eivalToCluster, VectorType& permutation)
205 DynVectorType indexNextEntry = blockStart;
206 permutation.resize(eivalToCluster.rows());
207 for (
Index i = 0; i < eivalToCluster.rows(); i++) {
208 Index cluster = eivalToCluster[i];
209 permutation[i] = indexNextEntry[cluster];
210 ++indexNextEntry[cluster];
215template <
typename VectorType,
typename MatrixType>
216void matrix_function_permute_schur(VectorType& permutation, MatrixType& U, MatrixType& T)
218 for (
Index i = 0; i < permutation.rows() - 1; i++) {
220 for (j = i; j < permutation.rows(); j++) {
221 if (permutation(j) == i)
break;
223 eigen_assert(permutation(j) == i);
224 for (
Index k = j-1; k >= i; k--) {
225 JacobiRotation<typename MatrixType::Scalar> rotation;
226 rotation.makeGivens(T(k, k+1), T(k+1, k+1) - T(k, k));
227 T.applyOnTheLeft(k, k+1, rotation.adjoint());
228 T.applyOnTheRight(k, k+1, rotation);
229 U.applyOnTheRight(k, k+1, rotation);
230 std::swap(permutation.coeffRef(k), permutation.coeffRef(k+1));
241template <
typename MatrixType,
typename AtomicType,
typename VectorType>
242void matrix_function_compute_block_atomic(
const MatrixType& T, AtomicType& atomic,
const VectorType& blockStart,
const VectorType& clusterSize, MatrixType& fT)
244 fT.setZero(T.rows(), T.cols());
245 for (
Index i = 0; i < clusterSize.rows(); ++i) {
246 fT.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i))
247 = atomic.compute(T.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i)));
273template <
typename MatrixType>
274MatrixType matrix_function_solve_triangular_sylvester(
const MatrixType& A,
const MatrixType& B,
const MatrixType& C)
276 eigen_assert(A.rows() == A.cols());
277 eigen_assert(A.isUpperTriangular());
278 eigen_assert(B.rows() == B.cols());
279 eigen_assert(B.isUpperTriangular());
280 eigen_assert(C.rows() == A.rows());
281 eigen_assert(C.cols() == B.rows());
283 typedef typename MatrixType::Scalar Scalar;
289 for (
Index i = m - 1; i >= 0; --i) {
290 for (
Index j = 0; j < n; ++j) {
297 Matrix<Scalar,1,1> AXmatrix = A.row(i).tail(m-1-i) * X.col(j).tail(m-1-i);
306 Matrix<Scalar,1,1> XBmatrix = X.row(i).head(j) * B.col(j).head(j);
310 X(i,j) = (C(i,j) - AX - XB) / (A(i,i) + B(j,j));
322template <
typename MatrixType,
typename VectorType>
323void matrix_function_compute_above_diagonal(
const MatrixType& T,
const VectorType& blockStart,
const VectorType& clusterSize, MatrixType& fT)
325 typedef internal::traits<MatrixType> Traits;
326 typedef typename MatrixType::Scalar Scalar;
327 static const int Options = MatrixType::Options;
328 typedef Matrix<Scalar, Dynamic, Dynamic, Options, Traits::RowsAtCompileTime, Traits::ColsAtCompileTime> DynMatrixType;
330 for (
Index k = 1; k < clusterSize.rows(); k++) {
331 for (
Index i = 0; i < clusterSize.rows() - k; i++) {
333 DynMatrixType A = T.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i));
334 DynMatrixType B = -T.block(blockStart(i+k), blockStart(i+k), clusterSize(i+k), clusterSize(i+k));
335 DynMatrixType C = fT.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i))
336 * T.block(blockStart(i), blockStart(i+k), clusterSize(i), clusterSize(i+k));
337 C -= T.block(blockStart(i), blockStart(i+k), clusterSize(i), clusterSize(i+k))
338 * fT.block(blockStart(i+k), blockStart(i+k), clusterSize(i+k), clusterSize(i+k));
339 for (
Index m = i + 1; m < i + k; m++) {
340 C += fT.block(blockStart(i), blockStart(m), clusterSize(i), clusterSize(m))
341 * T.block(blockStart(m), blockStart(i+k), clusterSize(m), clusterSize(i+k));
342 C -= T.block(blockStart(i), blockStart(m), clusterSize(i), clusterSize(m))
343 * fT.block(blockStart(m), blockStart(i+k), clusterSize(m), clusterSize(i+k));
345 fT.block(blockStart(i), blockStart(i+k), clusterSize(i), clusterSize(i+k))
346 = matrix_function_solve_triangular_sylvester(A, B, C);
366template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
367struct matrix_function_compute
379 template <
typename AtomicType,
typename ResultType>
380 static void run(
const MatrixType& A, AtomicType& atomic, ResultType &result);
389template <
typename MatrixType>
390struct matrix_function_compute<MatrixType, 0>
392 template <
typename MatA,
typename AtomicType,
typename ResultType>
393 static void run(
const MatA& A, AtomicType& atomic, ResultType &result)
395 typedef internal::traits<MatrixType> Traits;
396 typedef typename Traits::Scalar Scalar;
397 static const int Rows = Traits::RowsAtCompileTime, Cols = Traits::ColsAtCompileTime;
398 static const int MaxRows = Traits::MaxRowsAtCompileTime, MaxCols = Traits::MaxColsAtCompileTime;
400 typedef std::complex<Scalar> ComplexScalar;
401 typedef Matrix<ComplexScalar, Rows, Cols, 0, MaxRows, MaxCols> ComplexMatrix;
403 ComplexMatrix CA = A.template cast<ComplexScalar>();
404 ComplexMatrix Cresult;
405 matrix_function_compute<ComplexMatrix>::run(CA, atomic, Cresult);
406 result = Cresult.real();
413template <
typename MatrixType>
414struct matrix_function_compute<MatrixType, 1>
416 template <
typename MatA,
typename AtomicType,
typename ResultType>
417 static void run(
const MatA& A, AtomicType& atomic, ResultType &result)
419 typedef internal::traits<MatrixType> Traits;
422 const ComplexSchur<MatrixType> schurOfA(A);
423 eigen_assert(schurOfA.info()==
Success);
424 MatrixType T = schurOfA.matrixT();
425 MatrixType U = schurOfA.matrixU();
428 std::list<std::list<Index> > clusters;
429 matrix_function_partition_eigenvalues(T.diagonal(), clusters);
432 Matrix<Index, Dynamic, 1> clusterSize;
433 matrix_function_compute_cluster_size(clusters, clusterSize);
436 Matrix<Index, Dynamic, 1> blockStart;
437 matrix_function_compute_block_start(clusterSize, blockStart);
440 Matrix<Index, Dynamic, 1> eivalToCluster;
441 matrix_function_compute_map(T.diagonal(), clusters, eivalToCluster);
444 Matrix<Index, Traits::RowsAtCompileTime, 1> permutation;
445 matrix_function_compute_permutation(blockStart, eivalToCluster, permutation);
448 matrix_function_permute_schur(permutation, U, T);
452 matrix_function_compute_block_atomic(T, atomic, blockStart, clusterSize, fT);
453 matrix_function_compute_above_diagonal(T, blockStart, clusterSize, fT);
454 result = U * (fT.template triangularView<Upper>() * U.adjoint());
471:
public ReturnByValue<MatrixFunctionReturnValue<Derived> >
474 typedef typename Derived::Scalar Scalar;
475 typedef typename internal::stem_function<Scalar>::type StemFunction;
478 typedef typename internal::ref_selector<Derived>::type DerivedNested;
493 template <
typename ResultType>
494 inline void evalTo(ResultType& result)
const
496 typedef typename internal::nested_eval<Derived, 10>::type NestedEvalType;
497 typedef typename internal::remove_all<NestedEvalType>::type NestedEvalTypeClean;
498 typedef internal::traits<NestedEvalTypeClean> Traits;
499 typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
502 typedef internal::MatrixFunctionAtomic<DynMatrixType> AtomicType;
503 AtomicType atomic(m_f);
505 internal::matrix_function_compute<typename NestedEvalTypeClean::PlainObject>::run(m_A, atomic, result);
508 Index rows()
const {
return m_A.rows(); }
509 Index cols()
const {
return m_A.cols(); }
512 const DerivedNested m_A;
517template<
typename Derived>
518struct traits<MatrixFunctionReturnValue<Derived> >
520 typedef typename Derived::PlainObject ReturnType;
528template <
typename Derived>
531 eigen_assert(rows() == cols());
532 return MatrixFunctionReturnValue<Derived>(derived(), f);
535template <
typename Derived>
538 eigen_assert(rows() == cols());
539 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
540 return MatrixFunctionReturnValue<Derived>(derived(), internal::stem_function_sin<ComplexScalar>);
543template <
typename Derived>
546 eigen_assert(rows() == cols());
547 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
548 return MatrixFunctionReturnValue<Derived>(derived(), internal::stem_function_cos<ComplexScalar>);
551template <
typename Derived>
554 eigen_assert(rows() == cols());
555 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
556 return MatrixFunctionReturnValue<Derived>(derived(), internal::stem_function_sinh<ComplexScalar>);
559template <
typename Derived>
562 eigen_assert(rows() == cols());
563 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
564 return MatrixFunctionReturnValue<Derived>(derived(), internal::stem_function_cosh<ComplexScalar>);
const MatrixFunctionReturnValue< Derived > matrixFunction(StemFunction f) const
Definition: MatrixFunction.h:529
const MatrixFunctionReturnValue< Derived > cosh() const
Definition: MatrixFunction.h:560
const MatrixFunctionReturnValue< Derived > sin() const
Definition: MatrixFunction.h:536
const MatrixFunctionReturnValue< Derived > sinh() const
Definition: MatrixFunction.h:552
const MatrixFunctionReturnValue< Derived > cos() const
Definition: MatrixFunction.h:544
Proxy for the matrix function of some matrix (expression).
Definition: MatrixFunction.h:472
void evalTo(ResultType &result) const
Compute the matrix function.
Definition: MatrixFunction.h:494
MatrixFunctionReturnValue(const Derived &A, StemFunction f)
Constructor.
Definition: MatrixFunction.h:487
Namespace containing all symbols from the Eigen library.
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index