GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-05-06 10:44:05 +08:00
Code Issues Packages Projects Releases Wiki Activity

Remove useless code (abort specialization for complex exponent temporarily)

Browse Source
This commit is contained in:
jdh8 2012-08-15 20:56:03 +08:00
parent 4be172d84f
commit 2337ea430b
1 changed files with 0 additions and 241 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

241
unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
View File

@ -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
*