mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-18 14:34:17 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix doc and tidy up

Browse Source
This commit is contained in:
Chen-Pang He 2012-09-28 02:08:14 +08:00
parent 3b88216d42
commit ed18d6f2ad
3 changed files with 80 additions and 72 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -67,7 +67,6 @@ class MatrixPower : public MatrixPowerBase<MatrixPower<MatrixType>,MatrixType>
*
* \param[in] b a matrix with the same rows as A.
* \param[in] p exponent, a real scalar.
* \param[in] noalias
* \param[out] res \f$ A^p b \f$, where A is specified in the
* constructor.
*/
@ -75,16 +74,8 @@ class MatrixPower : public MatrixPowerBase<MatrixPower<MatrixType>,MatrixType>
void compute(const Derived& b, ResultType& res, RealScalar p);
private:
using Base::m_A;
static const int Rows = MatrixType::RowsAtCompileTime;
static const int Cols = MatrixType::ColsAtCompileTime;
static const int Options = MatrixType::Options;
static const int MaxRows = MatrixType::MaxRowsAtCompileTime;
static const int MaxCols = MatrixType::MaxColsAtCompileTime;
typedef Matrix<std::complex<RealScalar>,Rows,Cols,Options,MaxRows,MaxCols> ComplexMatrix;
typedef Matrix<std::complex<RealScalar>, RowsAtCompileTime, ColsAtCompileTime,
Options,MaxRowsAtCompileTime,MaxColsAtCompileTime> ComplexMatrix;
MatrixType m_tmp1, m_tmp2;
ComplexMatrix m_T, m_U, m_fT;
bool m_init;
@ -158,7 +149,7 @@ typename MatrixPower<MatrixType>::Base::RealScalar MatrixPower<MatrixType>::modf
m_U = schurOfA.matrixU();
m_init = true;
const Array<RealScalar,Rows,1,ColMajor,MaxRows> absTdiag = m_T.diagonal().array().abs();
const RealArray absTdiag = m_T.diagonal().array().abs();
maxAbsEival = absTdiag.maxCoeff();
minAbsEival = absTdiag.minCoeff();
}
@ -251,8 +242,8 @@ void MatrixPower<MatrixType>::computeFracPower(ResultType& res, RealScalar p)
{
if (p) {
MatrixPowerTriangularAtomic<ComplexMatrix>(m_T).compute(m_fT, p);
internal::recompose_complex_schur<NumTraits<Scalar>::IsComplex>::run(m_tmp1, m_fT, m_U);
res = m_tmp1 * res;
internal::recompose_complex_schur<NumTraits<Scalar>::IsComplex>::run(m_tmp2, m_fT, m_U);
res = m_tmp2 * res;
}
}
@ -347,9 +338,9 @@ class MatrixPowerReturnValue : public ReturnByValue<MatrixPowerReturnValue<Deriv
};
namespace internal {
template<typename MatrixType, typename Derived>
struct nested<MatrixPowerMatrixProduct<MatrixType,Derived> >
{ typedef typename MatrixPowerMatrixProduct<MatrixType,Derived>::PlainObject const& type; };
template<typename Lhs, typename Rhs>
struct nested<MatrixPowerMatrixProduct<Lhs,Rhs> >
{ typedef typename MatrixPowerMatrixProduct<Lhs,Rhs>::PlainObject const& type; };
template<typename Derived>
struct traits<MatrixPowerReturnValue<Derived> >

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

@ -97,17 +97,20 @@ inline int matrix_power_get_pade_degree(long double normIminusT)
}
} // namespace internal
template<typename MatrixType, int UpLo=Upper>
template<typename MatrixType, unsigned int Mode=Upper>
class MatrixPowerTriangularAtomic
{
private:
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
};
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef Array<Scalar,
EIGEN_SIZE_MIN_PREFER_FIXED(MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime),
1,ColMajor,
EIGEN_SIZE_MIN_PREFER_FIXED(MatrixType::MaxRowsAtCompileTime,MatrixType::MaxColsAtCompileTime)> ArrayType;
typedef Array<Scalar,RowsAtCompileTime,1,ColMajor,MaxRowsAtCompileTime> ArrayType;
const MatrixType& m_T;
const MatrixType m_Id;
void computePade(int degree, const MatrixType& IminusT, MatrixType& res, RealScalar p) const;
void compute2x2(MatrixType& res, RealScalar p) const;
@ -118,13 +121,14 @@ class MatrixPowerTriangularAtomic
void compute(MatrixType& res, RealScalar p) const;
};
template<typename MatrixType, int UpLo>
MatrixPowerTriangularAtomic<MatrixType,UpLo>::MatrixPowerTriangularAtomic(const MatrixType& T) :
m_T(T)
{ eigen_assert(T.rows() == T.cols()); }
template<typename MatrixType, unsigned int Mode>
MatrixPowerTriangularAtomic<MatrixType,Mode>::MatrixPowerTriangularAtomic(const MatrixType& T) :
m_T(T),
m_Id(MatrixType::Identity(T.rows(), T.cols()))
{ /* empty body */ }
template<typename MatrixType, int UpLo>
void MatrixPowerTriangularAtomic<MatrixType,UpLo>::compute(MatrixType& res, RealScalar p) const
template<typename MatrixType, unsigned int Mode>
void MatrixPowerTriangularAtomic<MatrixType,Mode>::compute(MatrixType& res, RealScalar p) const
{
switch (m_T.rows()) {
case 0:
@ -140,21 +144,21 @@ void MatrixPowerTriangularAtomic<MatrixType,UpLo>::compute(MatrixType& res, Real
}
}
template<typename MatrixType, int UpLo>
void MatrixPowerTriangularAtomic<MatrixType,UpLo>::computePade(int degree, const MatrixType& IminusT, MatrixType& res,
template<typename MatrixType, unsigned int Mode>
void MatrixPowerTriangularAtomic<MatrixType,Mode>::computePade(int degree, const MatrixType& IminusT, MatrixType& res,
RealScalar p) const
{
int i = degree<<1;
res = (p-degree) / ((i-1)<<1) * IminusT;
for (--i; i; --i) {
res = (MatrixType::Identity(m_T.rows(), m_T.cols()) + res).template triangularView<UpLo>()
.solve((i==1 ? -p : i&1 ? (-p-(i>>1))/(i<<1) : (p-(i>>1))/((i-1)<<1)) * IminusT).eval();
res = (m_Id + res).template triangularView<Mode>().solve((i==1 ? -p : i&1 ? (-p-(i>>1))/(i<<1) :
(p-(i>>1))/((i-1)<<1)) * IminusT).eval();
}
res += MatrixType::Identity(m_T.rows(), m_T.cols());
res += m_Id;
}
template<typename MatrixType, int UpLo>
void MatrixPowerTriangularAtomic<MatrixType,UpLo>::compute2x2(MatrixType& res, RealScalar p) const
template<typename MatrixType, unsigned int Mode>
void MatrixPowerTriangularAtomic<MatrixType,Mode>::compute2x2(MatrixType& res, RealScalar p) const
{
using std::abs;
using std::pow;
@ -180,8 +184,8 @@ void MatrixPowerTriangularAtomic<MatrixType,UpLo>::compute2x2(MatrixType& res, R
}
}
template<typename MatrixType, int UpLo>
void MatrixPowerTriangularAtomic<MatrixType,UpLo>::computeBig(MatrixType& res, RealScalar p) const
template<typename MatrixType, unsigned int Mode>
void MatrixPowerTriangularAtomic<MatrixType,Mode>::computeBig(MatrixType& res, RealScalar p) const
{
const int digits = std::numeric_limits<RealScalar>::digits;
const RealScalar maxNormForPade = digits <= 24? 4.3386528e-1f: // sigle precision
@ -217,9 +221,16 @@ void MatrixPowerTriangularAtomic<MatrixType,UpLo>::computeBig(MatrixType& res, R
}
#define EIGEN_MATRIX_POWER_PUBLIC_INTERFACE(Derived) \
typedef MatrixPowerBase<Derived<MatrixType>,MatrixType> Base; \
using typename Base::Scalar; \
using typename Base::RealScalar;
typedef MatrixPowerBase<Derived<MatrixType>, MatrixType> Base; \
using Base::RowsAtCompileTime; \
using Base::ColsAtCompileTime; \
using Base::Options; \
using Base::MaxRowsAtCompileTime; \
using Base::MaxColsAtCompileTime; \
typedef typename Base::Scalar Scalar; \
typedef typename Base::RealScalar RealScalar; \
typedef typename Base::RealArray RealArray; \
using Base::m_A;
#define EIGEN_MATRIX_POWER_PRODUCT_PUBLIC_INTERFACE(Derived) \
typedef MatrixPowerProductBase<Derived, Lhs, Rhs> Base; \
@ -254,25 +265,24 @@ struct traits<MatrixPowerProductBase<Derived,_Lhs,_Rhs> >
template<typename Derived, typename MatrixType>
class MatrixPowerBase
{
protected:
public:
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
Options = MatrixType::Options,
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef typename MatrixType::Index Index;
const MatrixType& m_A;
const bool m_del; // whether to delete the pointer at destruction
public:
explicit MatrixPowerBase(const MatrixType& A) :
m_A(A),
m_del(false)
{ /* empty body */ }
explicit MatrixPowerBase(const MatrixType& A)
: m_A(A), m_del(false) { }
template<typename OtherDerived>
explicit MatrixPowerBase(const MatrixBase<OtherDerived>& A) :
m_A(*new MatrixType(A)),
m_del(true)
{ /* empty body */ }
explicit MatrixPowerBase(const MatrixBase<OtherDerived>& A)
: m_A(*new MatrixType(A)), m_del(true) { }
~MatrixPowerBase()
{ if (m_del) delete &m_A; }
@ -286,6 +296,13 @@ class MatrixPowerBase
Index rows() const { return m_A.rows(); }
Index cols() const { return m_A.cols(); }
protected:
typedef Array<RealScalar,RowsAtCompileTime,1,ColMajor,MaxRowsAtCompileTime> RealArray;
const MatrixType& m_A;
private:
const bool m_del; // whether to delete the pointer at destruction
};
template<typename Derived, typename Lhs, typename Rhs>

38
unsupported/test/matrix_power.cpp
View File

@ -87,7 +87,7 @@ void testExponentLaws(const MatrixType& m, double tol)
}
template<typename MatrixType, typename VectorType>
void testMatrixVectorProduct(const MatrixType& m, const VectorType& v, double tol)
void testProduct(const MatrixType& m, const VectorType& v, double tol)
{
typedef typename MatrixType::RealScalar RealScalar;
MatrixType m1;
@ -108,9 +108,18 @@ void testMatrixVectorProduct(const MatrixType& m, const VectorType& v, double to
}
}
template<typename MatrixType, typename VectorType>
void testMatrixVector(const MatrixType& m, const VectorType& v, double tol)
{
testExponentLaws(m,tol);
testProduct(m,v,tol);
}
void test_matrix_power()
{
typedef Matrix<double,3,3,RowMajor> Matrix3dRowMajor;
typedef Matrix<long double,Dynamic,Dynamic> MatrixXe;
typedef Matrix<long double,Dynamic,1> VectorXe;
CALL_SUBTEST_2(test2dRotation<double>(1e-13));
CALL_SUBTEST_1(test2dRotation<float>(2e-5)); // was 1e-5, relaxed for clang 2.8 / linux / x86-64
@ -119,22 +128,13 @@ void test_matrix_power()
CALL_SUBTEST_1(test2dHyperbolicRotation<float>(1e-5));
CALL_SUBTEST_9(test2dHyperbolicRotation<long double>(1e-14));
CALL_SUBTEST_2(testExponentLaws(Matrix2d(), 1e-13));
CALL_SUBTEST_7(testExponentLaws(Matrix<double,3,3,RowMajor>(), 1e-13));
CALL_SUBTEST_3(testExponentLaws(Matrix4cd(), 1e-13));
CALL_SUBTEST_4(testExponentLaws(MatrixXd(8,8), 1e-13));
CALL_SUBTEST_1(testExponentLaws(Matrix2f(), 1e-4));
CALL_SUBTEST_5(testExponentLaws(Matrix3cf(), 1e-4));
CALL_SUBTEST_8(testExponentLaws(Matrix4f(), 1e-4));
CALL_SUBTEST_6(testExponentLaws(MatrixXf(8,8), 1e-4));
CALL_SUBTEST_2(testMatrixVectorProduct(Matrix2d(), Vector2d(), 1e-13));
CALL_SUBTEST_7(testMatrixVectorProduct(Matrix<double,3,3,RowMajor>(), Vector3d(), 1e-13));
CALL_SUBTEST_3(testMatrixVectorProduct(Matrix4cd(), Vector4cd(), 1e-13));
CALL_SUBTEST_4(testMatrixVectorProduct(MatrixXd(8,8), MatrixXd(8,2), 1e-13));
CALL_SUBTEST_1(testMatrixVectorProduct(Matrix2f(), Vector2f(), 1e-4));
CALL_SUBTEST_5(testMatrixVectorProduct(Matrix3cf(), Vector3cf(), 1e-4));
CALL_SUBTEST_8(testMatrixVectorProduct(Matrix4f(), Vector4f(), 1e-4));
CALL_SUBTEST_6(testMatrixVectorProduct(MatrixXf(8,8), VectorXf(8), 1e-4));
CALL_SUBTEST_9(testMatrixVectorProduct(MatrixXe(7,7), MatrixXe(7,9), 1e-13));
CALL_SUBTEST_2(testMatrixVector(Matrix2d(), Vector2d(), 1e-13));
CALL_SUBTEST_7(testMatrixVector(Matrix3dRowMajor(), MatrixXd(3,5), 1e-13));
CALL_SUBTEST_3(testMatrixVector(Matrix4cd(), Vector4cd(), 1e-13));
CALL_SUBTEST_4(testMatrixVector(MatrixXd(8,8), VectorXd(8), 1e-13));
CALL_SUBTEST_1(testMatrixVector(Matrix2f(), Vector2f(), 1e-4));
CALL_SUBTEST_5(testMatrixVector(Matrix3cf(), Vector3cf(), 1e-4));
CALL_SUBTEST_8(testMatrixVector(Matrix4f(), Vector4f(), 1e-4));
CALL_SUBTEST_6(testMatrixVector(MatrixXf(8,8), VectorXf(8), 1e-4));
CALL_SUBTEST_9(testMatrixVector(MatrixXe(7,7), VectorXe(7), 1e-13));
}