![]() |
Eigen-unsupported
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2009-2011, 2013 Jitse Niesen <jitse@maths.leeds.ac.uk> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_MATRIX_FUNCTION 00011 #define EIGEN_MATRIX_FUNCTION 00012 00013 #include "StemFunction.h" 00014 00015 00016 namespace Eigen { 00017 00018 namespace internal { 00019 00021 static const float matrix_function_separation = 0.1f; 00022 00029 template <typename MatrixType> 00030 class MatrixFunctionAtomic 00031 { 00032 public: 00033 00034 typedef typename MatrixType::Scalar Scalar; 00035 typedef typename stem_function<Scalar>::type StemFunction; 00036 00040 MatrixFunctionAtomic(StemFunction f) : m_f(f) { } 00041 00046 MatrixType compute(const MatrixType& A); 00047 00048 private: 00049 StemFunction* m_f; 00050 }; 00051 00052 template <typename MatrixType> 00053 typename NumTraits<typename MatrixType::Scalar>::Real matrix_function_compute_mu(const MatrixType& A) 00054 { 00055 typedef typename plain_col_type<MatrixType>::type VectorType; 00056 typename MatrixType::Index rows = A.rows(); 00057 const MatrixType N = MatrixType::Identity(rows, rows) - A; 00058 VectorType e = VectorType::Ones(rows); 00059 N.template triangularView<Upper>().solveInPlace(e); 00060 return e.cwiseAbs().maxCoeff(); 00061 } 00062 00063 template <typename MatrixType> 00064 MatrixType MatrixFunctionAtomic<MatrixType>::compute(const MatrixType& A) 00065 { 00066 // TODO: Use that A is upper triangular 00067 typedef typename NumTraits<Scalar>::Real RealScalar; 00068 typedef typename MatrixType::Index Index; 00069 Index rows = A.rows(); 00070 Scalar avgEival = A.trace() / Scalar(RealScalar(rows)); 00071 MatrixType Ashifted = A - avgEival * MatrixType::Identity(rows, rows); 00072 RealScalar mu = matrix_function_compute_mu(Ashifted); 00073 MatrixType F = m_f(avgEival, 0) * MatrixType::Identity(rows, rows); 00074 MatrixType P = Ashifted; 00075 MatrixType Fincr; 00076 for (Index s = 1; s < 1.1 * rows + 10; s++) { // upper limit is fairly arbitrary 00077 Fincr = m_f(avgEival, static_cast<int>(s)) * P; 00078 F += Fincr; 00079 P = Scalar(RealScalar(1.0/(s + 1))) * P * Ashifted; 00080 00081 // test whether Taylor series converged 00082 const RealScalar F_norm = F.cwiseAbs().rowwise().sum().maxCoeff(); 00083 const RealScalar Fincr_norm = Fincr.cwiseAbs().rowwise().sum().maxCoeff(); 00084 if (Fincr_norm < NumTraits<Scalar>::epsilon() * F_norm) { 00085 RealScalar delta = 0; 00086 RealScalar rfactorial = 1; 00087 for (Index r = 0; r < rows; r++) { 00088 RealScalar mx = 0; 00089 for (Index i = 0; i < rows; i++) 00090 mx = (std::max)(mx, std::abs(m_f(Ashifted(i, i) + avgEival, static_cast<int>(s+r)))); 00091 if (r != 0) 00092 rfactorial *= RealScalar(r); 00093 delta = (std::max)(delta, mx / rfactorial); 00094 } 00095 const RealScalar P_norm = P.cwiseAbs().rowwise().sum().maxCoeff(); 00096 if (mu * delta * P_norm < NumTraits<Scalar>::epsilon() * F_norm) // series converged 00097 break; 00098 } 00099 } 00100 return F; 00101 } 00102 00108 template <typename Index, typename ListOfClusters> 00109 typename ListOfClusters::iterator matrix_function_find_cluster(Index key, ListOfClusters& clusters) 00110 { 00111 typename std::list<Index>::iterator j; 00112 for (typename ListOfClusters::iterator i = clusters.begin(); i != clusters.end(); ++i) { 00113 j = std::find(i->begin(), i->end(), key); 00114 if (j != i->end()) 00115 return i; 00116 } 00117 return clusters.end(); 00118 } 00119 00131 template <typename EivalsType, typename Cluster> 00132 void matrix_function_partition_eigenvalues(const EivalsType& eivals, std::list<Cluster>& clusters) 00133 { 00134 typedef typename EivalsType::Index Index; 00135 typedef typename EivalsType::RealScalar RealScalar; 00136 for (Index i=0; i<eivals.rows(); ++i) { 00137 // Find cluster containing i-th ei'val, adding a new cluster if necessary 00138 typename std::list<Cluster>::iterator qi = matrix_function_find_cluster(i, clusters); 00139 if (qi == clusters.end()) { 00140 Cluster l; 00141 l.push_back(i); 00142 clusters.push_back(l); 00143 qi = clusters.end(); 00144 --qi; 00145 } 00146 00147 // Look for other element to add to the set 00148 for (Index j=i+1; j<eivals.rows(); ++j) { 00149 if (abs(eivals(j) - eivals(i)) <= RealScalar(matrix_function_separation) 00150 && std::find(qi->begin(), qi->end(), j) == qi->end()) { 00151 typename std::list<Cluster>::iterator qj = matrix_function_find_cluster(j, clusters); 00152 if (qj == clusters.end()) { 00153 qi->push_back(j); 00154 } else { 00155 qi->insert(qi->end(), qj->begin(), qj->end()); 00156 clusters.erase(qj); 00157 } 00158 } 00159 } 00160 } 00161 } 00162 00164 template <typename ListOfClusters, typename Index> 00165 void matrix_function_compute_cluster_size(const ListOfClusters& clusters, Matrix<Index, Dynamic, 1>& clusterSize) 00166 { 00167 const Index numClusters = static_cast<Index>(clusters.size()); 00168 clusterSize.setZero(numClusters); 00169 Index clusterIndex = 0; 00170 for (typename ListOfClusters::const_iterator cluster = clusters.begin(); cluster != clusters.end(); ++cluster) { 00171 clusterSize[clusterIndex] = cluster->size(); 00172 ++clusterIndex; 00173 } 00174 } 00175 00177 template <typename VectorType> 00178 void matrix_function_compute_block_start(const VectorType& clusterSize, VectorType& blockStart) 00179 { 00180 blockStart.resize(clusterSize.rows()); 00181 blockStart(0) = 0; 00182 for (typename VectorType::Index i = 1; i < clusterSize.rows(); i++) { 00183 blockStart(i) = blockStart(i-1) + clusterSize(i-1); 00184 } 00185 } 00186 00188 template <typename EivalsType, typename ListOfClusters, typename VectorType> 00189 void matrix_function_compute_map(const EivalsType& eivals, const ListOfClusters& clusters, VectorType& eivalToCluster) 00190 { 00191 typedef typename EivalsType::Index Index; 00192 eivalToCluster.resize(eivals.rows()); 00193 Index clusterIndex = 0; 00194 for (typename ListOfClusters::const_iterator cluster = clusters.begin(); cluster != clusters.end(); ++cluster) { 00195 for (Index i = 0; i < eivals.rows(); ++i) { 00196 if (std::find(cluster->begin(), cluster->end(), i) != cluster->end()) { 00197 eivalToCluster[i] = clusterIndex; 00198 } 00199 } 00200 ++clusterIndex; 00201 } 00202 } 00203 00205 template <typename DynVectorType, typename VectorType> 00206 void matrix_function_compute_permutation(const DynVectorType& blockStart, const DynVectorType& eivalToCluster, VectorType& permutation) 00207 { 00208 typedef typename VectorType::Index Index; 00209 DynVectorType indexNextEntry = blockStart; 00210 permutation.resize(eivalToCluster.rows()); 00211 for (Index i = 0; i < eivalToCluster.rows(); i++) { 00212 Index cluster = eivalToCluster[i]; 00213 permutation[i] = indexNextEntry[cluster]; 00214 ++indexNextEntry[cluster]; 00215 } 00216 } 00217 00219 template <typename VectorType, typename MatrixType> 00220 void matrix_function_permute_schur(VectorType& permutation, MatrixType& U, MatrixType& T) 00221 { 00222 typedef typename VectorType::Index Index; 00223 for (Index i = 0; i < permutation.rows() - 1; i++) { 00224 Index j; 00225 for (j = i; j < permutation.rows(); j++) { 00226 if (permutation(j) == i) break; 00227 } 00228 eigen_assert(permutation(j) == i); 00229 for (Index k = j-1; k >= i; k--) { 00230 JacobiRotation<typename MatrixType::Scalar> rotation; 00231 rotation.makeGivens(T(k, k+1), T(k+1, k+1) - T(k, k)); 00232 T.applyOnTheLeft(k, k+1, rotation.adjoint()); 00233 T.applyOnTheRight(k, k+1, rotation); 00234 U.applyOnTheRight(k, k+1, rotation); 00235 std::swap(permutation.coeffRef(k), permutation.coeffRef(k+1)); 00236 } 00237 } 00238 } 00239 00246 template <typename MatrixType, typename AtomicType, typename VectorType> 00247 void matrix_function_compute_block_atomic(const MatrixType& T, AtomicType& atomic, const VectorType& blockStart, const VectorType& clusterSize, MatrixType& fT) 00248 { 00249 fT.setZero(T.rows(), T.cols()); 00250 for (typename VectorType::Index i = 0; i < clusterSize.rows(); ++i) { 00251 fT.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i)) 00252 = atomic.compute(T.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i))); 00253 } 00254 } 00255 00278 template <typename MatrixType> 00279 MatrixType matrix_function_solve_triangular_sylvester(const MatrixType& A, const MatrixType& B, const MatrixType& C) 00280 { 00281 eigen_assert(A.rows() == A.cols()); 00282 eigen_assert(A.isUpperTriangular()); 00283 eigen_assert(B.rows() == B.cols()); 00284 eigen_assert(B.isUpperTriangular()); 00285 eigen_assert(C.rows() == A.rows()); 00286 eigen_assert(C.cols() == B.rows()); 00287 00288 typedef typename MatrixType::Index Index; 00289 typedef typename MatrixType::Scalar Scalar; 00290 00291 Index m = A.rows(); 00292 Index n = B.rows(); 00293 MatrixType X(m, n); 00294 00295 for (Index i = m - 1; i >= 0; --i) { 00296 for (Index j = 0; j < n; ++j) { 00297 00298 // Compute AX = \sum_{k=i+1}^m A_{ik} X_{kj} 00299 Scalar AX; 00300 if (i == m - 1) { 00301 AX = 0; 00302 } else { 00303 Matrix<Scalar,1,1> AXmatrix = A.row(i).tail(m-1-i) * X.col(j).tail(m-1-i); 00304 AX = AXmatrix(0,0); 00305 } 00306 00307 // Compute XB = \sum_{k=1}^{j-1} X_{ik} B_{kj} 00308 Scalar XB; 00309 if (j == 0) { 00310 XB = 0; 00311 } else { 00312 Matrix<Scalar,1,1> XBmatrix = X.row(i).head(j) * B.col(j).head(j); 00313 XB = XBmatrix(0,0); 00314 } 00315 00316 X(i,j) = (C(i,j) - AX - XB) / (A(i,i) + B(j,j)); 00317 } 00318 } 00319 return X; 00320 } 00321 00328 template <typename MatrixType, typename VectorType> 00329 void matrix_function_compute_above_diagonal(const MatrixType& T, const VectorType& blockStart, const VectorType& clusterSize, MatrixType& fT) 00330 { 00331 typedef internal::traits<MatrixType> Traits; 00332 typedef typename MatrixType::Scalar Scalar; 00333 typedef typename MatrixType::Index Index; 00334 static const int RowsAtCompileTime = Traits::RowsAtCompileTime; 00335 static const int ColsAtCompileTime = Traits::ColsAtCompileTime; 00336 static const int Options = MatrixType::Options; 00337 typedef Matrix<Scalar, Dynamic, Dynamic, Options, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType; 00338 00339 for (Index k = 1; k < clusterSize.rows(); k++) { 00340 for (Index i = 0; i < clusterSize.rows() - k; i++) { 00341 // compute (i, i+k) block 00342 DynMatrixType A = T.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i)); 00343 DynMatrixType B = -T.block(blockStart(i+k), blockStart(i+k), clusterSize(i+k), clusterSize(i+k)); 00344 DynMatrixType C = fT.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i)) 00345 * T.block(blockStart(i), blockStart(i+k), clusterSize(i), clusterSize(i+k)); 00346 C -= T.block(blockStart(i), blockStart(i+k), clusterSize(i), clusterSize(i+k)) 00347 * fT.block(blockStart(i+k), blockStart(i+k), clusterSize(i+k), clusterSize(i+k)); 00348 for (Index m = i + 1; m < i + k; m++) { 00349 C += fT.block(blockStart(i), blockStart(m), clusterSize(i), clusterSize(m)) 00350 * T.block(blockStart(m), blockStart(i+k), clusterSize(m), clusterSize(i+k)); 00351 C -= T.block(blockStart(i), blockStart(m), clusterSize(i), clusterSize(m)) 00352 * fT.block(blockStart(m), blockStart(i+k), clusterSize(m), clusterSize(i+k)); 00353 } 00354 fT.block(blockStart(i), blockStart(i+k), clusterSize(i), clusterSize(i+k)) 00355 = matrix_function_solve_triangular_sylvester(A, B, C); 00356 } 00357 } 00358 } 00359 00375 template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex> 00376 struct matrix_function_compute 00377 { 00388 template <typename AtomicType, typename ResultType> 00389 static void run(const MatrixType& A, AtomicType& atomic, ResultType &result); 00390 }; 00391 00398 template <typename MatrixType> 00399 struct matrix_function_compute<MatrixType, 0> 00400 { 00401 template <typename AtomicType, typename ResultType> 00402 static void run(const MatrixType& A, AtomicType& atomic, ResultType &result) 00403 { 00404 typedef internal::traits<MatrixType> Traits; 00405 typedef typename Traits::Scalar Scalar; 00406 static const int Rows = Traits::RowsAtCompileTime, Cols = Traits::ColsAtCompileTime; 00407 static const int MaxRows = Traits::MaxRowsAtCompileTime, MaxCols = Traits::MaxColsAtCompileTime; 00408 00409 typedef std::complex<Scalar> ComplexScalar; 00410 typedef Matrix<ComplexScalar, Rows, Cols, 0, MaxRows, MaxCols> ComplexMatrix; 00411 00412 ComplexMatrix CA = A.template cast<ComplexScalar>(); 00413 ComplexMatrix Cresult; 00414 matrix_function_compute<ComplexMatrix>::run(CA, atomic, Cresult); 00415 result = Cresult.real(); 00416 } 00417 }; 00418 00422 template <typename MatrixType> 00423 struct matrix_function_compute<MatrixType, 1> 00424 { 00425 template <typename AtomicType, typename ResultType> 00426 static void run(const MatrixType& A, AtomicType& atomic, ResultType &result) 00427 { 00428 typedef internal::traits<MatrixType> Traits; 00429 typedef typename MatrixType::Index Index; 00430 00431 // compute Schur decomposition of A 00432 const ComplexSchur<MatrixType> schurOfA(A); 00433 MatrixType T = schurOfA.matrixT(); 00434 MatrixType U = schurOfA.matrixU(); 00435 00436 // partition eigenvalues into clusters of ei'vals "close" to each other 00437 std::list<std::list<Index> > clusters; 00438 matrix_function_partition_eigenvalues(T.diagonal(), clusters); 00439 00440 // compute size of each cluster 00441 Matrix<Index, Dynamic, 1> clusterSize; 00442 matrix_function_compute_cluster_size(clusters, clusterSize); 00443 00444 // blockStart[i] is row index at which block corresponding to i-th cluster starts 00445 Matrix<Index, Dynamic, 1> blockStart; 00446 matrix_function_compute_block_start(clusterSize, blockStart); 00447 00448 // compute map so that eivalToCluster[i] = j means that i-th ei'val is in j-th cluster 00449 Matrix<Index, Dynamic, 1> eivalToCluster; 00450 matrix_function_compute_map(T.diagonal(), clusters, eivalToCluster); 00451 00452 // compute permutation which groups ei'vals in same cluster together 00453 Matrix<Index, Traits::RowsAtCompileTime, 1> permutation; 00454 matrix_function_compute_permutation(blockStart, eivalToCluster, permutation); 00455 00456 // permute Schur decomposition 00457 matrix_function_permute_schur(permutation, U, T); 00458 00459 // compute result 00460 MatrixType fT; // matrix function applied to T 00461 matrix_function_compute_block_atomic(T, atomic, blockStart, clusterSize, fT); 00462 matrix_function_compute_above_diagonal(T, blockStart, clusterSize, fT); 00463 result = U * (fT.template triangularView<Upper>() * U.adjoint()); 00464 } 00465 }; 00466 00467 } // end of namespace internal 00468 00479 template<typename Derived> class MatrixFunctionReturnValue 00480 : public ReturnByValue<MatrixFunctionReturnValue<Derived> > 00481 { 00482 public: 00483 typedef typename Derived::Scalar Scalar; 00484 typedef typename Derived::Index Index; 00485 typedef typename internal::stem_function<Scalar>::type StemFunction; 00486 00487 protected: 00488 typedef typename internal::ref_selector<Derived>::type DerivedNested; 00489 00490 public: 00491 00497 MatrixFunctionReturnValue(const Derived& A, StemFunction f) : m_A(A), m_f(f) { } 00498 00503 template <typename ResultType> 00504 inline void evalTo(ResultType& result) const 00505 { 00506 typedef typename internal::nested_eval<Derived, 10>::type NestedEvalType; 00507 typedef typename internal::remove_all<NestedEvalType>::type NestedEvalTypeClean; 00508 typedef internal::traits<NestedEvalTypeClean> Traits; 00509 static const int RowsAtCompileTime = Traits::RowsAtCompileTime; 00510 static const int ColsAtCompileTime = Traits::ColsAtCompileTime; 00511 typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar; 00512 typedef Matrix<ComplexScalar, Dynamic, Dynamic, 0, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType; 00513 00514 typedef internal::MatrixFunctionAtomic<DynMatrixType> AtomicType; 00515 AtomicType atomic(m_f); 00516 00517 internal::matrix_function_compute<NestedEvalTypeClean>::run(m_A, atomic, result); 00518 } 00519 00520 Index rows() const { return m_A.rows(); } 00521 Index cols() const { return m_A.cols(); } 00522 00523 private: 00524 const DerivedNested m_A; 00525 StemFunction *m_f; 00526 }; 00527 00528 namespace internal { 00529 template<typename Derived> 00530 struct traits<MatrixFunctionReturnValue<Derived> > 00531 { 00532 typedef typename Derived::PlainObject ReturnType; 00533 }; 00534 } 00535 00536 00537 /********** MatrixBase methods **********/ 00538 00539 00540 template <typename Derived> 00541 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename internal::stem_function<typename internal::traits<Derived>::Scalar>::type f) const 00542 { 00543 eigen_assert(rows() == cols()); 00544 return MatrixFunctionReturnValue<Derived>(derived(), f); 00545 } 00546 00547 template <typename Derived> 00548 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const 00549 { 00550 eigen_assert(rows() == cols()); 00551 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar; 00552 return MatrixFunctionReturnValue<Derived>(derived(), internal::stem_function_sin<ComplexScalar>); 00553 } 00554 00555 template <typename Derived> 00556 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const 00557 { 00558 eigen_assert(rows() == cols()); 00559 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar; 00560 return MatrixFunctionReturnValue<Derived>(derived(), internal::stem_function_cos<ComplexScalar>); 00561 } 00562 00563 template <typename Derived> 00564 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const 00565 { 00566 eigen_assert(rows() == cols()); 00567 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar; 00568 return MatrixFunctionReturnValue<Derived>(derived(), internal::stem_function_sinh<ComplexScalar>); 00569 } 00570 00571 template <typename Derived> 00572 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const 00573 { 00574 eigen_assert(rows() == cols()); 00575 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar; 00576 return MatrixFunctionReturnValue<Derived>(derived(), internal::stem_function_cosh<ComplexScalar>); 00577 } 00578 00579 } // end namespace Eigen 00580 00581 #endif // EIGEN_MATRIX_FUNCTION