Remove useless code (abort specialization for complex exponent temporarily)

1 files changed, 0 insertions, 241 deletions

diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
index 69c4000cf..c4dd8ab29 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
@@ -231,89 +231,6 @@ class MatrixPower<MatrixType, IntExponent, PlainObject, 1>
     void partialPivLuSolve(IntExponent p, ResultType& result);
 };
 
-/**
- * \internal \ingroup MatrixFunctions_Module
- * \brief Partial specialization for complex matrices raised to complex exponents.
- */
-template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger>
-class MatrixPower<MatrixType, std::complex<RealScalar>, PlainObject, IsInteger>
-{
-  private:
-    typedef internal::traits<MatrixType> Traits;
-    static const int Rows = Traits::RowsAtCompileTime;
-    static const int Cols = Traits::ColsAtCompileTime;
-    static const int Options = Traits::Options;
-    static const int MaxRows = Traits::MaxRowsAtCompileTime;
-    static const int MaxCols = Traits::MaxColsAtCompileTime;
-
-    typedef typename MatrixType::Scalar Scalar;
-    typedef typename MatrixType::Index Index;
-    typedef Array<Scalar, Rows, 1, ColMajor, MaxRows> ArrayType;
-
-  public:
-    /**
-     * \brief Constructor.
-     *
-     * \param[in] A  the base of the matrix power.
-     * \param[in] p  the exponent of the matrix power.
-     * \param[in] b  the multiplier.
-     */
-    MatrixPower(const MatrixType& A, const Scalar& p, const PlainObject& b) :
-      m_A(A),
-      m_p(p),
-      m_b(b),
-      m_dimA(A.cols()),
-      m_dimb(b.cols())
-    { EIGEN_STATIC_ASSERT(false, COMPLEX_POWER_OF_A_MATRIX_IS_UNDER_CONSTRUCTION) }
-
-    /**
-     * \brief Compute the matrix power.
-     *
-     * \param[out] result  \f$ A^p b \f$, as specified in the constructor.
-     */
-    template <typename ResultType> void compute(ResultType& result);
- 
-  private:
-    /** \brief Compute Schur decomposition of #m_A. */
-    void computeSchurDecomposition();
-
-    /** \brief Compute atanh (inverse hyperbolic tangent). */
-    Scalar atanh(const Scalar& x);
-
-    /** \brief Compute power of 2x2 triangular matrix. */
-    void compute2x2(const Scalar& p);
-
-    /**
-     * \brief Compute power of triangular matrices with size > 2. 
-     * \details This uses a Schur-Pad&eacute; algorithm.
-     */
-    void computeBig();
-
-    /** \brief Get suitable degree for Pade approximation. (specialized for \c RealScalar = \c double) */
-    inline int getPadeDegree(double);
-/* TODO
- *  inline int getPadeDegree(float);
- *
- *  inline int getPadeDegree(long double);
- */
-    /** \brief Compute Pad&eacute; approximation to matrix fractional power. */
-    void computePade(int degree, const MatrixType& IminusT);
-
-    /** \brief Get a certain coefficient of the Pad&eacute; approximation. */
-    inline Scalar coeff(int degree);
-
-    const MatrixType& m_A;  ///< \brief Reference to the matrix base.
-    const Scalar& m_p;      ///< \brief Reference to the real exponent.
-    const PlainObject& m_b; ///< \brief Reference to the multiplier.
-    const Index m_dimA;     ///< \brief The dimension of #m_A, equivalent to %m_A.cols().
-    const Index m_dimb;     ///< \brief The dimension of #m_b, equivalent to %m_b.cols().
-    MatrixType m_tmp;       ///< \brief Used for temporary storage.
-    MatrixType m_T;         ///< \brief Triangular part of Schur decomposition.
-    MatrixType m_U;         ///< \brief Unitary part of Schur decomposition.
-    MatrixType m_fT;        ///< \brief #m_T to the power of #m_pfrac.
-    ArrayType m_logTdiag;   ///< \brief Logarithm of the main diagonal of #m_T.
-};
-
 /******* Specialized for real exponents *******/
 
 template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger>
@@ -642,164 +559,6 @@ void MatrixPower<MatrixType,IntExponent,PlainObject,1>::computeChainProduct(Resu
     result = m_tmp * result;
 }
 
-/******* Specialized for complex exponents *******/
-
-template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger>
-template <typename ResultType>
-void MatrixPower<MatrixType,std::complex<RealScalar>,PlainObject,IsInteger>::compute(ResultType& result)
-{
-  using std::floor;
-  using std::pow;
-
-  if (m_dimA == 1)
-    result = pow(m_A(0,0), m_p) * m_b;
-  else {
-    computeSchurDecomposition();
-    if (m_dimA == 2)
-      compute2x2(m_p);
-    else
-      computeBig();
-    result = m_U * m_fT * m_U.adjoint();
-  }
-}
-
-template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger>
-void MatrixPower<MatrixType,std::complex<RealScalar>,PlainObject,IsInteger>::computeSchurDecomposition()
-{
-  const ComplexSchur<MatrixType> schurOfA(m_A);
-  m_T = schurOfA.matrixT();
-  m_U = schurOfA.matrixU();
-  m_logTdiag = m_T.diagonal().array().log();
-}
-
-template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger>
-typename MatrixType::Scalar MatrixPower<MatrixType,std::complex<RealScalar>,PlainObject,IsInteger>::atanh(const Scalar& x)
-{
-  using std::abs;
-  using std::log;
-  using std::sqrt;
-
-  if (abs(x) > sqrt(NumTraits<RealScalar>::epsilon()))
-    return RealScalar(0.5) * log((RealScalar(1) + x) / (RealScalar(1) - x));
-  else
-    return x + x*x*x / RealScalar(3);
-}
-
-template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger>
-void MatrixPower<MatrixType,std::complex<RealScalar>,PlainObject,IsInteger>::compute2x2(const Scalar& p)
-{
-  using std::abs;
-  using std::ceil;
-  using std::exp;
-  using std::imag;
-  using std::ldexp;
-  using std::log;
-  using std::pow;
-  using std::sinh;
-
-  int i, j, unwindingNumber;
-  Scalar w;
-
-  m_fT(0,0) = pow(m_T(0,0), p);
-
-  for (j = 1; j < m_dimA; j++) {
-    i = j - 1;
-    m_fT(j,j) = pow(m_T(j,j), p);
-
-    if (m_T(i,i) == m_T(j,j))
-      m_fT(i,j) = p * pow(m_T(i,j), p - RealScalar(1));
-    else if (abs(m_T(i,i)) < ldexp(abs(m_T(j,j)), -1) || abs(m_T(j,j)) < ldexp(abs(m_T(i,i)), -1))
-      m_fT(i,j) = m_T(i,j) * (m_fT(j,j) - m_fT(i,i)) / (m_T(j,j) - m_T(i,i));
-    else {
-      // computation in previous branch is inaccurate if abs(m_T(j,j)) \approx abs(m_T(i,i))
-      unwindingNumber = static_cast<int>(ceil((imag(m_logTdiag[j] - m_logTdiag[i]) - M_PI) / (2 * M_PI)));
-      w = atanh((m_T(j,j) - m_T(i,i)) / (m_T(j,j) + m_T(i,i))) + Scalar(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));
-    }
-  }
-}
-
-template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger>
-void MatrixPower<MatrixType,std::complex<RealScalar>,PlainObject,IsInteger>::computeBig()
-{
-  using std::abs;
-  using std::ceil;
-  using std::frexp;
-  using std::ldexp;
-
-  const RealScalar maxNormForPade = 2.787629930862099e-1;
-  int degree, degree2, numberOfSquareRoots, numberOfExtraSquareRoots = 0;
-  MatrixType IminusT, sqrtT, T = m_T;
-  RealScalar normIminusT;
-  Scalar p;
-/*
-  frexp(abs(m_p), &numberOfSquareRoots);
-  if (numberOfSquareRoots > 0)
-    p = m_p * ldexp(RealScalar(1), -numberOfSquareRoots);
-  else {
-    p = m_p;
-    numberOfSquareRoots = 0;
-  }
-*/
-  while (true) {
-    IminusT = MatrixType::Identity(m_A.rows(), m_A.cols()) - T;
-    normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff();
-    if (normIminusT < maxNormForPade) {
-      degree = getPadeDegree(normIminusT);
-      degree2 = getPadeDegree(normIminusT * RealScalar(0.5));
-      if (degree - degree2 <= 1 || numberOfExtraSquareRoots)
-	break;
-      numberOfExtraSquareRoots++;
-    }
-    MatrixSquareRootTriangular<MatrixType>(T).compute(sqrtT);
-    T = sqrtT;
-    numberOfSquareRoots++;
-  }
-  computePade(degree, IminusT);
-
-  for (; numberOfSquareRoots; numberOfSquareRoots--) {
-    compute2x2(p * ldexp(RealScalar(1), -numberOfSquareRoots));
-    m_fT *= m_fT;
-  }
-  compute2x2(p);
-}
-
-template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger>
-inline int MatrixPower<MatrixType,std::complex<RealScalar>,PlainObject,IsInteger>::getPadeDegree(double normIminusT)
-{
-  const double maxNormForPade[] = { 1.882832775783885e-2 /* degree = 3 */ , 6.036100693089764e-2,
-      1.239372725584911e-1, 1.998030690604271e-1, 2.787629930862099e-1 };
-  for (int degree = 3; degree <= 7; degree++)
-    if (normIminusT <= maxNormForPade[degree - 3])
-      return degree;
-  assert(false); // this line should never be reached
-}
-
-template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger>
-void MatrixPower<MatrixType,std::complex<RealScalar>,PlainObject,IsInteger>::computePade(int degree, const MatrixType& IminusT)
-{
-  degree <<= 1;
-  m_fT = coeff(degree) * IminusT;
-
-  for (int i = degree - 1; i; i--) {
-    m_fT = (MatrixType::Identity(m_A.rows(), m_A.cols()) + m_fT).template triangularView<Upper>()
-	.solve(coeff(i) * IminusT).eval();
-  }
-  m_fT += MatrixType::Identity(m_A.rows(), m_A.cols());
-}
-
-template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger>
-inline typename MatrixType::Scalar MatrixPower<MatrixType,std::complex<RealScalar>,PlainObject,IsInteger>::coeff(int i)
-{
-  if (i == 1)
-    return -m_p;
-  else if (i & 1)
-    return (-m_p - RealScalar(i)) / RealScalar((i<<2) + 2);
-  else
-    return (m_p - RealScalar(i)) / RealScalar((i<<2) - 2);
-}
-
 /**
  * \ingroup MatrixFunctions_Module
  *