MatrixFunction.h
Go to the documentation of this file.
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


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:38:06