Avoid inefficient 2x2 LU. Move atanh to internal for maintainability.

Eigen/src/Core/MathFunctions.h

@@ -519,6 +519,53 @@ inline EIGEN_MATHFUNC_RETVAL(atan2, Scalar) atan2(const Scalar& x, const Scalar&
   return EIGEN_MATHFUNC_IMPL(atan2, Scalar)::run(x, y);
 }
 
+/****************************************************************************
+* Implementation of atanh2                                                *
+****************************************************************************/
+
+template<typename Scalar, bool IsInteger>
+struct atanh2_default_impl
+{
+  typedef Scalar retval;
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+  static inline Scalar run(const Scalar& x, const Scalar& y)
+  {
+    using std::abs;
+    using std::log;
+    using std::sqrt;
+    Scalar z = x / y;
+    if (abs(z) > sqrt(NumTraits<RealScalar>::epsilon()))
+      return RealScalar(0.5) * log((y + x) / (y - x));
+    else
+      return z + z*z*z / RealScalar(3);
+  }
+};
+
+template<typename Scalar>
+struct atanh2_default_impl<Scalar, true>
+{
+  static inline Scalar run(const Scalar&, const Scalar&)
+  {
+    EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
+    return Scalar(0);
+  }
+};
+
+template<typename Scalar>
+struct atanh2_impl : atanh2_default_impl<Scalar, NumTraits<Scalar>::IsInteger> {};
+
+template<typename Scalar>
+struct atanh2_retval
+{
+  typedef Scalar type;
+};
+
+template<typename Scalar>
+inline EIGEN_MATHFUNC_RETVAL(atanh2, Scalar) atanh2(const Scalar& x, const Scalar& y)
+{
+  return EIGEN_MATHFUNC_IMPL(atanh2, Scalar)::run(x, y);
+}
+
 /****************************************************************************
 * Implementation of pow                                                  *
 ****************************************************************************/

doc/I02_HiPerformance.dox

@@ -79,7 +79,7 @@ temp = m2 * m3;
 m1 += temp.adjoint(); \endcode</td>
 <td>\code
 m1.noalias() += m3.adjoint()
-*             * m2.adjoint(); \endcode</td>
+*              * m2.adjoint(); \endcode</td>
 <td>This is because the product expression has the EvalBeforeNesting bit which
     enforces the evaluation of the product by the Tranpose expression.</td>
 </tr>

unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h

@@ -51,7 +51,6 @@ private:
 
   void compute2x2(const MatrixType& A, MatrixType& result);
   void computeBig(const MatrixType& A, MatrixType& result);
-  static Scalar atanh2(Scalar y, Scalar x);
   int getPadeDegree(float normTminusI);
   int getPadeDegree(double normTminusI);
   int getPadeDegree(long double normTminusI);
@@ -93,20 +92,6 @@ MatrixType MatrixLogarithmAtomic<MatrixType>::compute(const MatrixType& A)
   return result;
 }
 
-/** \brief Compute atanh (inverse hyperbolic tangent) for \f$ y / x \f$. */
-template <typename MatrixType>
-typename MatrixType::Scalar MatrixLogarithmAtomic<MatrixType>::atanh2(Scalar y, Scalar x)
-{
-  using std::abs;
-  using std::sqrt;
-
-  Scalar z = y / x;
-  if (abs(z) > sqrt(NumTraits<Scalar>::epsilon()))
-    return Scalar(0.5) * log((x + y) / (x - y));
-  else
-    return z + z*z*z / Scalar(3);
-}
-
 /** \brief Compute logarithm of 2x2 triangular matrix. */
 template <typename MatrixType>
 void MatrixLogarithmAtomic<MatrixType>::compute2x2(const MatrixType& A, MatrixType& result)
@@ -131,7 +116,7 @@ void MatrixLogarithmAtomic<MatrixType>::compute2x2(const MatrixType& A, MatrixTy
     // computation in previous branch is inaccurate if A(1,1) \approx A(0,0)
     int unwindingNumber = static_cast<int>(ceil((imag(logA11 - logA00) - M_PI) / (2*M_PI)));
     Scalar y = A(1,1) - A(0,0), x = A(1,1) + A(0,0);
-    result(0,1) = A(0,1) * (Scalar(2) * atanh2(y,x) + Scalar(0,2*M_PI*unwindingNumber)) / y;
+    result(0,1) = A(0,1) * (Scalar(2) * internal::atanh2(y,x) + Scalar(0,2*M_PI*unwindingNumber)) / y;
   }
 }
 

unsupported/Eigen/src/MatrixFunctions/MatrixPower.h

@@ -111,9 +111,6 @@ class MatrixPower
      */
     void getFractionalExponent();
 
-    /** \brief Compute atanh (inverse hyperbolic tangent) for \f$ y / x \f$. */
-    static ComplexScalar atanh2(const ComplexScalar& y, const ComplexScalar& x);
-
     /** \brief Compute power of 2x2 triangular matrix. */
     void compute2x2(RealScalar p);
 
@@ -223,7 +220,7 @@ void MatrixPower<MatrixType,PlainObject>::computeChainProduct(ResultType& result
   int cost = computeCost(p);
 
   if (m_pInt < RealScalar(0)) {
-    if (p * m_dimb <= cost * m_dimA) {
+    if (p * m_dimb <= cost * m_dimA && m_dimA > 2) {
       partialPivLuSolve(result, p);
       return;
     } else {
@@ -296,21 +293,6 @@ void MatrixPower<MatrixType,PlainObject>::getFractionalExponent()
   }
 }
 
-template<typename MatrixType, typename PlainObject>
-std::complex<typename MatrixType::RealScalar>
-MatrixPower<MatrixType,PlainObject>::atanh2(const ComplexScalar& y, const ComplexScalar& x)
-{
-  using std::abs;
-  using std::log;
-  using std::sqrt;
-  const ComplexScalar z = y / x;
-
-  if (abs(z) > sqrt(NumTraits<RealScalar>::epsilon()))
-    return RealScalar(0.5) * log((x + y) / (x - y));
-  else
-    return z + z*z*z / RealScalar(3);
-}
-
 template<typename MatrixType, typename PlainObject>
 void MatrixPower<MatrixType,PlainObject>::compute2x2(RealScalar p)
 {
@@ -337,7 +319,7 @@ void MatrixPower<MatrixType,PlainObject>::compute2x2(RealScalar p)
     } else {
       // computation in previous branch is inaccurate if abs(m_T(j,j)) \approx abs(m_T(i,i))
       unwindingNumber = ceil((imag(m_logTdiag[j] - m_logTdiag[i]) - M_PI) / (2 * M_PI));
-      w = atanh2(m_T(j,j) - m_T(i,i), m_T(j,j) + m_T(i,i)) + ComplexScalar(0, M_PI * unwindingNumber);
+      w = internal::atanh2(m_T(j,j) - m_T(i,i), m_T(j,j) + m_T(i,i)) + ComplexScalar(0, M_PI * unwindingNumber);
       m_fT(i,j) = m_T(i,j) * RealScalar(2) * exp(RealScalar(0.5) * p * (m_logTdiag[j] + m_logTdiag[i])) *
 	  sinh(p * w) / (m_T(j,j) - m_T(i,i));
     }