summaryrefslogtreecommitdiffhomepage
path: root/eigen/unsupported/test/matrix_power.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'eigen/unsupported/test/matrix_power.cpp')
-rw-r--r--eigen/unsupported/test/matrix_power.cpp171
1 files changed, 121 insertions, 50 deletions
diff --git a/eigen/unsupported/test/matrix_power.cpp b/eigen/unsupported/test/matrix_power.cpp
index b9d513b..7ccfacf 100644
--- a/eigen/unsupported/test/matrix_power.cpp
+++ b/eigen/unsupported/test/matrix_power.cpp
@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
-// Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net>
+// Copyright (C) 2012, 2013 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
@@ -9,35 +9,8 @@
#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)
+void test2dRotation(const T& tol)
{
Matrix<T,2,2> A, B, C;
T angle, c, s;
@@ -46,19 +19,19 @@ void test2dRotation(double tol)
MatrixPower<Matrix<T,2,2> > Apow(A);
for (int i=0; i<=20; ++i) {
- angle = pow(10, (i-10) / 5.);
+ angle = std::pow(T(10), (i-10) / T(5.));
c = std::cos(angle);
s = std::sin(angle);
B << c, s, -s, c;
- C = Apow(std::ldexp(angle,1) / M_PI);
+ C = Apow(std::ldexp(angle,1) / T(EIGEN_PI));
std::cout << "test2dRotation: i = " << i << " error powerm = " << relerr(C,B) << '\n';
- VERIFY(C.isApprox(B, static_cast<T>(tol)));
+ VERIFY(C.isApprox(B, tol));
}
}
template<typename T>
-void test2dHyperbolicRotation(double tol)
+void test2dHyperbolicRotation(const T& tol)
{
Matrix<std::complex<T>,2,2> A, B, C;
T angle, ch = std::cosh((T)1);
@@ -75,12 +48,26 @@ void test2dHyperbolicRotation(double tol)
C = Apow(angle);
std::cout << "test2dHyperbolicRotation: i = " << i << " error powerm = " << relerr(C,B) << '\n';
- VERIFY(C.isApprox(B, static_cast<T>(tol)));
+ VERIFY(C.isApprox(B, tol));
+ }
+}
+
+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), (i-10) / T(5.));
+ VERIFY(AngleAxis<T>(angle, v).matrix().isApprox(AngleAxis<T>(1,v).matrix().pow(angle), tol));
}
}
template<typename MatrixType>
-void testExponentLaws(const MatrixType& m, double tol)
+void testGeneral(const MatrixType& m, const typename MatrixType::RealScalar& tol)
{
typedef typename MatrixType::RealScalar RealScalar;
MatrixType m1, m2, m3, m4, m5;
@@ -97,37 +84,121 @@ void testExponentLaws(const MatrixType& m, double tol)
m4 = mpow(x+y);
m5.noalias() = m2 * m3;
- VERIFY(m4.isApprox(m5, static_cast<RealScalar>(tol)));
+ VERIFY(m4.isApprox(m5, tol));
m4 = mpow(x*y);
m5 = m2.pow(y);
- VERIFY(m4.isApprox(m5, static_cast<RealScalar>(tol)));
+ VERIFY(m4.isApprox(m5, tol));
m4 = (std::abs(x) * m1).pow(y);
m5 = std::pow(std::abs(x), y) * m3;
- VERIFY(m4.isApprox(m5, static_cast<RealScalar>(tol)));
+ VERIFY(m4.isApprox(m5, tol));
+ }
+}
+
+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 typename internal::conditional<IsComplex, TriangularView<MatrixType,Upper>, const MatrixType&>::type TriangularType;
+ typename internal::conditional< IsComplex, ComplexSchur<MatrixType>, RealSchur<MatrixType> >::type 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));
+ }
+}
+
+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));
}
}
typedef Matrix<double,3,3,RowMajor> Matrix3dRowMajor;
+typedef Matrix<long double,3,3> Matrix3e;
typedef Matrix<long double,Dynamic,Dynamic> MatrixXe;
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_9(test2dRotation<long double>(1e-13));
+ CALL_SUBTEST_9(test2dRotation<long double>(1e-13L));
CALL_SUBTEST_2(test2dHyperbolicRotation<double>(1e-14));
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(Matrix3dRowMajor(), 1e-13));
- CALL_SUBTEST_3(testExponentLaws(Matrix4cd(), 1e-13));
- CALL_SUBTEST_4(testExponentLaws(MatrixXd(8,8), 2e-12));
- 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(2,2), 1e-3)); // see bug 614
- CALL_SUBTEST_9(testExponentLaws(MatrixXe(7,7), 1e-13));
+ CALL_SUBTEST_9(test2dHyperbolicRotation<long double>(1e-14L));
+
+ CALL_SUBTEST_10(test3dRotation<double>(1e-13));
+ CALL_SUBTEST_11(test3dRotation<float>(1e-5));
+ 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-4));
+ CALL_SUBTEST_5(testGeneral(Matrix3cf(), 1e-4));
+ CALL_SUBTEST_8(testGeneral(Matrix4f(), 1e-4));
+ CALL_SUBTEST_6(testGeneral(MatrixXf(2,2), 1e-3)); // see bug 614
+ CALL_SUBTEST_9(testGeneral(MatrixXe(7,7), 1e-13L));
+ CALL_SUBTEST_10(testGeneral(Matrix3d(), 1e-13));
+ CALL_SUBTEST_11(testGeneral(Matrix3f(), 1e-4));
+ 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-4));
+ CALL_SUBTEST_5(testSingular(Matrix3cf(), 1e-4));
+ CALL_SUBTEST_8(testSingular(Matrix4f(), 1e-4));
+ CALL_SUBTEST_6(testSingular(MatrixXf(2,2), 1e-3));
+ CALL_SUBTEST_9(testSingular(MatrixXe(7,7), 1e-13L));
+ CALL_SUBTEST_10(testSingular(Matrix3d(), 1e-13));
+ CALL_SUBTEST_11(testSingular(Matrix3f(), 1e-4));
+ 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-4));
+ CALL_SUBTEST_5(testLogThenExp(Matrix3cf(), 1e-4));
+ CALL_SUBTEST_8(testLogThenExp(Matrix4f(), 1e-4));
+ CALL_SUBTEST_6(testLogThenExp(MatrixXf(2,2), 1e-3));
+ CALL_SUBTEST_9(testLogThenExp(MatrixXe(7,7), 1e-13L));
+ CALL_SUBTEST_10(testLogThenExp(Matrix3d(), 1e-13));
+ CALL_SUBTEST_11(testLogThenExp(Matrix3f(), 1e-4));
+ CALL_SUBTEST_12(testLogThenExp(Matrix3e(), 1e-13L));
}
generated by cgit v1.2.3 (git 2.39.1) at 2025-04-17 05:41:05 +0000