10 #ifndef EIGEN_MATRIX_FUNCTION 11 #define EIGEN_MATRIX_FUNCTION 13 #include "StemFunction.h" 21 static const float matrix_function_separation = 0.1f;
29 template <
typename MatrixType>
30 class 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);
52 template <
typename MatrixType>
53 typename NumTraits<typename MatrixType::Scalar>::Real matrix_function_compute_mu(
const MatrixType& A)
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();
63 template <
typename MatrixType>
64 MatrixType MatrixFunctionAtomic<MatrixType>::compute(
const MatrixType& A)
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;
76 for (Index s = 1; s < 1.1 * rows + 10; s++) {
77 Fincr = m_f(avgEival, static_cast<int>(s)) * P;
79 P = Scalar(RealScalar(1.0/(s + 1))) * P * Ashifted;
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) {
86 RealScalar rfactorial = 1;
87 for (Index r = 0; r < rows; r++) {
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))));
92 rfactorial *= RealScalar(r);
93 delta = (std::max)(delta, mx / rfactorial);
95 const RealScalar P_norm = P.cwiseAbs().rowwise().sum().maxCoeff();
96 if (mu * delta * P_norm < NumTraits<Scalar>::epsilon() * F_norm)
108 template <
typename Index,
typename ListOfClusters>
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);
117 return clusters.end();
131 template <
typename EivalsType,
typename Cluster>
134 typedef typename EivalsType::Index Index;
135 typedef typename EivalsType::RealScalar RealScalar;
136 for (Index i=0; i<eivals.rows(); ++i) {
139 if (qi == clusters.end()) {
142 clusters.push_back(l);
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()) {
152 if (qj == clusters.end()) {
155 qi->insert(qi->end(), qj->begin(), qj->end());
164 template <
typename ListOfClusters,
typename Index>
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();
177 template <
typename VectorType>
180 blockStart.resize(clusterSize.rows());
182 for (
typename VectorType::Index i = 1; i < clusterSize.rows(); i++) {
183 blockStart(i) = blockStart(i-1) + clusterSize(i-1);
188 template <
typename EivalsType,
typename ListOfClusters,
typename VectorType>
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;
205 template <
typename DynVectorType,
typename VectorType>
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];
219 template <
typename VectorType,
typename MatrixType>
222 typedef typename VectorType::Index Index;
223 for (Index i = 0; i < permutation.rows() - 1; i++) {
225 for (j = i; j < permutation.rows(); j++) {
226 if (permutation(j) == i)
break;
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));
246 template <
typename MatrixType,
typename AtomicType,
typename VectorType>
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)));
278 template <
typename MatrixType>
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());
288 typedef typename MatrixType::Index Index;
289 typedef typename MatrixType::Scalar Scalar;
295 for (Index i = m - 1; i >= 0; --i) {
296 for (Index j = 0; j < n; ++j) {
303 Matrix<Scalar,1,1> AXmatrix = A.row(i).tail(m-1-i) * X.col(j).tail(m-1-i);
312 Matrix<Scalar,1,1> XBmatrix = X.row(i).head(j) * B.col(j).head(j);
316 X(i,j) = (C(i,j) - AX - XB) / (A(i,i) + B(j,j));
328 template <
typename MatrixType,
typename VectorType>
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;
339 for (Index k = 1; k < clusterSize.rows(); k++) {
340 for (Index i = 0; i < clusterSize.rows() - k; i++) {
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));
354 fT.block(blockStart(i), blockStart(i+k), clusterSize(i), clusterSize(i+k))
375 template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
376 struct matrix_function_compute
388 template <
typename AtomicType,
typename ResultType>
389 static void run(
const MatrixType& A, AtomicType& atomic, ResultType &result);
398 template <
typename MatrixType>
399 struct matrix_function_compute<MatrixType, 0>
401 template <
typename MatA,
typename AtomicType,
typename ResultType>
402 static void run(
const MatA& A, AtomicType& atomic, ResultType &result)
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;
409 typedef std::complex<Scalar> ComplexScalar;
410 typedef Matrix<ComplexScalar, Rows, Cols, 0, MaxRows, MaxCols> ComplexMatrix;
412 ComplexMatrix CA = A.template cast<ComplexScalar>();
413 ComplexMatrix Cresult;
414 matrix_function_compute<ComplexMatrix>::run(CA, atomic, Cresult);
415 result = Cresult.real();
422 template <
typename MatrixType>
423 struct matrix_function_compute<MatrixType, 1>
425 template <
typename MatA,
typename AtomicType,
typename ResultType>
426 static void run(
const MatA& A, AtomicType& atomic, ResultType &result)
428 typedef internal::traits<MatrixType> Traits;
431 const ComplexSchur<MatrixType> schurOfA(A);
432 MatrixType T = schurOfA.matrixT();
433 MatrixType U = schurOfA.matrixU();
436 std::list<std::list<Index> > clusters;
440 Matrix<Index, Dynamic, 1> clusterSize;
444 Matrix<Index, Dynamic, 1> blockStart;
448 Matrix<Index, Dynamic, 1> eivalToCluster;
452 Matrix<Index, Traits::RowsAtCompileTime, 1> permutation;
462 result = U * (fT.template triangularView<Upper>() * U.adjoint());
479 :
public ReturnByValue<MatrixFunctionReturnValue<Derived> >
482 typedef typename Derived::Scalar Scalar;
483 typedef typename Derived::Index Index;
484 typedef typename internal::stem_function<Scalar>::type StemFunction;
487 typedef typename internal::ref_selector<Derived>::type DerivedNested;
502 template <
typename ResultType>
503 inline void evalTo(ResultType& result)
const 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;
513 typedef internal::MatrixFunctionAtomic<DynMatrixType> AtomicType;
514 AtomicType atomic(m_f);
516 internal::matrix_function_compute<typename NestedEvalTypeClean::PlainObject>::run(m_A, atomic, result);
519 Index rows()
const {
return m_A.rows(); }
520 Index cols()
const {
return m_A.cols(); }
523 const DerivedNested m_A;
528 template<
typename Derived>
529 struct traits<MatrixFunctionReturnValue<Derived> >
531 typedef typename Derived::PlainObject ReturnType;
539 template <
typename Derived>
540 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(
typename internal::stem_function<
typename internal::traits<Derived>::Scalar>::type f)
const 542 eigen_assert(rows() == cols());
546 template <
typename Derived>
549 eigen_assert(rows() == cols());
550 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
554 template <
typename Derived>
557 eigen_assert(rows() == cols());
558 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
562 template <
typename Derived>
565 eigen_assert(rows() == cols());
566 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
570 template <
typename Derived>
573 eigen_assert(rows() == cols());
574 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
580 #endif // EIGEN_MATRIX_FUNCTION void matrix_function_compute_permutation(const DynVectorType &blockStart, const DynVectorType &eivalToCluster, VectorType &permutation)
Compute permutation which groups ei'vals in same cluster together.
Definition: MatrixFunction.h:206
void matrix_function_compute_map(const EivalsType &eivals, const ListOfClusters &clusters, VectorType &eivalToCluster)
Compute mapping of eigenvalue indices to cluster indices.
Definition: MatrixFunction.h:189
Namespace containing all symbols from the Eigen library.
Definition: AdolcForward:45
Proxy for the matrix function of some matrix (expression).
Definition: MatrixFunction.h:478
void matrix_function_compute_block_start(const VectorType &clusterSize, VectorType &blockStart)
Compute start of each block using clusterSize.
Definition: MatrixFunction.h:178
MatrixFunctionReturnValue(const Derived &A, StemFunction f)
Constructor.
Definition: MatrixFunction.h:496
void matrix_function_compute_block_atomic(const MatrixType &T, AtomicType &atomic, const VectorType &blockStart, const VectorType &clusterSize, MatrixType &fT)
Compute block diagonal part of matrix function.
Definition: MatrixFunction.h:247
void matrix_function_partition_eigenvalues(const EivalsType &eivals, std::list< Cluster > &clusters)
Partition eigenvalues in clusters of ei'vals close to each other.
Definition: MatrixFunction.h:132
ListOfClusters::iterator matrix_function_find_cluster(Index key, ListOfClusters &clusters)
Find cluster in clusters containing some value.
Definition: MatrixFunction.h:109
void matrix_function_compute_cluster_size(const ListOfClusters &clusters, Matrix< Index, Dynamic, 1 > &clusterSize)
Compute size of each cluster given a partitioning.
Definition: MatrixFunction.h:165
void matrix_function_permute_schur(VectorType &permutation, MatrixType &U, MatrixType &T)
Permute Schur decomposition in U and T according to permutation.
Definition: MatrixFunction.h:220
void matrix_function_compute_above_diagonal(const MatrixType &T, const VectorType &blockStart, const VectorType &clusterSize, MatrixType &fT)
Compute part of matrix function above block diagonal.
Definition: MatrixFunction.h:329
void evalTo(ResultType &result) const
Compute the matrix function.
Definition: MatrixFunction.h:503
MatrixType matrix_function_solve_triangular_sylvester(const MatrixType &A, const MatrixType &B, const MatrixType &C)
Solve a triangular Sylvester equation AX + XB = C.
Definition: MatrixFunction.h:279