mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-11-21 03:11:25 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add test for real MatrixPowerTriangular.

Browse Source
This commit is contained in:
Chen-Pang He 2012-09-30 19:21:53 +08:00
parent eb33d307af
commit e92fe88159
3 changed files with 88 additions and 43 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -20,7 +20,7 @@ namespace Eigen {
* \tparam MatrixType type of the base, expected to be an instantiation
* of the Matrix class template.
*
* This class is capable of computing complex upper triangular matrices raised
* This class is capable of computing upper triangular matrices raised
* to an arbitrary real power.
*/
template<typename MatrixType>
@ -37,7 +37,7 @@ class MatrixPowerTriangular : public MatrixPowerBase<MatrixPowerTriangular<Matri
* The class stores a reference to A, so it should not be changed
* (or destroyed) before evaluation.
*/
explicit MatrixPowerTriangular(const MatrixType& A) : Base(A,0), m_T(Base::m_A)
explicit MatrixPowerTriangular(const MatrixType& A) : Base(A), m_T(Base::m_A)
{ }
#ifdef EIGEN_PARSED_BY_DOXYGEN
@ -262,7 +262,7 @@ class MatrixPower : public MatrixPowerBase<MatrixPower<MatrixType>,MatrixType>
* The class stores a reference to A, so it should not be changed
* (or destroyed) before evaluation.
*/
explicit MatrixPower(const MatrixType& A) : Base(A,0)
explicit MatrixPower(const MatrixType& A) : Base(A)
{ }
#ifdef EIGEN_PARSED_BY_DOXYGEN

69
unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
View File

@ -75,10 +75,10 @@ class MatrixPowerBase
typedef typename MatrixType::RealScalar RealScalar;
typedef typename MatrixType::Index Index;
explicit MatrixPowerBase(const MatrixType& A, RealScalar cond) :
explicit MatrixPowerBase(const MatrixType& A) :
m_A(A),
m_Id(MatrixType::Identity(A.rows(),A.cols())),
m_conditionNumber(cond)
m_Id(MatrixType::Identity(A.rows(), A.cols())),
m_conditionNumber(0)
{ eigen_assert(A.rows() == A.cols()); }
#ifndef EIGEN_PARSED_BY_DOXYGEN
@ -173,10 +173,8 @@ struct traits<MatrixPowerProduct<Derived,_Lhs,_Rhs> >
};
};
template<bool IsComplex> struct recompose_complex_schur;
template<>
struct recompose_complex_schur<true>
template<int IsComplex>
struct recompose_complex_schur
{
template<typename ResultType, typename MatrixType>
static inline void run(ResultType& res, const MatrixType& T, const MatrixType& U)
@ -184,14 +182,14 @@ struct recompose_complex_schur<true>
};
template<>
struct recompose_complex_schur<false>
struct recompose_complex_schur<0>
{
template<typename ResultType, typename MatrixType>
static inline void run(ResultType& res, const MatrixType& T, const MatrixType& U)
{ res.noalias() = (U * (T.template triangularView<Upper>() * U.adjoint())).real(); }
};
template<typename Scalar, int IsComplex=NumTraits<Scalar>::IsComplex>
template<typename Scalar, int IsComplex = NumTraits<Scalar>::IsComplex>
struct matrix_power_unwinder
{
static inline Scalar run(const Scalar& eival, const Scalar& eival0, int unwindingNumber)
@ -274,11 +272,8 @@ inline int matrix_power_get_pade_degree(long double normIminusT)
} // namespace internal
template<typename MatrixType, bool IsComplex=NumTraits<typename MatrixType::RealScalar>::IsComplex>
class MatrixPowerTriangularAtomic;
template<typename MatrixType>
class MatrixPowerTriangularAtomic<MatrixType,true>
class MatrixPowerTriangularAtomic
{
private:
enum {
@ -289,7 +284,7 @@ class MatrixPowerTriangularAtomic<MatrixType,true>
typedef typename MatrixType::RealScalar RealScalar;
typedef Array<Scalar,RowsAtCompileTime,1,ColMajor,MaxRowsAtCompileTime> ArrayType;
const MatrixType& m_T;
const MatrixType& m_A;
const MatrixType m_Id;
void computePade(int degree, const MatrixType& IminusT, MatrixType& res, RealScalar p) const;
@ -302,19 +297,19 @@ class MatrixPowerTriangularAtomic<MatrixType,true>
};
template<typename MatrixType>
MatrixPowerTriangularAtomic<MatrixType,true>::MatrixPowerTriangularAtomic(const MatrixType& T) :
m_T(T),
MatrixPowerTriangularAtomic<MatrixType>::MatrixPowerTriangularAtomic(const MatrixType& T) :
m_A(T),
m_Id(MatrixType::Identity(T.rows(), T.cols()))
{ eigen_assert(T.rows() == T.cols()); }
template<typename MatrixType>
void MatrixPowerTriangularAtomic<MatrixType,true>::compute(MatrixType& res, RealScalar p) const
void MatrixPowerTriangularAtomic<MatrixType>::compute(MatrixType& res, RealScalar p) const
{
switch (m_T.rows()) {
switch (m_A.rows()) {
case 0:
break;
case 1:
res(0,0) = std::pow(m_T(0,0), p);
res(0,0) = std::pow(m_A(0,0), p);
break;
case 2:
compute2x2(res, p);
@ -325,7 +320,7 @@ void MatrixPowerTriangularAtomic<MatrixType,true>::compute(MatrixType& res, Real
}
template<typename MatrixType>
void MatrixPowerTriangularAtomic<MatrixType,true>::computePade(int degree, const MatrixType& IminusT, MatrixType& res,
void MatrixPowerTriangularAtomic<MatrixType>::computePade(int degree, const MatrixType& IminusT, MatrixType& res,
RealScalar p) const
{
int i = degree<<1;
@ -338,33 +333,33 @@ void MatrixPowerTriangularAtomic<MatrixType,true>::computePade(int degree, const
}
template<typename MatrixType>
void MatrixPowerTriangularAtomic<MatrixType,true>::compute2x2(MatrixType& res, RealScalar p) const
void MatrixPowerTriangularAtomic<MatrixType>::compute2x2(MatrixType& res, RealScalar p) const
{
using std::abs;
using std::pow;
ArrayType logTdiag = m_T.diagonal().array().log();
res.coeffRef(0,0) = pow(m_T.coeff(0,0), p);
ArrayType logTdiag = m_A.diagonal().array().log();
res.coeffRef(0,0) = pow(m_A.coeff(0,0), p);
for (int i=1; i < m_T.cols(); ++i) {
res.coeffRef(i,i) = pow(m_T.coeff(i,i), p);
if (m_T.coeff(i-1,i-1) == m_T.coeff(i,i)) {
res.coeffRef(i-1,i) = p * pow(m_T.coeff(i-1,i), p-1);
for (int i=1; i < m_A.cols(); ++i) {
res.coeffRef(i,i) = pow(m_A.coeff(i,i), p);
if (m_A.coeff(i-1,i-1) == m_A.coeff(i,i)) {
res.coeffRef(i-1,i) = p * pow(m_A.coeff(i-1,i), p-1);
}
else if (2*abs(m_T.coeff(i-1,i-1)) < abs(m_T.coeff(i,i)) || 2*abs(m_T.coeff(i,i)) < abs(m_T.coeff(i-1,i-1))) {
res.coeffRef(i-1,i) = m_T.coeff(i-1,i) * (res.coeff(i,i)-res.coeff(i-1,i-1)) / (m_T.coeff(i,i)-m_T.coeff(i-1,i-1));
else if (2*abs(m_A.coeff(i-1,i-1)) < abs(m_A.coeff(i,i)) || 2*abs(m_A.coeff(i,i)) < abs(m_A.coeff(i-1,i-1))) {
res.coeffRef(i-1,i) = m_A.coeff(i-1,i) * (res.coeff(i,i)-res.coeff(i-1,i-1)) / (m_A.coeff(i,i)-m_A.coeff(i-1,i-1));
}
else {
int unwindingNumber = std::ceil((internal::imag(logTdiag[i]-logTdiag[i-1]) - M_PI) / (2*M_PI));
Scalar w = internal::matrix_power_unwinder<Scalar>::run(m_T.coeff(i,i), m_T.coeff(i-1,i-1), unwindingNumber);
res.coeffRef(i-1,i) = m_T.coeff(i-1,i) * RealScalar(2) * std::exp(RealScalar(0.5)*p*(logTdiag[i]+logTdiag[i-1])) *
std::sinh(p * w) / (m_T.coeff(i,i) - m_T.coeff(i-1,i-1));
Scalar w = internal::matrix_power_unwinder<Scalar>::run(m_A.coeff(i,i), m_A.coeff(i-1,i-1), unwindingNumber);
res.coeffRef(i-1,i) = m_A.coeff(i-1,i) * RealScalar(2) * std::exp(RealScalar(0.5)*p*(logTdiag[i]+logTdiag[i-1])) *
std::sinh(p * w) / (m_A.coeff(i,i) - m_A.coeff(i-1,i-1));
}
}
}
template<typename MatrixType>
void MatrixPowerTriangularAtomic<MatrixType,true>::computeBig(MatrixType& res, RealScalar p) const
void MatrixPowerTriangularAtomic<MatrixType>::computeBig(MatrixType& res, RealScalar p) const
{
const int digits = std::numeric_limits<RealScalar>::digits;
const RealScalar maxNormForPade = digits <= 24? 4.3386528e-1f: // sigle precision
@ -372,10 +367,10 @@ void MatrixPowerTriangularAtomic<MatrixType,true>::computeBig(MatrixType& res, R
digits <= 64? 2.4471944416607995472e-1L: // extended precision
digits <= 106? 1.1016843812851143391275867258512e-1L: // double-double
9.134603732914548552537150753385375e-2L; // quadruple precision
MatrixType IminusT, sqrtT, T=m_T;
MatrixType IminusT, sqrtT, T = m_A.template triangularView<Upper>();
RealScalar normIminusT;
int degree, degree2, numberOfSquareRoots=0;
bool hasExtraSquareRoot=false;
int degree, degree2, numberOfSquareRoots = 0;
bool hasExtraSquareRoot = false;
while (true) {
IminusT = m_Id - T;
@ -388,7 +383,7 @@ void MatrixPowerTriangularAtomic<MatrixType,true>::computeBig(MatrixType& res, R
hasExtraSquareRoot = true;
}
MatrixSquareRootTriangular<MatrixType>(T).compute(sqrtT);
T = sqrtT;
T = sqrtT.template triangularView<Upper>();
++numberOfSquareRoots;
}
computePade(degree, IminusT, res, p);

56
unsupported/test/matrix_power.cpp
View File

@ -9,6 +9,33 @@
#include "matrix_functions.h"
template <typename MatrixType, int IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex>
struct generateTriangularMatrix;
// for real matrices, make sure none of the eigenvalues are negative
template <typename MatrixType>
struct generateTriangularMatrix<MatrixType,0>
{
static void run(MatrixType& result, typename MatrixType::Index size)
{
result.resize(size, size);
result.template triangularView<Upper>() = MatrixType::Random(size, size);
for (typename MatrixType::Index i = 0; i < size; ++i)
result.coeffRef(i,i) = std::abs(result.coeff(i,i));
}
};
// for complex matrices, any matrix is fine
template <typename MatrixType>
struct generateTriangularMatrix<MatrixType,1>
{
static void run(MatrixType& result, typename MatrixType::Index size)
{
result.resize(size, size);
result.template triangularView<Upper>() = MatrixType::Random(size, size);
}
};
template<typename T>
void test2dRotation(double tol)
{
@ -59,7 +86,7 @@ void testExponentLaws(const MatrixType& m, double tol)
MatrixType m1, m2, m3, m4, m5;
RealScalar x, y;
for (int i=0; i<g_repeat; ++i) {
for (int i=0; i < g_repeat; ++i) {
generateTestMatrix<MatrixType>::run(m1, m.rows());
MatrixPower<MatrixType> mpow(m1);
@ -90,7 +117,7 @@ void testProduct(const MatrixType& m, const VectorType& v, double tol)
VectorType v1, v2, v3;
RealScalar p;
for (int i=0; i<g_repeat; ++i) {
for (int i=0; i < g_repeat; ++i) {
generateTestMatrix<MatrixType>::run(m1, m.rows());
MatrixPower<MatrixType> mpow(m1);
@ -99,7 +126,29 @@ void testProduct(const MatrixType& m, const VectorType& v, double tol)
v2.noalias() = mpow(p) * v1;
v3.noalias() = mpow(p).eval() * v1;
std::cout << "testMatrixVectorProduct: error powerm = " << relerr(v2, v3) << '\n';
std::cout << "testProduct: error powerm = " << relerr(v2, v3) << '\n';
VERIFY(v2.isApprox(v3, static_cast<RealScalar>(tol)));
}
}
template<typename MatrixType, typename VectorType>
void testTriangularProduct(const MatrixType& m, const VectorType& v, double tol)
{
typedef typename MatrixType::RealScalar RealScalar;
MatrixType m1;
VectorType v1, v2, v3;
RealScalar p;
for (int i=0; i < g_repeat; ++i) {
generateTriangularMatrix<MatrixType>::run(m1, m.rows());
MatrixPowerTriangular<MatrixType> mpow(m1);
v1 = VectorType::Random(v.rows(), v.cols());
p = internal::random<RealScalar>();
v2.noalias() = mpow(p) * v1;
v3.noalias() = mpow(p).eval() * v1;
std::cout << "testTriangularProduct: error powerm = " << relerr(v2, v3) << '\n';
VERIFY(v2.isApprox(v3, static_cast<RealScalar>(tol)));
}
}
@ -109,6 +158,7 @@ void testMatrixVector(const MatrixType& m, const VectorType& v, double tol)
{
testExponentLaws(m,tol);
testProduct(m,v,tol);
testTriangularProduct(m,v,tol);
}
void test_matrix_power()