00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
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
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
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++) {
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 }
00145
00146 #endif // EIGEN_MATRIX_FUNCTION_ATOMIC