matrix_power.cpp File Reference

#include "matrix_functions.h"

Typedefs
typedef Matrix< double, 3, 3, RowMajor >	Matrix3dRowMajor
typedef Matrix< long double, 3, 3 >	Matrix3e
typedef Matrix< long double, Dynamic, Dynamic >	MatrixXe

Functions
template<typename T >
void	test2dRotation (const T &tol)
template<typename T >
void	test2dHyperbolicRotation (const T &tol)
template<typename T >
void	test3dRotation (const T &tol)
template<typename MatrixType >
void	testGeneral (const MatrixType &m, const typename MatrixType::RealScalar &tol)
template<typename MatrixType >
void	testSingular (const MatrixType &m_const, const typename MatrixType::RealScalar &tol)
template<typename MatrixType >
void	testLogThenExp (const MatrixType &m_const, const typename MatrixType::RealScalar &tol)
	EIGEN_DECLARE_TEST (matrix_power)

Typedef Documentation

Matrix3dRowMajor

typedef Matrix<double, 3, 3, RowMajor> Matrix3dRowMajor

Matrix3e

typedef Matrix<long double, 3, 3> Matrix3e

MatrixXe

typedef Matrix<long double, Dynamic, Dynamic> MatrixXe

Function Documentation

EIGEN_DECLARE_TEST()

EIGEN_DECLARE_TEST

(

matrix_power

)

                                 {
  CALL_SUBTEST_2(test2dRotation<double>(1e-13));
  CALL_SUBTEST_1(test2dRotation<float>(2e-5f));  // was 1e-5, relaxed for clang 2.8 / linux / x86-64
  CALL_SUBTEST_9(test2dRotation<long double>(1e-13L));
  CALL_SUBTEST_2(test2dHyperbolicRotation<double>(1e-14));
  CALL_SUBTEST_1(test2dHyperbolicRotation<float>(1e-5f));
  CALL_SUBTEST_9(test2dHyperbolicRotation<long double>(1e-14L));
 
  CALL_SUBTEST_10(test3dRotation<double>(1e-13));
  CALL_SUBTEST_11(test3dRotation<float>(1e-5f));
  CALL_SUBTEST_12(test3dRotation<long double>(1e-13L));
 
  CALL_SUBTEST_2(testGeneral(Matrix2d(), 1e-13));
  CALL_SUBTEST_7(testGeneral(Matrix3dRowMajor(), 1e-13));
  CALL_SUBTEST_3(testGeneral(Matrix4cd(), 1e-13));
  CALL_SUBTEST_4(testGeneral(MatrixXd(8, 8), 2e-12));
  CALL_SUBTEST_1(testGeneral(Matrix2f(), 1e-4f));
  CALL_SUBTEST_5(testGeneral(Matrix3cf(), 1e-4f));
  CALL_SUBTEST_8(testGeneral(Matrix4f(), 1e-4f));
  CALL_SUBTEST_6(testGeneral(MatrixXf(2, 2), 1e-3f));  // see bug 614
  CALL_SUBTEST_9(testGeneral(MatrixXe(7, 7), 1e-12L));
  CALL_SUBTEST_10(testGeneral(Matrix3d(), 1e-13));
  CALL_SUBTEST_11(testGeneral(Matrix3f(), 1e-4f));
  CALL_SUBTEST_12(testGeneral(Matrix3e(), 1e-13L));
 
  CALL_SUBTEST_2(testSingular(Matrix2d(), 1e-13));
  CALL_SUBTEST_7(testSingular(Matrix3dRowMajor(), 1e-13));
  CALL_SUBTEST_3(testSingular(Matrix4cd(), 1e-13));
  CALL_SUBTEST_4(testSingular(MatrixXd(8, 8), 2e-12));
  CALL_SUBTEST_1(testSingular(Matrix2f(), 1e-4f));
  CALL_SUBTEST_5(testSingular(Matrix3cf(), 1e-4f));
  CALL_SUBTEST_8(testSingular(Matrix4f(), 1e-4f));
  CALL_SUBTEST_6(testSingular(MatrixXf(2, 2), 1e-3f));
  CALL_SUBTEST_9(testSingular(MatrixXe(7, 7), 1e-12L));
  CALL_SUBTEST_10(testSingular(Matrix3d(), 1e-13));
  CALL_SUBTEST_11(testSingular(Matrix3f(), 1e-4f));
  CALL_SUBTEST_12(testSingular(Matrix3e(), 1e-13L));
 
  CALL_SUBTEST_2(testLogThenExp(Matrix2d(), 1e-13));
  CALL_SUBTEST_7(testLogThenExp(Matrix3dRowMajor(), 1e-13));
  CALL_SUBTEST_3(testLogThenExp(Matrix4cd(), 1e-13));
  CALL_SUBTEST_4(testLogThenExp(MatrixXd(8, 8), 2e-12));
  CALL_SUBTEST_1(testLogThenExp(Matrix2f(), 1e-4f));
  CALL_SUBTEST_5(testLogThenExp(Matrix3cf(), 1e-4f));
  CALL_SUBTEST_8(testLogThenExp(Matrix4f(), 1e-4f));
  CALL_SUBTEST_6(testLogThenExp(MatrixXf(2, 2), 1e-3f));
  CALL_SUBTEST_9(testLogThenExp(MatrixXe(7, 7), 1e-12L));
  CALL_SUBTEST_10(testLogThenExp(Matrix3d(), 1e-13));
  CALL_SUBTEST_11(testLogThenExp(Matrix3f(), 1e-4f));
  CALL_SUBTEST_12(testLogThenExp(Matrix3e(), 1e-13L));
}

References CALL_SUBTEST_1, CALL_SUBTEST_10, CALL_SUBTEST_11, CALL_SUBTEST_12, CALL_SUBTEST_2, CALL_SUBTEST_3, CALL_SUBTEST_4, CALL_SUBTEST_5, CALL_SUBTEST_6, CALL_SUBTEST_7, CALL_SUBTEST_8, CALL_SUBTEST_9, e(), testGeneral(), testLogThenExp(), and testSingular().

test2dHyperbolicRotation()

template<typename T >

void test2dHyperbolicRotation

(

const T &

tol

)

                                            {
  Matrix<std::complex<T>, 2, 2> A, B, C;
  T angle, ch = std::cosh((T)1);
  std::complex<T> ish(0, std::sinh((T)1));
 
  A << ch, ish, -ish, ch;
  MatrixPower<Matrix<std::complex<T>, 2, 2> > Apow(A);
 
  for (int i = 0; i <= 20; ++i) {
    angle = std::ldexp(static_cast<T>(i - 10), -1);
    ch = std::cosh(angle);
    ish = std::complex<T>(0, std::sinh(angle));
    B << ch, ish, -ish, ch;
 
    C = Apow(angle);
    std::cout << "test2dHyperbolicRotation: i = " << i << "   error powerm = " << relerr(C, B) << '\n';
    VERIFY(C.isApprox(B, tol));
  }
}

References Jeffery_Solution::angle(), Eigen::bfloat16_impl::cosh(), i, relerr(), Eigen::bfloat16_impl::sinh(), and VERIFY.

test2dRotation()

template<typename T >

void test2dRotation

(

const T &

tol

)

                                  {
  Matrix<T, 2, 2> A, B, C;
  T angle, c, s;
 
  A << 0, 1, -1, 0;
  MatrixPower<Matrix<T, 2, 2> > Apow(A);
 
  for (int i = 0; i <= 20; ++i) {
    angle = std::pow(T(10), T(i - 10) / T(5.));
    c = std::cos(angle);
    s = std::sin(angle);
    B << c, s, -s, c;
 
    C = Apow(std::ldexp(angle, 1) / T(EIGEN_PI));
    std::cout << "test2dRotation: i = " << i << "   error powerm = " << relerr(C, B) << '\n';
    VERIFY(C.isApprox(B, tol));
  }
}

References Jeffery_Solution::angle(), calibrate::c, cos(), EIGEN_PI, i, Eigen::bfloat16_impl::pow(), relerr(), s, sin(), and VERIFY.

test3dRotation()

template<typename T >

void test3dRotation

(

const T &

tol

)

                                  {
  Matrix<T, 3, 1> v;
  T angle;
 
  for (int i = 0; i <= 20; ++i) {
    v = Matrix<T, 3, 1>::Random();
    v.normalize();
    angle = std::pow(T(10), T(i - 10) / T(5.));
    VERIFY(AngleAxis<T>(angle, v).matrix().isApprox(AngleAxis<T>(1, v).matrix().pow(angle), tol));
  }
}

References Jeffery_Solution::angle(), i, Eigen::internal::isApprox(), matrix(), Eigen::bfloat16_impl::pow(), Eigen::ArrayBase< Derived >::pow(), v, and VERIFY.

testGeneral()

template<typename MatrixType >

void testGeneral	(	const MatrixType &	m,
		const typename MatrixType::RealScalar &	tol
	)

                                                                                {
  typedef typename MatrixType::RealScalar RealScalar;
  MatrixType m1, m2, m3, m4, m5;
  RealScalar x, y;
 
  for (int i = 0; i < g_repeat; ++i) {
    generateTestMatrix<MatrixType>::run(m1, m.rows());
    MatrixPower<MatrixType> mpow(m1);
 
    x = internal::random<RealScalar>();
    y = internal::random<RealScalar>();
    m2 = mpow(x);
    m3 = mpow(y);
 
    m4 = mpow(x + y);
    m5.noalias() = m2 * m3;
    VERIFY(m4.isApprox(m5, tol));
 
    m4 = mpow(x * y);
    m5 = m2.pow(y);
    VERIFY(m4.isApprox(m5, tol));
 
    m4 = (std::abs(x) * m1).pow(y);
    m5 = std::pow(std::abs(x), y) * m3;
    VERIFY(m4.isApprox(m5, tol));
  }
}

References abs(), Eigen::g_repeat, i, m, m1, m2(), Eigen::bfloat16_impl::pow(), Eigen::ArrayBase< Derived >::pow(), run(), VERIFY, plotDoE::x, and y.

Referenced by EIGEN_DECLARE_TEST().

testLogThenExp()

template<typename MatrixType >

void testLogThenExp	(	const MatrixType &	m_const,
		const typename MatrixType::RealScalar &	tol
	)

                                                                                         {
  // we need to pass by reference in order to prevent errors with
  // MSVC for aligned data types ...
  MatrixType& m = const_cast<MatrixType&>(m_const);
 
  typedef typename MatrixType::Scalar Scalar;
  Scalar x;
 
  for (int i = 0; i < g_repeat; ++i) {
    generateTestMatrix<MatrixType>::run(m, m.rows());
    x = internal::random<Scalar>();
    VERIFY(m.pow(x).isApprox((x * m.log()).exp(), tol));
  }
}

References Eigen::g_repeat, i, m, run(), VERIFY, and plotDoE::x.

Referenced by EIGEN_DECLARE_TEST().

testSingular()

template<typename MatrixType >

void testSingular	(	const MatrixType &	m_const,
		const typename MatrixType::RealScalar &	tol
	)

                                                                                       {
  // we need to pass by reference in order to prevent errors with
  // MSVC for aligned data types ...
  MatrixType& m = const_cast<MatrixType&>(m_const);
 
  const int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex;
  typedef std::conditional_t<IsComplex, TriangularView<MatrixType, Upper>, const MatrixType&> TriangularType;
  std::conditional_t<IsComplex, ComplexSchur<MatrixType>, RealSchur<MatrixType> > schur;
  MatrixType T;
 
  for (int i = 0; i < g_repeat; ++i) {
    m.setRandom();
    m.col(0).fill(0);
 
    schur.compute(m);
    T = schur.matrixT();
    const MatrixType& U = schur.matrixU();
    processTriangularMatrix<MatrixType>::run(m, T, U);
    MatrixPower<MatrixType> mpow(m);
 
    T = T.sqrt();
    VERIFY(mpow(0.5L).isApprox(U * (TriangularType(T) * U.adjoint()), tol));
 
    T = T.sqrt();
    VERIFY(mpow(0.25L).isApprox(U * (TriangularType(T) * U.adjoint()), tol));
 
    T = T.sqrt();
    VERIFY(mpow(0.125L).isApprox(U * (TriangularType(T) * U.adjoint()), tol));
  }
}

References Eigen::g_repeat, i, Eigen::internal::isApprox(), IsComplex, m, processTriangularMatrix< MatrixType, IsComplex >::run(), schur(), RachelsAdvectionDiffusion::U, and VERIFY.

Referenced by EIGEN_DECLARE_TEST().