diff --git a/unsupported/test/matrix_exponential.cpp b/unsupported/test/matrix_exponential.cpp index 86e942edb..61f30334d 100644 --- a/unsupported/test/matrix_exponential.cpp +++ b/unsupported/test/matrix_exponential.cpp @@ -133,7 +133,7 @@ void randomTest(const MatrixType& m, double tol) m1 = MatrixType::Random(rows, cols); m2 = ei_matrix_function(m1, expfn) * ei_matrix_function(-m1, expfn); - std::cout << "randomTest: error funm = " << relerr(identity, m2 * m3); + std::cout << "randomTest: error funm = " << relerr(identity, m2); VERIFY(identity.isApprox(m2, static_cast(tol))); m2 = ei_matrix_exponential(m1) * ei_matrix_exponential(-m1); diff --git a/unsupported/test/matrix_function.cpp b/unsupported/test/matrix_function.cpp index 7a1501da2..e40af7b4e 100644 --- a/unsupported/test/matrix_function.cpp +++ b/unsupported/test/matrix_function.cpp @@ -25,6 +25,17 @@ #include "main.h" #include +// Variant of VERIFY_IS_APPROX which uses absolute error instead of +// relative error. +#define VERIFY_IS_APPROX_ABS(a, b) VERIFY(test_isApprox_abs(a, b)) + +template +inline bool test_isApprox_abs(const Type1& a, const Type2& b) +{ + return ((a-b).array().abs() < test_precision()).all(); +} + + // Returns a matrix with eigenvalues clustered around 0, 1 and 2. template MatrixType randomMatrixWithRealEivals(const int size) @@ -37,7 +48,8 @@ MatrixType randomMatrixWithRealEivals(const int size) + ei_random() * Scalar(RealScalar(0.01)); } MatrixType A = MatrixType::Random(size, size); - return A.inverse() * diag * A; + HouseholderQR QRofA(A); + return QRofA.householderQ().inverse() * diag * QRofA.householderQ(); } template ::Scalar>::IsComplex> @@ -69,7 +81,8 @@ struct randomMatrixWithImagEivals } } MatrixType A = MatrixType::Random(size, size); - return A.inverse() * diag * A; + HouseholderQR QRofA(A); + return QRofA.householderQ().inverse() * diag * QRofA.householderQ(); } }; @@ -88,10 +101,12 @@ struct randomMatrixWithImagEivals + ei_random() * Scalar(RealScalar(0.01)); } MatrixType A = MatrixType::Random(size, size); - return A.inverse() * diag * A; + HouseholderQR QRofA(A); + return QRofA.householderQ().inverse() * diag * QRofA.householderQ(); } }; + template void testMatrixExponential(const MatrixType& A) { @@ -99,50 +114,45 @@ void testMatrixExponential(const MatrixType& A) typedef typename NumTraits::Real RealScalar; typedef std::complex ComplexScalar; - for (int i = 0; i < g_repeat; i++) { - VERIFY_IS_APPROX(ei_matrix_exponential(A), - ei_matrix_function(A, StdStemFunctions::exp)); - } + VERIFY_IS_APPROX(ei_matrix_exponential(A), + ei_matrix_function(A, StdStemFunctions::exp)); } template void testHyperbolicFunctions(const MatrixType& A) { - for (int i = 0; i < g_repeat; i++) { - MatrixType expA = ei_matrix_exponential(A); - MatrixType expmA = ei_matrix_exponential(-A); - VERIFY_IS_APPROX(ei_matrix_sinh(A), (expA - expmA) / 2); - VERIFY_IS_APPROX(ei_matrix_cosh(A), (expA + expmA) / 2); - } + // Need to use absolute error because of possible cancellation when + // adding/subtracting expA and expmA. + MatrixType expA = ei_matrix_exponential(A); + MatrixType expmA = ei_matrix_exponential(-A); + VERIFY_IS_APPROX_ABS(ei_matrix_sinh(A), (expA - expmA) / 2); + VERIFY_IS_APPROX_ABS(ei_matrix_cosh(A), (expA + expmA) / 2); } template void testGonioFunctions(const MatrixType& A) { - typedef ei_traits Traits; - typedef typename Traits::Scalar Scalar; + typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; typedef std::complex ComplexScalar; - typedef Matrix ComplexMatrix; + typedef Matrix ComplexMatrix; ComplexScalar imagUnit(0,1); ComplexScalar two(2,0); - for (int i = 0; i < g_repeat; i++) { - ComplexMatrix Ac = A.template cast(); - - ComplexMatrix exp_iA = ei_matrix_exponential(imagUnit * Ac); - ComplexMatrix exp_miA = ei_matrix_exponential(-imagUnit * Ac); - - MatrixType sinA = ei_matrix_sin(A); - ComplexMatrix sinAc = sinA.template cast(); - VERIFY_IS_APPROX(sinAc, (exp_iA - exp_miA) / (two*imagUnit)); - - MatrixType cosA = ei_matrix_cos(A); - ComplexMatrix cosAc = cosA.template cast(); - VERIFY_IS_APPROX(cosAc, (exp_iA + exp_miA) / 2); - } + ComplexMatrix Ac = A.template cast(); + + ComplexMatrix exp_iA = ei_matrix_exponential(imagUnit * Ac); + ComplexMatrix exp_miA = ei_matrix_exponential(-imagUnit * Ac); + + MatrixType sinA = ei_matrix_sin(A); + ComplexMatrix sinAc = sinA.template cast(); + VERIFY_IS_APPROX_ABS(sinAc, (exp_iA - exp_miA) / (two*imagUnit)); + + MatrixType cosA = ei_matrix_cos(A); + ComplexMatrix cosAc = cosA.template cast(); + VERIFY_IS_APPROX_ABS(cosAc, (exp_iA + exp_miA) / 2); } template