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

merge

Browse Source
This commit is contained in:
Jitse Niesen 2012-08-27 22:57:39 +01:00
commit edc7a09ee7
8 changed files with 234 additions and 114 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
unsupported/Eigen/MatrixFunctions
View File

@ -228,15 +228,15 @@ const MatrixPowerReturnValue<Derived, ExponentType> MatrixBase<Derived>::pow(con
\endcode \endcode
\param[in] M base of the matrix power, should be a square matrix. \param[in] M base of the matrix power, should be a square matrix.
\param[in] p exponent of the matrix power, should be an integer or \param[in] p exponent of the matrix power, should be real.
the same type as the real scalar in \p M.
The matrix power \f$ M^p \f$ is defined as \f$ \exp(p \log(M)) \f$, The matrix power \f$ M^p \f$ is defined as \f$ \exp(p \log(M)) \f$,
where exp denotes the matrix exponential, and log denotes the matrix where exp denotes the matrix exponential, and log denotes the matrix
logarithm. logarithm.
The matrix \f$ M \f$ should meet the conditions to be an argument of The matrix \f$ M \f$ should meet the conditions to be an argument of
matrix logarithm. matrix logarithm. If \p p is neither an integer nor the real scalar
type of \p M, it is casted into the real scalar type of \p M.
This function computes the matrix logarithm using the This function computes the matrix logarithm using the
Schur-Pad&eacute; algorithm as implemented by MatrixBase::pow(). Schur-Pad&eacute; algorithm as implemented by MatrixBase::pow().

18
unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h
View File

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

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

@ -21,16 +21,31 @@ namespace Eigen {
* *
* \brief Class for computing matrix powers. * \brief Class for computing matrix powers.
* *
* \tparam MatrixType type of the base, expected to be an instantiation * \tparam MatrixType type of the base, expected to be an instantiation
* of the Matrix class template. * of the Matrix class template.
* \tparam RealScalar type of the exponent, a real scalar. * \tparam ExponentType type of the exponent, a real scalar.
* \tparam PlainObject type of the multiplier. * \tparam PlainObject type of the multiplier.
* \tparam IsInteger used internally to select correct specialization. * \tparam IsInteger used internally to select correct specialization.
*/ */
template <typename MatrixType, typename RealScalar, typename PlainObject = MatrixType, template <typename MatrixType, typename ExponentType, typename PlainObject = MatrixType,
int IsInteger = NumTraits<RealScalar>::IsInteger> int IsInteger = NumTraits<ExponentType>::IsInteger>
class MatrixPower class MatrixPower
{ {
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::RealScalar RealScalar;
typedef std::complex<RealScalar> ComplexScalar;
typedef typename MatrixType::Index Index;
typedef Matrix<ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols> ComplexMatrix;
typedef Array<ComplexScalar, Rows, 1, ColMajor, MaxRows> ComplexArray;
public: public:
/** /**
* \brief Constructor. * \brief Constructor.
@ -39,7 +54,7 @@ class MatrixPower
* \param[in] p the exponent of the matrix power. * \param[in] p the exponent of the matrix power.
* \param[in] b the multiplier. * \param[in] b the multiplier.
*/ */
MatrixPower(const MatrixType& A, const RealScalar& p, const PlainObject& b) : MatrixPower(const MatrixType& A, RealScalar p, const PlainObject& b) :
m_A(A), m_A(A),
m_p(p), m_p(p),
m_b(b), m_b(b),
@ -55,19 +70,6 @@ class MatrixPower
template <typename ResultType> void compute(ResultType& result); template <typename ResultType> void compute(ResultType& result);
private: 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 std::complex<RealScalar> ComplexScalar;
typedef typename MatrixType::Index Index;
typedef Matrix<ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols> ComplexMatrix;
typedef Array<ComplexScalar, Rows, 1, ColMajor, MaxRows> ComplexArray;
/** /**
* \brief Compute the matrix power. * \brief Compute the matrix power.
* *
@ -112,8 +114,8 @@ class MatrixPower
*/ */
void getFractionalExponent(); void getFractionalExponent();
/** \brief Compute atanh (inverse hyperbolic tangent). */ /** \brief Compute atanh (inverse hyperbolic tangent) for \f$ y / x \f$. */
ComplexScalar atanh(const ComplexScalar& x); ComplexScalar atanh2(const ComplexScalar& y, const ComplexScalar& x);
/** \brief Compute power of 2x2 triangular matrix. */ /** \brief Compute power of 2x2 triangular matrix. */
void compute2x2(const RealScalar& p); void compute2x2(const RealScalar& p);
@ -237,9 +239,9 @@ class MatrixPower<MatrixType, IntExponent, PlainObject, 1>
/******* Specialized for real exponents *******/ /******* Specialized for real exponents *******/
template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger> template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger>
template <typename ResultType> template <typename ResultType>
void MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::compute(ResultType& result) void MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::compute(ResultType& result)
{ {
using std::floor; using std::floor;
using std::pow; using std::pow;
@ -264,9 +266,9 @@ void MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::compute(ResultTyp
} }
} }
template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger> template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger>
template <typename ResultType> template <typename ResultType>
void MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::computeIntPower(ResultType& result) void MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::computeIntPower(ResultType& result)
{ {
if (m_dimb > m_dimA) { if (m_dimb > m_dimA) {
MatrixType tmp = MatrixType::Identity(m_A.rows(), m_A.cols()); MatrixType tmp = MatrixType::Identity(m_A.rows(), m_A.cols());
@ -278,9 +280,9 @@ void MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::computeIntPower(R
} }
} }
template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger> template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger>
template <typename ResultType> template <typename ResultType>
void MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::computeChainProduct(ResultType& result) void MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::computeChainProduct(ResultType& result)
{ {
using std::frexp; using std::frexp;
using std::ldexp; using std::ldexp;
@ -312,8 +314,8 @@ void MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::computeChainProdu
result = m_tmp * result; result = m_tmp * result;
} }
template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger> template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger>
int MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::computeCost(RealScalar p) int MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::computeCost(RealScalar p)
{ {
using std::frexp; using std::frexp;
using std::ldexp; using std::ldexp;
@ -326,25 +328,25 @@ int MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::computeCost(RealSc
return cost; return cost;
} }
template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger> template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger>
template <typename ResultType> template <typename ResultType>
void MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::partialPivLuSolve(ResultType& result, RealScalar p) void MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::partialPivLuSolve(ResultType& result, RealScalar p)
{ {
const PartialPivLU<MatrixType> Asolver(m_A); const PartialPivLU<MatrixType> Asolver(m_A);
for (; p >= RealScalar(1); p--) for (; p >= RealScalar(1); p--)
result = Asolver.solve(result); result = Asolver.solve(result);
} }
template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger> template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger>
void MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::computeSchurDecomposition() void MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::computeSchurDecomposition()
{ {
const ComplexSchur<MatrixType> schurOfA(m_A); const ComplexSchur<MatrixType> schurOfA(m_A);
m_T = schurOfA.matrixT(); m_T = schurOfA.matrixT();
m_U = schurOfA.matrixU(); m_U = schurOfA.matrixU();
} }
template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger> template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger>
void MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::getFractionalExponent() void MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::getFractionalExponent()
{ {
using std::pow; using std::pow;
@ -373,21 +375,24 @@ void MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::getFractionalExpo
} }
} }
template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger> template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger>
std::complex<RealScalar> MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::atanh(const ComplexScalar& x) std::complex<typename MatrixType::RealScalar>
MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::atanh2(const ComplexScalar& y, const ComplexScalar& x)
{ {
using std::abs; using std::abs;
using std::log; using std::log;
using std::sqrt; using std::sqrt;
if (abs(x) > sqrt(NumTraits<RealScalar>::epsilon())) const ComplexScalar z = y / x;
return RealScalar(0.5) * log((RealScalar(1) + x) / (RealScalar(1) - x));
if (abs(z) > sqrt(NumTraits<RealScalar>::epsilon()))
return RealScalar(0.5) * log((x + y) / (x - y));
else else
return x + x*x*x / RealScalar(3); return z + z*z*z / RealScalar(3);
} }
template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger> template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger>
void MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::compute2x2(const RealScalar& p) void MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::compute2x2(const RealScalar& p)
{ {
using std::abs; using std::abs;
using std::ceil; using std::ceil;
@ -414,15 +419,15 @@ void MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::compute2x2(const
else { else {
// computation in previous branch is inaccurate if abs(m_T(j,j)) \approx abs(m_T(i,i)) // 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))); 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))) + ComplexScalar(0, M_PI * unwindingNumber); w = 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])) * 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)); sinh(p * w) / (m_T(j,j) - m_T(i,i));
} }
} }
} }
template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger> template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger>
void MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::computeBig() void MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::computeBig()
{ {
using std::ldexp; using std::ldexp;
const int digits = std::numeric_limits<RealScalar>::digits; const int digits = std::numeric_limits<RealScalar>::digits;
@ -458,8 +463,8 @@ void MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::computeBig()
compute2x2(m_pfrac); compute2x2(m_pfrac);
} }
template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger> template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger>
inline int MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::getPadeDegree(float normIminusT) inline int MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::getPadeDegree(float normIminusT)
{ {
const float maxNormForPade[] = { 2.7996156e-1f /* degree = 3 */ , 4.3268868e-1f }; const float maxNormForPade[] = { 2.7996156e-1f /* degree = 3 */ , 4.3268868e-1f };
int degree = 3; int degree = 3;
@ -469,8 +474,8 @@ inline int MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::getPadeDegr
return degree; return degree;
} }
template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger> template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger>
inline int MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::getPadeDegree(double normIminusT) inline int MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::getPadeDegree(double normIminusT)
{ {
const double maxNormForPade[] = { 1.882832775783710e-2 /* degree = 3 */ , 6.036100693089536e-2, const double maxNormForPade[] = { 1.882832775783710e-2 /* degree = 3 */ , 6.036100693089536e-2,
1.239372725584857e-1, 1.998030690604104e-1, 2.787629930861592e-1 }; 1.239372725584857e-1, 1.998030690604104e-1, 2.787629930861592e-1 };
@ -481,8 +486,8 @@ inline int MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::getPadeDegr
return degree; return degree;
} }
template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger> template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger>
inline int MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::getPadeDegree(long double normIminusT) inline int MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::getPadeDegree(long double normIminusT)
{ {
#if LDBL_MANT_DIG == 53 #if LDBL_MANT_DIG == 53
const int maxPadeDegree = 7; const int maxPadeDegree = 7;
@ -514,8 +519,8 @@ inline int MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::getPadeDegr
break; break;
return degree; return degree;
} }
template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger> template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger>
void MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::computePade(const int& degree, const ComplexMatrix& IminusT) void MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::computePade(const int& degree, const ComplexMatrix& IminusT)
{ {
int i = degree << 1; int i = degree << 1;
m_fT = coeff(i) * IminusT; m_fT = coeff(i) * IminusT;
@ -526,8 +531,8 @@ void MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::computePade(const
m_fT += ComplexMatrix::Identity(m_A.rows(), m_A.cols()); m_fT += ComplexMatrix::Identity(m_A.rows(), m_A.cols());
} }
template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger> template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger>
inline RealScalar MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::coeff(const int& i) inline typename MatrixType::RealScalar MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::coeff(const int& i)
{ {
if (i == 1) if (i == 1)
return -m_pfrac; return -m_pfrac;
@ -537,13 +542,13 @@ inline RealScalar MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::coef
return (m_pfrac - RealScalar(i >> 1)) / RealScalar(i-1 << 1); return (m_pfrac - RealScalar(i >> 1)) / RealScalar(i-1 << 1);
} }
template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger> template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger>
void MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::computeTmp(RealScalar) void MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::computeTmp(RealScalar)
{ m_tmp = (m_U * m_fT * m_U.adjoint()).real(); } { m_tmp = (m_U * m_fT * m_U.adjoint()).real(); }
template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger> template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger>
void MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::computeTmp(ComplexScalar) void MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::computeTmp(ComplexScalar)
{ m_tmp = (m_U * m_fT * m_U.adjoint()).eval(); } { m_tmp = m_U * m_fT * m_U.adjoint(); }
/******* Specialized for integral exponents *******/ /******* Specialized for integral exponents *******/

1
unsupported/test/CMakeLists.txt
View File

@ -33,6 +33,7 @@ endif()
ei_add_test(matrix_exponential) ei_add_test(matrix_exponential)
ei_add_test(matrix_function) ei_add_test(matrix_function)
ei_add_test(matrix_power)
ei_add_test(matrix_square_root) ei_add_test(matrix_square_root)
ei_add_test(alignedvector3) ei_add_test(alignedvector3)
ei_add_test(FFT) ei_add_test(FFT)

12
unsupported/test/matrix_exponential.cpp
View File

@ -7,8 +7,7 @@
// Public License v. 2.0. If a copy of the MPL was not distributed // Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#include "main.h" #include "matrix_functions.h"
#include <unsupported/Eigen/MatrixFunctions>
double binom(int n, int k) double binom(int n, int k)
{ {
@ -18,12 +17,6 @@ double binom(int n, int k)
return res; return res;
} }
template <typename Derived, typename OtherDerived>
double relerr(const MatrixBase<Derived>& A, const MatrixBase<OtherDerived>& B)
{
return std::sqrt((A - B).cwiseAbs2().sum() / (std::min)(A.cwiseAbs2().sum(), B.cwiseAbs2().sum()));
}
template <typename T> template <typename T>
T expfn(T x, int) T expfn(T x, int)
{ {
@ -109,8 +102,7 @@ void randomTest(const MatrixType& m, double tol)
*/ */
typename MatrixType::Index rows = m.rows(); typename MatrixType::Index rows = m.rows();
typename MatrixType::Index cols = m.cols(); typename MatrixType::Index cols = m.cols();
MatrixType m1(rows, cols), m2(rows, cols), m3(rows, cols), MatrixType m1(rows, cols), m2(rows, cols), identity = MatrixType::Identity(rows, cols);
identity = MatrixType::Identity(rows, rows);
typedef typename NumTraits<typename internal::traits<MatrixType>::Scalar>::Real RealScalar; typedef typename NumTraits<typename internal::traits<MatrixType>::Scalar>::Real RealScalar;

47
unsupported/test/matrix_functions.h Normal file
View File

@ -0,0 +1,47 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009-2011 Jitse Niesen <jitse@maths.leeds.ac.uk>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#include "main.h"
#include <unsupported/Eigen/MatrixFunctions>
template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
struct generateTestMatrix;
// for real matrices, make sure none of the eigenvalues are negative
template <typename MatrixType>
struct generateTestMatrix<MatrixType,0>
{
static void run(MatrixType& result, typename MatrixType::Index size)
{
MatrixType mat = MatrixType::Random(size, size);
EigenSolver<MatrixType> es(mat);
typename EigenSolver<MatrixType>::EigenvalueType eivals = es.eigenvalues();
for (typename MatrixType::Index i = 0; i < size; ++i) {
if (eivals(i).imag() == 0 && eivals(i).real() < 0)
eivals(i) = -eivals(i);
}
result = (es.eigenvectors() * eivals.asDiagonal() * es.eigenvectors().inverse()).real();
}
};
// for complex matrices, any matrix is fine
template <typename MatrixType>
struct generateTestMatrix<MatrixType,1>
{
static void run(MatrixType& result, typename MatrixType::Index size)
{
result = MatrixType::Random(size, size);
}
};
template <typename Derived, typename OtherDerived>
double relerr(const MatrixBase<Derived>& A, const MatrixBase<OtherDerived>& B)
{
return std::sqrt((A - B).cwiseAbs2().sum() / (std::min)(A.cwiseAbs2().sum(), B.cwiseAbs2().sum()));
}

104
unsupported/test/matrix_power.cpp Normal file
View File

@ -0,0 +1,104 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#include "matrix_functions.h"
template <typename T>
void test2dRotation(double tol)
{
Matrix<T,2,2> A, B, C;
T angle, c, s;
A << 0, 1, -1, 0;
for (int i = 0; i <= 20; i++) {
angle = pow(10, (i-10) / 5.);
c = std::cos(angle);
s = std::sin(angle);
B << c, s, -s, c;
C = A.pow(std::ldexp(angle, 1) / M_PI);
std::cout << "test2dRotation: i = " << i << " error powerm = " << relerr(C, B) << "\n";
VERIFY(C.isApprox(B, T(tol)));
}
}
template <typename T>
void test2dHyperbolicRotation(double tol)
{
Matrix<std::complex<T>,2,2> A, B, C;
T angle, ch = std::cosh(1);
std::complex<T> ish(0, std::sinh(1));
A << ch, ish, -ish, ch;
for (int i = 0; i <= 20; i++) {
angle = std::ldexp(T(i-10), -1);
ch = std::cosh(angle);
ish = std::complex<T>(0, std::sinh(angle));
B << ch, ish, -ish, ch;
C = A.pow(angle);
std::cout << "test2dHyperbolicRotation: i = " << i << " error powerm = " << relerr(C, B) << "\n";
VERIFY(C.isApprox(B, T(tol)));
}
}
template <typename MatrixType>
void testExponentLaws(const MatrixType& m, double tol)
{
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
typename MatrixType::Index rows = m.rows();
typename MatrixType::Index cols = m.cols();
MatrixType m1, m1x, m1y, m2, m3;
RealScalar x = internal::random<RealScalar>(), y = internal::random<RealScalar>();
double err[3];
for(int i = 0; i < g_repeat; i++) {
generateTestMatrix<MatrixType>::run(m1, m.rows());
m1x = m1.pow(x);
m1y = m1.pow(y);
m2 = m1.pow(x + y);
m3 = m1x * m1y;
err[0] = relerr(m2, m3);
VERIFY(m2.isApprox(m3, static_cast<RealScalar>(tol)));
m2 = m1.pow(x * y);
m3 = m1x.pow(y);
err[1] = relerr(m2, m3);
VERIFY(m2.isApprox(m3, static_cast<RealScalar>(tol)));
m2 = (std::abs(x) * m1).pow(y);
m3 = std::pow(std::abs(x), y) * m1y;
err[2] = relerr(m2, m3);
VERIFY(m2.isApprox(m3, static_cast<RealScalar>(tol)));
std::cout << "testExponentLaws: error powerm = " << err[0] << " " << err[1] << " " << err[2] << "\n";
}
}
void test_matrix_power()
{
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
CALL_SUBTEST_8(test2dRotation<long double>(1e-13));
CALL_SUBTEST_2(test2dHyperbolicRotation<double>(1e-14));
CALL_SUBTEST_1(test2dHyperbolicRotation<float>(1e-5));
CALL_SUBTEST_8(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_1(testExponentLaws(Matrix4f(), 1e-4));
CALL_SUBTEST_6(testExponentLaws(MatrixXf(8,8), 1e-4));
CALL_SUBTEST_9(testExponentLaws(Matrix<long double,Dynamic,Dynamic>(7,7), 1e-13));
}

33
unsupported/test/matrix_square_root.cpp
View File

@ -7,38 +7,7 @@
// Public License v. 2.0. If a copy of the MPL was not distributed // Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#include "main.h" #include "matrix_functions.h"
#include <unsupported/Eigen/MatrixFunctions>
template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
struct generateTestMatrix;
// for real matrices, make sure none of the eigenvalues are negative
template <typename MatrixType>
struct generateTestMatrix<MatrixType,0>
{
static void run(MatrixType& result, typename MatrixType::Index size)
{
MatrixType mat = MatrixType::Random(size, size);
EigenSolver<MatrixType> es(mat);
typename EigenSolver<MatrixType>::EigenvalueType eivals = es.eigenvalues();
for (typename MatrixType::Index i = 0; i < size; ++i) {
if (eivals(i).imag() == 0 && eivals(i).real() < 0)
eivals(i) = -eivals(i);
}
result = (es.eigenvectors() * eivals.asDiagonal() * es.eigenvectors().inverse()).real();
}
};
// for complex matrices, any matrix is fine
template <typename MatrixType>
struct generateTestMatrix<MatrixType,1>
{
static void run(MatrixType& result, typename MatrixType::Index size)
{
result = MatrixType::Random(size, size);
}
};
template<typename MatrixType> template<typename MatrixType>
void testMatrixSqrt(const MatrixType& m) void testMatrixSqrt(const MatrixType& m)