$treeview $search $mathjax
Eigen-unsupported
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2009-2011 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 #include "MatrixFunctionAtomic.h" 00015 00016 00017 namespace Eigen { 00018 00034 template <typename MatrixType, 00035 typename AtomicType, 00036 int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex> 00037 class MatrixFunction 00038 { 00039 public: 00040 00049 MatrixFunction(const MatrixType& A, AtomicType& atomic); 00050 00059 template <typename ResultType> 00060 void compute(ResultType &result); 00061 }; 00062 00063 00067 template <typename MatrixType, typename AtomicType> 00068 class MatrixFunction<MatrixType, AtomicType, 0> 00069 { 00070 private: 00071 00072 typedef internal::traits<MatrixType> Traits; 00073 typedef typename Traits::Scalar Scalar; 00074 static const int Rows = Traits::RowsAtCompileTime; 00075 static const int Cols = Traits::ColsAtCompileTime; 00076 static const int Options = MatrixType::Options; 00077 static const int MaxRows = Traits::MaxRowsAtCompileTime; 00078 static const int MaxCols = Traits::MaxColsAtCompileTime; 00079 00080 typedef std::complex<Scalar> ComplexScalar; 00081 typedef Matrix<ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols> ComplexMatrix; 00082 00083 public: 00084 00090 MatrixFunction(const MatrixType& A, AtomicType& atomic) : m_A(A), m_atomic(atomic) { } 00091 00101 template <typename ResultType> 00102 void compute(ResultType& result) 00103 { 00104 ComplexMatrix CA = m_A.template cast<ComplexScalar>(); 00105 ComplexMatrix Cresult; 00106 MatrixFunction<ComplexMatrix, AtomicType> mf(CA, m_atomic); 00107 mf.compute(Cresult); 00108 result = Cresult.real(); 00109 } 00110 00111 private: 00112 typename internal::nested<MatrixType>::type m_A; 00113 AtomicType& m_atomic; 00115 MatrixFunction& operator=(const MatrixFunction&); 00116 }; 00117 00118 00122 template <typename MatrixType, typename AtomicType> 00123 class MatrixFunction<MatrixType, AtomicType, 1> 00124 { 00125 private: 00126 00127 typedef internal::traits<MatrixType> Traits; 00128 typedef typename MatrixType::Scalar Scalar; 00129 typedef typename MatrixType::Index Index; 00130 static const int RowsAtCompileTime = Traits::RowsAtCompileTime; 00131 static const int ColsAtCompileTime = Traits::ColsAtCompileTime; 00132 static const int Options = MatrixType::Options; 00133 typedef typename NumTraits<Scalar>::Real RealScalar; 00134 typedef Matrix<Scalar, Traits::RowsAtCompileTime, 1> VectorType; 00135 typedef Matrix<Index, Traits::RowsAtCompileTime, 1> IntVectorType; 00136 typedef Matrix<Index, Dynamic, 1> DynamicIntVectorType; 00137 typedef std::list<Scalar> Cluster; 00138 typedef std::list<Cluster> ListOfClusters; 00139 typedef Matrix<Scalar, Dynamic, Dynamic, Options, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType; 00140 00141 public: 00142 00143 MatrixFunction(const MatrixType& A, AtomicType& atomic); 00144 template <typename ResultType> void compute(ResultType& result); 00145 00146 private: 00147 00148 void computeSchurDecomposition(); 00149 void partitionEigenvalues(); 00150 typename ListOfClusters::iterator findCluster(Scalar key); 00151 void computeClusterSize(); 00152 void computeBlockStart(); 00153 void constructPermutation(); 00154 void permuteSchur(); 00155 void swapEntriesInSchur(Index index); 00156 void computeBlockAtomic(); 00157 Block<MatrixType> block(MatrixType& A, Index i, Index j); 00158 void computeOffDiagonal(); 00159 DynMatrixType solveTriangularSylvester(const DynMatrixType& A, const DynMatrixType& B, const DynMatrixType& C); 00160 00161 typename internal::nested<MatrixType>::type m_A; 00162 AtomicType& m_atomic; 00163 MatrixType m_T; 00164 MatrixType m_U; 00165 MatrixType m_fT; 00166 ListOfClusters m_clusters; 00167 DynamicIntVectorType m_eivalToCluster; 00168 DynamicIntVectorType m_clusterSize; 00169 DynamicIntVectorType m_blockStart; 00170 IntVectorType m_permutation; 00178 static const RealScalar separation() { return static_cast<RealScalar>(0.1); } 00179 00180 MatrixFunction& operator=(const MatrixFunction&); 00181 }; 00182 00188 template <typename MatrixType, typename AtomicType> 00189 MatrixFunction<MatrixType,AtomicType,1>::MatrixFunction(const MatrixType& A, AtomicType& atomic) 00190 : m_A(A), m_atomic(atomic) 00191 { 00192 /* empty body */ 00193 } 00194 00200 template <typename MatrixType, typename AtomicType> 00201 template <typename ResultType> 00202 void MatrixFunction<MatrixType,AtomicType,1>::compute(ResultType& result) 00203 { 00204 computeSchurDecomposition(); 00205 partitionEigenvalues(); 00206 computeClusterSize(); 00207 computeBlockStart(); 00208 constructPermutation(); 00209 permuteSchur(); 00210 computeBlockAtomic(); 00211 computeOffDiagonal(); 00212 result = m_U * (m_fT.template triangularView<Upper>() * m_U.adjoint()); 00213 } 00214 00216 template <typename MatrixType, typename AtomicType> 00217 void MatrixFunction<MatrixType,AtomicType,1>::computeSchurDecomposition() 00218 { 00219 const ComplexSchur<MatrixType> schurOfA(m_A); 00220 m_T = schurOfA.matrixT(); 00221 m_U = schurOfA.matrixU(); 00222 } 00223 00235 template <typename MatrixType, typename AtomicType> 00236 void MatrixFunction<MatrixType,AtomicType,1>::partitionEigenvalues() 00237 { 00238 using std::abs; 00239 const Index rows = m_T.rows(); 00240 VectorType diag = m_T.diagonal(); // contains eigenvalues of A 00241 00242 for (Index i=0; i<rows; ++i) { 00243 // Find set containing diag(i), adding a new set if necessary 00244 typename ListOfClusters::iterator qi = findCluster(diag(i)); 00245 if (qi == m_clusters.end()) { 00246 Cluster l; 00247 l.push_back(diag(i)); 00248 m_clusters.push_back(l); 00249 qi = m_clusters.end(); 00250 --qi; 00251 } 00252 00253 // Look for other element to add to the set 00254 for (Index j=i+1; j<rows; ++j) { 00255 if (abs(diag(j) - diag(i)) <= separation() && std::find(qi->begin(), qi->end(), diag(j)) == qi->end()) { 00256 typename ListOfClusters::iterator qj = findCluster(diag(j)); 00257 if (qj == m_clusters.end()) { 00258 qi->push_back(diag(j)); 00259 } else { 00260 qi->insert(qi->end(), qj->begin(), qj->end()); 00261 m_clusters.erase(qj); 00262 } 00263 } 00264 } 00265 } 00266 } 00267 00273 template <typename MatrixType, typename AtomicType> 00274 typename MatrixFunction<MatrixType,AtomicType,1>::ListOfClusters::iterator MatrixFunction<MatrixType,AtomicType,1>::findCluster(Scalar key) 00275 { 00276 typename Cluster::iterator j; 00277 for (typename ListOfClusters::iterator i = m_clusters.begin(); i != m_clusters.end(); ++i) { 00278 j = std::find(i->begin(), i->end(), key); 00279 if (j != i->end()) 00280 return i; 00281 } 00282 return m_clusters.end(); 00283 } 00284 00286 template <typename MatrixType, typename AtomicType> 00287 void MatrixFunction<MatrixType,AtomicType,1>::computeClusterSize() 00288 { 00289 const Index rows = m_T.rows(); 00290 VectorType diag = m_T.diagonal(); 00291 const Index numClusters = static_cast<Index>(m_clusters.size()); 00292 00293 m_clusterSize.setZero(numClusters); 00294 m_eivalToCluster.resize(rows); 00295 Index clusterIndex = 0; 00296 for (typename ListOfClusters::const_iterator cluster = m_clusters.begin(); cluster != m_clusters.end(); ++cluster) { 00297 for (Index i = 0; i < diag.rows(); ++i) { 00298 if (std::find(cluster->begin(), cluster->end(), diag(i)) != cluster->end()) { 00299 ++m_clusterSize[clusterIndex]; 00300 m_eivalToCluster[i] = clusterIndex; 00301 } 00302 } 00303 ++clusterIndex; 00304 } 00305 } 00306 00308 template <typename MatrixType, typename AtomicType> 00309 void MatrixFunction<MatrixType,AtomicType,1>::computeBlockStart() 00310 { 00311 m_blockStart.resize(m_clusterSize.rows()); 00312 m_blockStart(0) = 0; 00313 for (Index i = 1; i < m_clusterSize.rows(); i++) { 00314 m_blockStart(i) = m_blockStart(i-1) + m_clusterSize(i-1); 00315 } 00316 } 00317 00319 template <typename MatrixType, typename AtomicType> 00320 void MatrixFunction<MatrixType,AtomicType,1>::constructPermutation() 00321 { 00322 DynamicIntVectorType indexNextEntry = m_blockStart; 00323 m_permutation.resize(m_T.rows()); 00324 for (Index i = 0; i < m_T.rows(); i++) { 00325 Index cluster = m_eivalToCluster[i]; 00326 m_permutation[i] = indexNextEntry[cluster]; 00327 ++indexNextEntry[cluster]; 00328 } 00329 } 00330 00332 template <typename MatrixType, typename AtomicType> 00333 void MatrixFunction<MatrixType,AtomicType,1>::permuteSchur() 00334 { 00335 IntVectorType p = m_permutation; 00336 for (Index i = 0; i < p.rows() - 1; i++) { 00337 Index j; 00338 for (j = i; j < p.rows(); j++) { 00339 if (p(j) == i) break; 00340 } 00341 eigen_assert(p(j) == i); 00342 for (Index k = j-1; k >= i; k--) { 00343 swapEntriesInSchur(k); 00344 std::swap(p.coeffRef(k), p.coeffRef(k+1)); 00345 } 00346 } 00347 } 00348 00350 template <typename MatrixType, typename AtomicType> 00351 void MatrixFunction<MatrixType,AtomicType,1>::swapEntriesInSchur(Index index) 00352 { 00353 JacobiRotation<Scalar> rotation; 00354 rotation.makeGivens(m_T(index, index+1), m_T(index+1, index+1) - m_T(index, index)); 00355 m_T.applyOnTheLeft(index, index+1, rotation.adjoint()); 00356 m_T.applyOnTheRight(index, index+1, rotation); 00357 m_U.applyOnTheRight(index, index+1, rotation); 00358 } 00359 00366 template <typename MatrixType, typename AtomicType> 00367 void MatrixFunction<MatrixType,AtomicType,1>::computeBlockAtomic() 00368 { 00369 m_fT.resize(m_T.rows(), m_T.cols()); 00370 m_fT.setZero(); 00371 for (Index i = 0; i < m_clusterSize.rows(); ++i) { 00372 block(m_fT, i, i) = m_atomic.compute(block(m_T, i, i)); 00373 } 00374 } 00375 00377 template <typename MatrixType, typename AtomicType> 00378 Block<MatrixType> MatrixFunction<MatrixType,AtomicType,1>::block(MatrixType& A, Index i, Index j) 00379 { 00380 return A.block(m_blockStart(i), m_blockStart(j), m_clusterSize(i), m_clusterSize(j)); 00381 } 00382 00390 template <typename MatrixType, typename AtomicType> 00391 void MatrixFunction<MatrixType,AtomicType,1>::computeOffDiagonal() 00392 { 00393 for (Index diagIndex = 1; diagIndex < m_clusterSize.rows(); diagIndex++) { 00394 for (Index blockIndex = 0; blockIndex < m_clusterSize.rows() - diagIndex; blockIndex++) { 00395 // compute (blockIndex, blockIndex+diagIndex) block 00396 DynMatrixType A = block(m_T, blockIndex, blockIndex); 00397 DynMatrixType B = -block(m_T, blockIndex+diagIndex, blockIndex+diagIndex); 00398 DynMatrixType C = block(m_fT, blockIndex, blockIndex) * block(m_T, blockIndex, blockIndex+diagIndex); 00399 C -= block(m_T, blockIndex, blockIndex+diagIndex) * block(m_fT, blockIndex+diagIndex, blockIndex+diagIndex); 00400 for (Index k = blockIndex + 1; k < blockIndex + diagIndex; k++) { 00401 C += block(m_fT, blockIndex, k) * block(m_T, k, blockIndex+diagIndex); 00402 C -= block(m_T, blockIndex, k) * block(m_fT, k, blockIndex+diagIndex); 00403 } 00404 block(m_fT, blockIndex, blockIndex+diagIndex) = solveTriangularSylvester(A, B, C); 00405 } 00406 } 00407 } 00408 00432 template <typename MatrixType, typename AtomicType> 00433 typename MatrixFunction<MatrixType,AtomicType,1>::DynMatrixType MatrixFunction<MatrixType,AtomicType,1>::solveTriangularSylvester( 00434 const DynMatrixType& A, 00435 const DynMatrixType& B, 00436 const DynMatrixType& C) 00437 { 00438 eigen_assert(A.rows() == A.cols()); 00439 eigen_assert(A.isUpperTriangular()); 00440 eigen_assert(B.rows() == B.cols()); 00441 eigen_assert(B.isUpperTriangular()); 00442 eigen_assert(C.rows() == A.rows()); 00443 eigen_assert(C.cols() == B.rows()); 00444 00445 Index m = A.rows(); 00446 Index n = B.rows(); 00447 DynMatrixType X(m, n); 00448 00449 for (Index i = m - 1; i >= 0; --i) { 00450 for (Index j = 0; j < n; ++j) { 00451 00452 // Compute AX = \sum_{k=i+1}^m A_{ik} X_{kj} 00453 Scalar AX; 00454 if (i == m - 1) { 00455 AX = 0; 00456 } else { 00457 Matrix<Scalar,1,1> AXmatrix = A.row(i).tail(m-1-i) * X.col(j).tail(m-1-i); 00458 AX = AXmatrix(0,0); 00459 } 00460 00461 // Compute XB = \sum_{k=1}^{j-1} X_{ik} B_{kj} 00462 Scalar XB; 00463 if (j == 0) { 00464 XB = 0; 00465 } else { 00466 Matrix<Scalar,1,1> XBmatrix = X.row(i).head(j) * B.col(j).head(j); 00467 XB = XBmatrix(0,0); 00468 } 00469 00470 X(i,j) = (C(i,j) - AX - XB) / (A(i,i) + B(j,j)); 00471 } 00472 } 00473 return X; 00474 } 00475 00488 template<typename Derived> class MatrixFunctionReturnValue 00489 : public ReturnByValue<MatrixFunctionReturnValue<Derived> > 00490 { 00491 public: 00492 00493 typedef typename Derived::Scalar Scalar; 00494 typedef typename Derived::Index Index; 00495 typedef typename internal::stem_function<Scalar>::type StemFunction; 00496 00503 MatrixFunctionReturnValue(const Derived& A, StemFunction f) : m_A(A), m_f(f) { } 00504 00510 template <typename ResultType> 00511 inline void evalTo(ResultType& result) const 00512 { 00513 typedef typename Derived::PlainObject PlainObject; 00514 typedef internal::traits<PlainObject> Traits; 00515 static const int RowsAtCompileTime = Traits::RowsAtCompileTime; 00516 static const int ColsAtCompileTime = Traits::ColsAtCompileTime; 00517 static const int Options = PlainObject::Options; 00518 typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar; 00519 typedef Matrix<ComplexScalar, Dynamic, Dynamic, Options, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType; 00520 typedef MatrixFunctionAtomic<DynMatrixType> AtomicType; 00521 AtomicType atomic(m_f); 00522 00523 const PlainObject Aevaluated = m_A.eval(); 00524 MatrixFunction<PlainObject, AtomicType> mf(Aevaluated, atomic); 00525 mf.compute(result); 00526 } 00527 00528 Index rows() const { return m_A.rows(); } 00529 Index cols() const { return m_A.cols(); } 00530 00531 private: 00532 typename internal::nested<Derived>::type m_A; 00533 StemFunction *m_f; 00534 00535 MatrixFunctionReturnValue& operator=(const MatrixFunctionReturnValue&); 00536 }; 00537 00538 namespace internal { 00539 template<typename Derived> 00540 struct traits<MatrixFunctionReturnValue<Derived> > 00541 { 00542 typedef typename Derived::PlainObject ReturnType; 00543 }; 00544 } 00545 00546 00547 /********** MatrixBase methods **********/ 00548 00549 00550 template <typename Derived> 00551 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename internal::stem_function<typename internal::traits<Derived>::Scalar>::type f) const 00552 { 00553 eigen_assert(rows() == cols()); 00554 return MatrixFunctionReturnValue<Derived>(derived(), f); 00555 } 00556 00557 template <typename Derived> 00558 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const 00559 { 00560 eigen_assert(rows() == cols()); 00561 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar; 00562 return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::sin); 00563 } 00564 00565 template <typename Derived> 00566 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const 00567 { 00568 eigen_assert(rows() == cols()); 00569 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar; 00570 return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::cos); 00571 } 00572 00573 template <typename Derived> 00574 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const 00575 { 00576 eigen_assert(rows() == cols()); 00577 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar; 00578 return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::sinh); 00579 } 00580 00581 template <typename Derived> 00582 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const 00583 { 00584 eigen_assert(rows() == cols()); 00585 typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar; 00586 return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::cosh); 00587 } 00588 00589 } // end namespace Eigen 00590 00591 #endif // EIGEN_MATRIX_FUNCTION