MatrixFunction.h
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 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_MATRIX_FUNCTION
00026 #define EIGEN_MATRIX_FUNCTION
00027 
00028 #include "StemFunction.h"
00029 #include "MatrixFunctionAtomic.h"
00030 
00031 
00032 namespace Eigen { 
00033 
00049 template <typename MatrixType, 
00050           typename AtomicType,  
00051           int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
00052 class MatrixFunction
00053 {  
00054   public:
00055 
00064     MatrixFunction(const MatrixType& A, AtomicType& atomic);
00065 
00074     template <typename ResultType> 
00075     void compute(ResultType &result);    
00076 };
00077 
00078 
00082 template <typename MatrixType, typename AtomicType>
00083 class MatrixFunction<MatrixType, AtomicType, 0>
00084 {  
00085   private:
00086 
00087     typedef internal::traits<MatrixType> Traits;
00088     typedef typename Traits::Scalar Scalar;
00089     static const int Rows = Traits::RowsAtCompileTime;
00090     static const int Cols = Traits::ColsAtCompileTime;
00091     static const int Options = MatrixType::Options;
00092     static const int MaxRows = Traits::MaxRowsAtCompileTime;
00093     static const int MaxCols = Traits::MaxColsAtCompileTime;
00094 
00095     typedef std::complex<Scalar> ComplexScalar;
00096     typedef Matrix<ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols> ComplexMatrix;
00097 
00098   public:
00099 
00105     MatrixFunction(const MatrixType& A, AtomicType& atomic) : m_A(A), m_atomic(atomic) { }
00106 
00116     template <typename ResultType>
00117     void compute(ResultType& result) 
00118     {
00119       ComplexMatrix CA = m_A.template cast<ComplexScalar>();
00120       ComplexMatrix Cresult;
00121       MatrixFunction<ComplexMatrix, AtomicType> mf(CA, m_atomic);
00122       mf.compute(Cresult);
00123       result = Cresult.real();
00124     }
00125 
00126   private:
00127     typename internal::nested<MatrixType>::type m_A; 
00128     AtomicType& m_atomic; 
00130     MatrixFunction& operator=(const MatrixFunction&);
00131 };
00132 
00133       
00137 template <typename MatrixType, typename AtomicType>
00138 class MatrixFunction<MatrixType, AtomicType, 1>
00139 {
00140   private:
00141 
00142     typedef internal::traits<MatrixType> Traits;
00143     typedef typename MatrixType::Scalar Scalar;
00144     typedef typename MatrixType::Index Index;
00145     static const int RowsAtCompileTime = Traits::RowsAtCompileTime;
00146     static const int ColsAtCompileTime = Traits::ColsAtCompileTime;
00147     static const int Options = MatrixType::Options;
00148     typedef typename NumTraits<Scalar>::Real RealScalar;
00149     typedef Matrix<Scalar, Traits::RowsAtCompileTime, 1> VectorType;
00150     typedef Matrix<Index, Traits::RowsAtCompileTime, 1> IntVectorType;
00151     typedef Matrix<Index, Dynamic, 1> DynamicIntVectorType;
00152     typedef std::list<Scalar> Cluster;
00153     typedef std::list<Cluster> ListOfClusters;
00154     typedef Matrix<Scalar, Dynamic, Dynamic, Options, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType;
00155 
00156   public:
00157 
00158     MatrixFunction(const MatrixType& A, AtomicType& atomic);
00159     template <typename ResultType> void compute(ResultType& result);
00160 
00161   private:
00162 
00163     void computeSchurDecomposition();
00164     void partitionEigenvalues();
00165     typename ListOfClusters::iterator findCluster(Scalar key);
00166     void computeClusterSize();
00167     void computeBlockStart();
00168     void constructPermutation();
00169     void permuteSchur();
00170     void swapEntriesInSchur(Index index);
00171     void computeBlockAtomic();
00172     Block<MatrixType> block(MatrixType& A, Index i, Index j);
00173     void computeOffDiagonal();
00174     DynMatrixType solveTriangularSylvester(const DynMatrixType& A, const DynMatrixType& B, const DynMatrixType& C);
00175 
00176     typename internal::nested<MatrixType>::type m_A; 
00177     AtomicType& m_atomic; 
00178     MatrixType m_T; 
00179     MatrixType m_U; 
00180     MatrixType m_fT; 
00181     ListOfClusters m_clusters; 
00182     DynamicIntVectorType m_eivalToCluster; 
00183     DynamicIntVectorType m_clusterSize; 
00184     DynamicIntVectorType m_blockStart; 
00185     IntVectorType m_permutation; 
00193     static const RealScalar separation() { return static_cast<RealScalar>(0.1); }
00194 
00195     MatrixFunction& operator=(const MatrixFunction&);
00196 };
00197 
00203 template <typename MatrixType, typename AtomicType>
00204 MatrixFunction<MatrixType,AtomicType,1>::MatrixFunction(const MatrixType& A, AtomicType& atomic)
00205   : m_A(A), m_atomic(atomic)
00206 {
00207   /* empty body */
00208 }
00209 
00215 template <typename MatrixType, typename AtomicType>
00216 template <typename ResultType>
00217 void MatrixFunction<MatrixType,AtomicType,1>::compute(ResultType& result) 
00218 {
00219   computeSchurDecomposition();
00220   partitionEigenvalues();
00221   computeClusterSize();
00222   computeBlockStart();
00223   constructPermutation();
00224   permuteSchur();
00225   computeBlockAtomic();
00226   computeOffDiagonal();
00227   result = m_U * m_fT * m_U.adjoint();
00228 }
00229 
00231 template <typename MatrixType, typename AtomicType>
00232 void MatrixFunction<MatrixType,AtomicType,1>::computeSchurDecomposition()
00233 {
00234   const ComplexSchur<MatrixType> schurOfA(m_A);  
00235   m_T = schurOfA.matrixT();
00236   m_U = schurOfA.matrixU();
00237 }
00238 
00250 template <typename MatrixType, typename AtomicType>
00251 void MatrixFunction<MatrixType,AtomicType,1>::partitionEigenvalues()
00252 {
00253   const Index rows = m_T.rows();
00254   VectorType diag = m_T.diagonal(); // contains eigenvalues of A
00255 
00256   for (Index i=0; i<rows; ++i) {
00257     // Find set containing diag(i), adding a new set if necessary
00258     typename ListOfClusters::iterator qi = findCluster(diag(i));
00259     if (qi == m_clusters.end()) {
00260       Cluster l;
00261       l.push_back(diag(i));
00262       m_clusters.push_back(l);
00263       qi = m_clusters.end();
00264       --qi;
00265     }
00266 
00267     // Look for other element to add to the set
00268     for (Index j=i+1; j<rows; ++j) {
00269       if (internal::abs(diag(j) - diag(i)) <= separation() && std::find(qi->begin(), qi->end(), diag(j)) == qi->end()) {
00270         typename ListOfClusters::iterator qj = findCluster(diag(j));
00271         if (qj == m_clusters.end()) {
00272           qi->push_back(diag(j));
00273         } else {
00274           qi->insert(qi->end(), qj->begin(), qj->end());
00275           m_clusters.erase(qj);
00276         }
00277       }
00278     }
00279   }
00280 }
00281 
00287 template <typename MatrixType, typename AtomicType>
00288 typename MatrixFunction<MatrixType,AtomicType,1>::ListOfClusters::iterator MatrixFunction<MatrixType,AtomicType,1>::findCluster(Scalar key)
00289 {
00290   typename Cluster::iterator j;
00291   for (typename ListOfClusters::iterator i = m_clusters.begin(); i != m_clusters.end(); ++i) {
00292     j = std::find(i->begin(), i->end(), key);
00293     if (j != i->end())
00294       return i;
00295   }
00296   return m_clusters.end();
00297 }
00298 
00300 template <typename MatrixType, typename AtomicType>
00301 void MatrixFunction<MatrixType,AtomicType,1>::computeClusterSize()
00302 {
00303   const Index rows = m_T.rows();
00304   VectorType diag = m_T.diagonal(); 
00305   const Index numClusters = static_cast<Index>(m_clusters.size());
00306 
00307   m_clusterSize.setZero(numClusters);
00308   m_eivalToCluster.resize(rows);
00309   Index clusterIndex = 0;
00310   for (typename ListOfClusters::const_iterator cluster = m_clusters.begin(); cluster != m_clusters.end(); ++cluster) {
00311     for (Index i = 0; i < diag.rows(); ++i) {
00312       if (std::find(cluster->begin(), cluster->end(), diag(i)) != cluster->end()) {
00313         ++m_clusterSize[clusterIndex];
00314         m_eivalToCluster[i] = clusterIndex;
00315       }
00316     }
00317     ++clusterIndex;
00318   }
00319 }
00320 
00322 template <typename MatrixType, typename AtomicType>
00323 void MatrixFunction<MatrixType,AtomicType,1>::computeBlockStart()
00324 {
00325   m_blockStart.resize(m_clusterSize.rows());
00326   m_blockStart(0) = 0;
00327   for (Index i = 1; i < m_clusterSize.rows(); i++) {
00328     m_blockStart(i) = m_blockStart(i-1) + m_clusterSize(i-1);
00329   }
00330 }
00331 
00333 template <typename MatrixType, typename AtomicType>
00334 void MatrixFunction<MatrixType,AtomicType,1>::constructPermutation()
00335 {
00336   DynamicIntVectorType indexNextEntry = m_blockStart;
00337   m_permutation.resize(m_T.rows());
00338   for (Index i = 0; i < m_T.rows(); i++) {
00339     Index cluster = m_eivalToCluster[i];
00340     m_permutation[i] = indexNextEntry[cluster];
00341     ++indexNextEntry[cluster];
00342   }
00343 }  
00344 
00346 template <typename MatrixType, typename AtomicType>
00347 void MatrixFunction<MatrixType,AtomicType,1>::permuteSchur()
00348 {
00349   IntVectorType p = m_permutation;
00350   for (Index i = 0; i < p.rows() - 1; i++) {
00351     Index j;
00352     for (j = i; j < p.rows(); j++) {
00353       if (p(j) == i) break;
00354     }
00355     eigen_assert(p(j) == i);
00356     for (Index k = j-1; k >= i; k--) {
00357       swapEntriesInSchur(k);
00358       std::swap(p.coeffRef(k), p.coeffRef(k+1));
00359     }
00360   }
00361 }
00362 
00364 template <typename MatrixType, typename AtomicType>
00365 void MatrixFunction<MatrixType,AtomicType,1>::swapEntriesInSchur(Index index)
00366 {
00367   JacobiRotation<Scalar> rotation;
00368   rotation.makeGivens(m_T(index, index+1), m_T(index+1, index+1) - m_T(index, index));
00369   m_T.applyOnTheLeft(index, index+1, rotation.adjoint());
00370   m_T.applyOnTheRight(index, index+1, rotation);
00371   m_U.applyOnTheRight(index, index+1, rotation);
00372 }  
00373 
00380 template <typename MatrixType, typename AtomicType>
00381 void MatrixFunction<MatrixType,AtomicType,1>::computeBlockAtomic()
00382 { 
00383   m_fT.resize(m_T.rows(), m_T.cols());
00384   m_fT.setZero();
00385   for (Index i = 0; i < m_clusterSize.rows(); ++i) {
00386     block(m_fT, i, i) = m_atomic.compute(block(m_T, i, i));
00387   }
00388 }
00389 
00391 template <typename MatrixType, typename AtomicType>
00392 Block<MatrixType> MatrixFunction<MatrixType,AtomicType,1>::block(MatrixType& A, Index i, Index j)
00393 {
00394   return A.block(m_blockStart(i), m_blockStart(j), m_clusterSize(i), m_clusterSize(j));
00395 }
00396 
00404 template <typename MatrixType, typename AtomicType>
00405 void MatrixFunction<MatrixType,AtomicType,1>::computeOffDiagonal()
00406 { 
00407   for (Index diagIndex = 1; diagIndex < m_clusterSize.rows(); diagIndex++) {
00408     for (Index blockIndex = 0; blockIndex < m_clusterSize.rows() - diagIndex; blockIndex++) {
00409       // compute (blockIndex, blockIndex+diagIndex) block
00410       DynMatrixType A = block(m_T, blockIndex, blockIndex);
00411       DynMatrixType B = -block(m_T, blockIndex+diagIndex, blockIndex+diagIndex);
00412       DynMatrixType C = block(m_fT, blockIndex, blockIndex) * block(m_T, blockIndex, blockIndex+diagIndex);
00413       C -= block(m_T, blockIndex, blockIndex+diagIndex) * block(m_fT, blockIndex+diagIndex, blockIndex+diagIndex);
00414       for (Index k = blockIndex + 1; k < blockIndex + diagIndex; k++) {
00415         C += block(m_fT, blockIndex, k) * block(m_T, k, blockIndex+diagIndex);
00416         C -= block(m_T, blockIndex, k) * block(m_fT, k, blockIndex+diagIndex);
00417       }
00418       block(m_fT, blockIndex, blockIndex+diagIndex) = solveTriangularSylvester(A, B, C);
00419     }
00420   }
00421 }
00422 
00446 template <typename MatrixType, typename AtomicType>
00447 typename MatrixFunction<MatrixType,AtomicType,1>::DynMatrixType MatrixFunction<MatrixType,AtomicType,1>::solveTriangularSylvester(
00448   const DynMatrixType& A, 
00449   const DynMatrixType& B, 
00450   const DynMatrixType& C)
00451 {
00452   eigen_assert(A.rows() == A.cols());
00453   eigen_assert(A.isUpperTriangular());
00454   eigen_assert(B.rows() == B.cols());
00455   eigen_assert(B.isUpperTriangular());
00456   eigen_assert(C.rows() == A.rows());
00457   eigen_assert(C.cols() == B.rows());
00458 
00459   Index m = A.rows();
00460   Index n = B.rows();
00461   DynMatrixType X(m, n);
00462 
00463   for (Index i = m - 1; i >= 0; --i) {
00464     for (Index j = 0; j < n; ++j) {
00465 
00466       // Compute AX = \sum_{k=i+1}^m A_{ik} X_{kj}
00467       Scalar AX;
00468       if (i == m - 1) {
00469         AX = 0; 
00470       } else {
00471         Matrix<Scalar,1,1> AXmatrix = A.row(i).tail(m-1-i) * X.col(j).tail(m-1-i);
00472         AX = AXmatrix(0,0);
00473       }
00474 
00475       // Compute XB = \sum_{k=1}^{j-1} X_{ik} B_{kj}
00476       Scalar XB;
00477       if (j == 0) {
00478         XB = 0; 
00479       } else {
00480         Matrix<Scalar,1,1> XBmatrix = X.row(i).head(j) * B.col(j).head(j);
00481         XB = XBmatrix(0,0);
00482       }
00483 
00484       X(i,j) = (C(i,j) - AX - XB) / (A(i,i) + B(j,j));
00485     }
00486   }
00487   return X;
00488 }
00489 
00502 template<typename Derived> class MatrixFunctionReturnValue
00503 : public ReturnByValue<MatrixFunctionReturnValue<Derived> >
00504 {
00505   public:
00506 
00507     typedef typename Derived::Scalar Scalar;
00508     typedef typename Derived::Index Index;
00509     typedef typename internal::stem_function<Scalar>::type StemFunction;
00510 
00517     MatrixFunctionReturnValue(const Derived& A, StemFunction f) : m_A(A), m_f(f) { }
00518 
00524     template <typename ResultType>
00525     inline void evalTo(ResultType& result) const
00526     {
00527       typedef typename Derived::PlainObject PlainObject;
00528       typedef internal::traits<PlainObject> Traits;
00529       static const int RowsAtCompileTime = Traits::RowsAtCompileTime;
00530       static const int ColsAtCompileTime = Traits::ColsAtCompileTime;
00531       static const int Options = PlainObject::Options;
00532       typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
00533       typedef Matrix<ComplexScalar, Dynamic, Dynamic, Options, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType;
00534       typedef MatrixFunctionAtomic<DynMatrixType> AtomicType;
00535       AtomicType atomic(m_f);
00536 
00537       const PlainObject Aevaluated = m_A.eval();
00538       MatrixFunction<PlainObject, AtomicType> mf(Aevaluated, atomic);
00539       mf.compute(result);
00540     }
00541 
00542     Index rows() const { return m_A.rows(); }
00543     Index cols() const { return m_A.cols(); }
00544 
00545   private:
00546     typename internal::nested<Derived>::type m_A;
00547     StemFunction *m_f;
00548 
00549     MatrixFunctionReturnValue& operator=(const MatrixFunctionReturnValue&);
00550 };
00551 
00552 namespace internal {
00553 template<typename Derived>
00554 struct traits<MatrixFunctionReturnValue<Derived> >
00555 {
00556   typedef typename Derived::PlainObject ReturnType;
00557 };
00558 }
00559 
00560 
00561 /********** MatrixBase methods **********/
00562 
00563 
00564 template <typename Derived>
00565 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename internal::stem_function<typename internal::traits<Derived>::Scalar>::type f) const
00566 {
00567   eigen_assert(rows() == cols());
00568   return MatrixFunctionReturnValue<Derived>(derived(), f);
00569 }
00570 
00571 template <typename Derived>
00572 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const
00573 {
00574   eigen_assert(rows() == cols());
00575   typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
00576   return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::sin);
00577 }
00578 
00579 template <typename Derived>
00580 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const
00581 {
00582   eigen_assert(rows() == cols());
00583   typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
00584   return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::cos);
00585 }
00586 
00587 template <typename Derived>
00588 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const
00589 {
00590   eigen_assert(rows() == cols());
00591   typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
00592   return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::sinh);
00593 }
00594 
00595 template <typename Derived>
00596 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const
00597 {
00598   eigen_assert(rows() == cols());
00599   typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
00600   return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::cosh);
00601 }
00602 
00603 } // end namespace Eigen
00604 
00605 #endif // EIGEN_MATRIX_FUNCTION