mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-04-06 19:10:36 +08:00
Code Issues Packages Projects Releases Wiki Activity

Make MatrixFunctions tests more robust.

Browse Source 
* Use absolute error instead of relative error.
* Test on well-conditioned matrices.
* Do not repeat the same test g_repeat times (bug fix).
* Correct diagnostic output in matrix_exponential.cpp .
This commit is contained in:
Jitse Niesen 2010-03-01 12:05:57 +00:00
parent f1f3c30ddc
commit 2d7bd1ec91
2 changed files with 42 additions and 32 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
unsupported/test/matrix_exponential.cpp
View File

@ -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<RealScalar>(tol)));
m2 = ei_matrix_exponential(m1) * ei_matrix_exponential(-m1);

72
unsupported/test/matrix_function.cpp
View File

@ -25,6 +25,17 @@
#include "main.h"
#include <unsupported/Eigen/MatrixFunctions>
// 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<typename Type1, typename Type2>
inline bool test_isApprox_abs(const Type1& a, const Type2& b)
{
return ((a-b).array().abs() < test_precision<typename Type1::RealScalar>()).all();
}
// Returns a matrix with eigenvalues clustered around 0, 1 and 2.
template<typename MatrixType>
MatrixType randomMatrixWithRealEivals(const int size)
@ -37,7 +48,8 @@ MatrixType randomMatrixWithRealEivals(const int size)
+ ei_random<Scalar>() * Scalar(RealScalar(0.01));
}
MatrixType A = MatrixType::Random(size, size);
return A.inverse() * diag * A;
HouseholderQR<MatrixType> QRofA(A);
return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
}
template <typename MatrixType, int IsComplex = NumTraits<typename ei_traits<MatrixType>::Scalar>::IsComplex>
@ -69,7 +81,8 @@ struct randomMatrixWithImagEivals<MatrixType, 0>
}
}
MatrixType A = MatrixType::Random(size, size);
return A.inverse() * diag * A;
HouseholderQR<MatrixType> QRofA(A);
return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
}
};
@ -88,10 +101,12 @@ struct randomMatrixWithImagEivals<MatrixType, 1>
+ ei_random<Scalar>() * Scalar(RealScalar(0.01));
}
MatrixType A = MatrixType::Random(size, size);
return A.inverse() * diag * A;
HouseholderQR<MatrixType> QRofA(A);
return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
}
};
template<typename MatrixType>
void testMatrixExponential(const MatrixType& A)
{
@ -99,50 +114,45 @@ void testMatrixExponential(const MatrixType& A)
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef std::complex<RealScalar> ComplexScalar;
for (int i = 0; i < g_repeat; i++) {
VERIFY_IS_APPROX(ei_matrix_exponential(A),
ei_matrix_function(A, StdStemFunctions<ComplexScalar>::exp));
}
VERIFY_IS_APPROX(ei_matrix_exponential(A),
ei_matrix_function(A, StdStemFunctions<ComplexScalar>::exp));
}
template<typename MatrixType>
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<typename MatrixType>
void testGonioFunctions(const MatrixType& A)
{
typedef ei_traits<MatrixType> Traits;
typedef typename Traits::Scalar Scalar;
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef std::complex<RealScalar> ComplexScalar;
typedef Matrix<ComplexScalar, Traits::RowsAtCompileTime,
Traits::ColsAtCompileTime, MatrixType::Options> ComplexMatrix;
typedef Matrix<ComplexScalar, MatrixType::RowsAtCompileTime,
MatrixType::ColsAtCompileTime, MatrixType::Options> ComplexMatrix;
ComplexScalar imagUnit(0,1);
ComplexScalar two(2,0);
for (int i = 0; i < g_repeat; i++) {
ComplexMatrix Ac = A.template cast<ComplexScalar>();
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<ComplexScalar>();
VERIFY_IS_APPROX(sinAc, (exp_iA - exp_miA) / (two*imagUnit));
MatrixType cosA = ei_matrix_cos(A);
ComplexMatrix cosAc = cosA.template cast<ComplexScalar>();
VERIFY_IS_APPROX(cosAc, (exp_iA + exp_miA) / 2);
}
ComplexMatrix Ac = A.template cast<ComplexScalar>();
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<ComplexScalar>();
VERIFY_IS_APPROX_ABS(sinAc, (exp_iA - exp_miA) / (two*imagUnit));
MatrixType cosA = ei_matrix_cos(A);
ComplexMatrix cosAc = cosA.template cast<ComplexScalar>();
VERIFY_IS_APPROX_ABS(cosAc, (exp_iA + exp_miA) / 2);
}
template<typename MatrixType>