MatrixFunctionAtomic.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2009 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_ATOMIC
00026 #define EIGEN_MATRIX_FUNCTION_ATOMIC
00027 
00028 namespace Eigen { 
00029 
00038 template <typename MatrixType>
00039 class MatrixFunctionAtomic
00040 {
00041   public:
00042 
00043     typedef typename MatrixType::Scalar Scalar;
00044     typedef typename MatrixType::Index Index;
00045     typedef typename NumTraits<Scalar>::Real RealScalar;
00046     typedef typename internal::stem_function<Scalar>::type StemFunction;
00047     typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
00048 
00052     MatrixFunctionAtomic(StemFunction f) : m_f(f) { }
00053 
00058     MatrixType compute(const MatrixType& A);
00059 
00060   private:
00061 
00062     // Prevent copying
00063     MatrixFunctionAtomic(const MatrixFunctionAtomic&);
00064     MatrixFunctionAtomic& operator=(const MatrixFunctionAtomic&);
00065 
00066     void computeMu();
00067     bool taylorConverged(Index s, const MatrixType& F, const MatrixType& Fincr, const MatrixType& P);
00068 
00070     StemFunction* m_f;
00071 
00073     Index m_Arows;
00074 
00076     Scalar m_avgEival;
00077 
00079     MatrixType m_Ashifted;
00080 
00082     RealScalar m_mu;
00083 };
00084 
00085 template <typename MatrixType>
00086 MatrixType MatrixFunctionAtomic<MatrixType>::compute(const MatrixType& A)
00087 {
00088   // TODO: Use that A is upper triangular
00089   m_Arows = A.rows();
00090   m_avgEival = A.trace() / Scalar(RealScalar(m_Arows));
00091   m_Ashifted = A - m_avgEival * MatrixType::Identity(m_Arows, m_Arows);
00092   computeMu();
00093   MatrixType F = m_f(m_avgEival, 0) * MatrixType::Identity(m_Arows, m_Arows);
00094   MatrixType P = m_Ashifted;
00095   MatrixType Fincr;
00096   for (Index s = 1; s < 1.1 * m_Arows + 10; s++) { // upper limit is fairly arbitrary
00097     Fincr = m_f(m_avgEival, static_cast<int>(s)) * P;
00098     F += Fincr;
00099     P = Scalar(RealScalar(1.0/(s + 1))) * P * m_Ashifted;
00100     if (taylorConverged(s, F, Fincr, P)) {
00101       return F;
00102     }
00103   }
00104   eigen_assert("Taylor series does not converge" && 0);
00105   return F;
00106 }
00107 
00109 template <typename MatrixType>
00110 void MatrixFunctionAtomic<MatrixType>::computeMu()
00111 {
00112   const MatrixType N = MatrixType::Identity(m_Arows, m_Arows) - m_Ashifted;
00113   VectorType e = VectorType::Ones(m_Arows);
00114   N.template triangularView<Upper>().solveInPlace(e);
00115   m_mu = e.cwiseAbs().maxCoeff();
00116 }
00117 
00119 template <typename MatrixType>
00120 bool MatrixFunctionAtomic<MatrixType>::taylorConverged(Index s, const MatrixType& F,
00121                                                        const MatrixType& Fincr, const MatrixType& P)
00122 {
00123   const Index n = F.rows();
00124   const RealScalar F_norm = F.cwiseAbs().rowwise().sum().maxCoeff();
00125   const RealScalar Fincr_norm = Fincr.cwiseAbs().rowwise().sum().maxCoeff();
00126   if (Fincr_norm < NumTraits<Scalar>::epsilon() * F_norm) {
00127     RealScalar delta = 0;
00128     RealScalar rfactorial = 1;
00129     for (Index r = 0; r < n; r++) {
00130       RealScalar mx = 0;
00131       for (Index i = 0; i < n; i++)
00132         mx = (std::max)(mx, std::abs(m_f(m_Ashifted(i, i) + m_avgEival, static_cast<int>(s+r))));
00133       if (r != 0)
00134         rfactorial *= RealScalar(r);
00135       delta = (std::max)(delta, mx / rfactorial);
00136     }
00137     const RealScalar P_norm = P.cwiseAbs().rowwise().sum().maxCoeff();
00138     if (m_mu * delta * P_norm < NumTraits<Scalar>::epsilon() * F_norm)
00139       return true;
00140   }
00141   return false;
00142 }
00143 
00144 } // end namespace Eigen
00145 
00146 #endif // EIGEN_MATRIX_FUNCTION_ATOMIC