From f0238cfb6997c4acfc2bd200de7295f3fa36968f Mon Sep 17 00:00:00 2001 From: Stanislaw Halik Date: Sun, 3 Mar 2019 21:09:10 +0100 Subject: don't index Eigen --- eigen/unsupported/test/matrix_function.cpp | 189 ----------------------------- 1 file changed, 189 deletions(-) delete mode 100644 eigen/unsupported/test/matrix_function.cpp (limited to 'eigen/unsupported/test/matrix_function.cpp') diff --git a/eigen/unsupported/test/matrix_function.cpp b/eigen/unsupported/test/matrix_function.cpp deleted file mode 100644 index 6a2b219..0000000 --- a/eigen/unsupported/test/matrix_function.cpp +++ /dev/null @@ -1,189 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2010 Jitse Niesen -// -// 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 - -// 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 typename MatrixType::Index size) -{ - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; - MatrixType diag = MatrixType::Zero(size, size); - for (Index i = 0; i < size; ++i) { - diag(i, i) = Scalar(RealScalar(internal::random(0,2))) - + internal::random() * Scalar(RealScalar(0.01)); - } - MatrixType A = MatrixType::Random(size, size); - HouseholderQR QRofA(A); - return QRofA.householderQ().inverse() * diag * QRofA.householderQ(); -} - -template ::Scalar>::IsComplex> -struct randomMatrixWithImagEivals -{ - // Returns a matrix with eigenvalues clustered around 0 and +/- i. - static MatrixType run(const typename MatrixType::Index size); -}; - -// Partial specialization for real matrices -template -struct randomMatrixWithImagEivals -{ - static MatrixType run(const typename MatrixType::Index size) - { - typedef typename MatrixType::Scalar Scalar; - MatrixType diag = MatrixType::Zero(size, size); - Index i = 0; - while (i < size) { - Index randomInt = internal::random(-1, 1); - if (randomInt == 0 || i == size-1) { - diag(i, i) = internal::random() * Scalar(0.01); - ++i; - } else { - Scalar alpha = Scalar(randomInt) + internal::random() * Scalar(0.01); - diag(i, i+1) = alpha; - diag(i+1, i) = -alpha; - i += 2; - } - } - MatrixType A = MatrixType::Random(size, size); - HouseholderQR QRofA(A); - return QRofA.householderQ().inverse() * diag * QRofA.householderQ(); - } -}; - -// Partial specialization for complex matrices -template -struct randomMatrixWithImagEivals -{ - static MatrixType run(const typename MatrixType::Index size) - { - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; - const Scalar imagUnit(0, 1); - MatrixType diag = MatrixType::Zero(size, size); - for (Index i = 0; i < size; ++i) { - diag(i, i) = Scalar(RealScalar(internal::random(-1, 1))) * imagUnit - + internal::random() * Scalar(RealScalar(0.01)); - } - MatrixType A = MatrixType::Random(size, size); - HouseholderQR QRofA(A); - return QRofA.householderQ().inverse() * diag * QRofA.householderQ(); - } -}; - - -template -void testMatrixExponential(const MatrixType& A) -{ - typedef typename internal::traits::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; - typedef std::complex ComplexScalar; - - VERIFY_IS_APPROX(A.exp(), A.matrixFunction(internal::stem_function_exp)); -} - -template -void testMatrixLogarithm(const MatrixType& A) -{ - typedef typename internal::traits::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; - - MatrixType scaledA; - RealScalar maxImagPartOfSpectrum = A.eigenvalues().imag().cwiseAbs().maxCoeff(); - if (maxImagPartOfSpectrum >= RealScalar(0.9L * EIGEN_PI)) - scaledA = A * RealScalar(0.9L * EIGEN_PI) / maxImagPartOfSpectrum; - else - scaledA = A; - - // identity X.exp().log() = X only holds if Im(lambda) < pi for all eigenvalues of X - MatrixType expA = scaledA.exp(); - MatrixType logExpA = expA.log(); - VERIFY_IS_APPROX(logExpA, scaledA); -} - -template -void testHyperbolicFunctions(const MatrixType& A) -{ - // Need to use absolute error because of possible cancellation when - // adding/subtracting expA and expmA. - VERIFY_IS_APPROX_ABS(A.sinh(), (A.exp() - (-A).exp()) / 2); - VERIFY_IS_APPROX_ABS(A.cosh(), (A.exp() + (-A).exp()) / 2); -} - -template -void testGonioFunctions(const MatrixType& A) -{ - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; - typedef std::complex ComplexScalar; - typedef Matrix ComplexMatrix; - - ComplexScalar imagUnit(0,1); - ComplexScalar two(2,0); - - ComplexMatrix Ac = A.template cast(); - - ComplexMatrix exp_iA = (imagUnit * Ac).exp(); - ComplexMatrix exp_miA = (-imagUnit * Ac).exp(); - - ComplexMatrix sinAc = A.sin().template cast(); - VERIFY_IS_APPROX_ABS(sinAc, (exp_iA - exp_miA) / (two*imagUnit)); - - ComplexMatrix cosAc = A.cos().template cast(); - VERIFY_IS_APPROX_ABS(cosAc, (exp_iA + exp_miA) / 2); -} - -template -void testMatrix(const MatrixType& A) -{ - testMatrixExponential(A); - testMatrixLogarithm(A); - testHyperbolicFunctions(A); - testGonioFunctions(A); -} - -template -void testMatrixType(const MatrixType& m) -{ - // Matrices with clustered eigenvalue lead to different code paths - // in MatrixFunction.h and are thus useful for testing. - - const Index size = m.rows(); - for (int i = 0; i < g_repeat; i++) { - testMatrix(MatrixType::Random(size, size).eval()); - testMatrix(randomMatrixWithRealEivals(size)); - testMatrix(randomMatrixWithImagEivals::run(size)); - } -} - -void test_matrix_function() -{ - CALL_SUBTEST_1(testMatrixType(Matrix())); - CALL_SUBTEST_2(testMatrixType(Matrix3cf())); - CALL_SUBTEST_3(testMatrixType(MatrixXf(8,8))); - CALL_SUBTEST_4(testMatrixType(Matrix2d())); - CALL_SUBTEST_5(testMatrixType(Matrix())); - CALL_SUBTEST_6(testMatrixType(Matrix4cd())); - CALL_SUBTEST_7(testMatrixType(MatrixXd(13,13))); -} -- cgit v1.2.3