MatrixFunctionAtomic.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 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 
00036 template <typename MatrixType>
00037 class MatrixFunctionAtomic
00038 {
00039   public:
00040 
00041     typedef typename MatrixType::Scalar Scalar;
00042     typedef typename MatrixType::Index Index;
00043     typedef typename NumTraits<Scalar>::Real RealScalar;
00044     typedef typename internal::stem_function<Scalar>::type StemFunction;
00045     typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
00046 
00050     MatrixFunctionAtomic(StemFunction f) : m_f(f) { }
00051 
00056     MatrixType compute(const MatrixType& A);
00057 
00058   private:
00059 
00060     // Prevent copying
00061     MatrixFunctionAtomic(const MatrixFunctionAtomic&);
00062     MatrixFunctionAtomic& operator=(const MatrixFunctionAtomic&);
00063 
00064     void computeMu();
00065     bool taylorConverged(Index s, const MatrixType& F, const MatrixType& Fincr, const MatrixType& P);
00066 
00068     StemFunction* m_f;
00069 
00071     Index m_Arows;
00072 
00074     Scalar m_avgEival;
00075 
00077     MatrixType m_Ashifted;
00078 
00080     RealScalar m_mu;
00081 };
00082 
00083 template <typename MatrixType>
00084 MatrixType MatrixFunctionAtomic<MatrixType>::compute(const MatrixType& A)
00085 {
00086   // TODO: Use that A is upper triangular
00087   m_Arows = A.rows();
00088   m_avgEival = A.trace() / Scalar(RealScalar(m_Arows));
00089   m_Ashifted = A - m_avgEival * MatrixType::Identity(m_Arows, m_Arows);
00090   computeMu();
00091   MatrixType F = m_f(m_avgEival, 0) * MatrixType::Identity(m_Arows, m_Arows);
00092   MatrixType P = m_Ashifted;
00093   MatrixType Fincr;
00094   for (Index s = 1; s < 1.1 * m_Arows + 10; s++) { // upper limit is fairly arbitrary
00095     Fincr = m_f(m_avgEival, static_cast<int>(s)) * P;
00096     F += Fincr;
00097     P = Scalar(RealScalar(1.0/(s + 1))) * P * m_Ashifted;
00098     if (taylorConverged(s, F, Fincr, P)) {
00099       return F;
00100     }
00101   }
00102   eigen_assert("Taylor series does not converge" && 0);
00103   return F;
00104 }
00105 
00107 template <typename MatrixType>
00108 void MatrixFunctionAtomic<MatrixType>::computeMu()
00109 {
00110   const MatrixType N = MatrixType::Identity(m_Arows, m_Arows) - m_Ashifted;
00111   VectorType e = VectorType::Ones(m_Arows);
00112   N.template triangularView<Upper>().solveInPlace(e);
00113   m_mu = e.cwiseAbs().maxCoeff();
00114 }
00115 
00117 template <typename MatrixType>
00118 bool MatrixFunctionAtomic<MatrixType>::taylorConverged(Index s, const MatrixType& F,
00119                                                        const MatrixType& Fincr, const MatrixType& P)
00120 {
00121   const Index n = F.rows();
00122   const RealScalar F_norm = F.cwiseAbs().rowwise().sum().maxCoeff();
00123   const RealScalar Fincr_norm = Fincr.cwiseAbs().rowwise().sum().maxCoeff();
00124   if (Fincr_norm < NumTraits<Scalar>::epsilon() * F_norm) {
00125     RealScalar delta = 0;
00126     RealScalar rfactorial = 1;
00127     for (Index r = 0; r < n; r++) {
00128       RealScalar mx = 0;
00129       for (Index i = 0; i < n; i++)
00130         mx = (std::max)(mx, std::abs(m_f(m_Ashifted(i, i) + m_avgEival, static_cast<int>(s+r))));
00131       if (r != 0)
00132         rfactorial *= RealScalar(r);
00133       delta = (std::max)(delta, mx / rfactorial);
00134     }
00135     const RealScalar P_norm = P.cwiseAbs().rowwise().sum().maxCoeff();
00136     if (m_mu * delta * P_norm < NumTraits<Scalar>::epsilon() * F_norm)
00137       return true;
00138   }
00139   return false;
00140 }
00141 
00142 #endif // EIGEN_MATRIX_FUNCTION_ATOMIC


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:07