You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

193 lines
6.5 KiB

  1. // This file is part of Eigen, a lightweight C++ template library
  2. // for linear algebra.
  3. //
  4. // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
  5. //
  6. // This Source Code Form is subject to the terms of the Mozilla
  7. // Public License v. 2.0. If a copy of the MPL was not distributed
  8. // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
  9. #include "main.h"
  10. #include <unsupported/StormEigen/MatrixFunctions>
  11. // Variant of VERIFY_IS_APPROX which uses absolute error instead of
  12. // relative error.
  13. #define VERIFY_IS_APPROX_ABS(a, b) VERIFY(test_isApprox_abs(a, b))
  14. template<typename Type1, typename Type2>
  15. inline bool test_isApprox_abs(const Type1& a, const Type2& b)
  16. {
  17. return ((a-b).array().abs() < test_precision<typename Type1::RealScalar>()).all();
  18. }
  19. // Returns a matrix with eigenvalues clustered around 0, 1 and 2.
  20. template<typename MatrixType>
  21. MatrixType randomMatrixWithRealEivals(const typename MatrixType::Index size)
  22. {
  23. typedef typename MatrixType::Index Index;
  24. typedef typename MatrixType::Scalar Scalar;
  25. typedef typename MatrixType::RealScalar RealScalar;
  26. MatrixType diag = MatrixType::Zero(size, size);
  27. for (Index i = 0; i < size; ++i) {
  28. diag(i, i) = Scalar(RealScalar(internal::random<int>(0,2)))
  29. + internal::random<Scalar>() * Scalar(RealScalar(0.01));
  30. }
  31. MatrixType A = MatrixType::Random(size, size);
  32. HouseholderQR<MatrixType> QRofA(A);
  33. return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
  34. }
  35. template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
  36. struct randomMatrixWithImagEivals
  37. {
  38. // Returns a matrix with eigenvalues clustered around 0 and +/- i.
  39. static MatrixType run(const typename MatrixType::Index size);
  40. };
  41. // Partial specialization for real matrices
  42. template<typename MatrixType>
  43. struct randomMatrixWithImagEivals<MatrixType, 0>
  44. {
  45. static MatrixType run(const typename MatrixType::Index size)
  46. {
  47. typedef typename MatrixType::Index Index;
  48. typedef typename MatrixType::Scalar Scalar;
  49. MatrixType diag = MatrixType::Zero(size, size);
  50. Index i = 0;
  51. while (i < size) {
  52. Index randomInt = internal::random<Index>(-1, 1);
  53. if (randomInt == 0 || i == size-1) {
  54. diag(i, i) = internal::random<Scalar>() * Scalar(0.01);
  55. ++i;
  56. } else {
  57. Scalar alpha = Scalar(randomInt) + internal::random<Scalar>() * Scalar(0.01);
  58. diag(i, i+1) = alpha;
  59. diag(i+1, i) = -alpha;
  60. i += 2;
  61. }
  62. }
  63. MatrixType A = MatrixType::Random(size, size);
  64. HouseholderQR<MatrixType> QRofA(A);
  65. return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
  66. }
  67. };
  68. // Partial specialization for complex matrices
  69. template<typename MatrixType>
  70. struct randomMatrixWithImagEivals<MatrixType, 1>
  71. {
  72. static MatrixType run(const typename MatrixType::Index size)
  73. {
  74. typedef typename MatrixType::Index Index;
  75. typedef typename MatrixType::Scalar Scalar;
  76. typedef typename MatrixType::RealScalar RealScalar;
  77. const Scalar imagUnit(0, 1);
  78. MatrixType diag = MatrixType::Zero(size, size);
  79. for (Index i = 0; i < size; ++i) {
  80. diag(i, i) = Scalar(RealScalar(internal::random<Index>(-1, 1))) * imagUnit
  81. + internal::random<Scalar>() * Scalar(RealScalar(0.01));
  82. }
  83. MatrixType A = MatrixType::Random(size, size);
  84. HouseholderQR<MatrixType> QRofA(A);
  85. return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
  86. }
  87. };
  88. template<typename MatrixType>
  89. void testMatrixExponential(const MatrixType& A)
  90. {
  91. typedef typename internal::traits<MatrixType>::Scalar Scalar;
  92. typedef typename NumTraits<Scalar>::Real RealScalar;
  93. typedef std::complex<RealScalar> ComplexScalar;
  94. VERIFY_IS_APPROX(A.exp(), A.matrixFunction(internal::stem_function_exp<ComplexScalar>));
  95. }
  96. template<typename MatrixType>
  97. void testMatrixLogarithm(const MatrixType& A)
  98. {
  99. typedef typename internal::traits<MatrixType>::Scalar Scalar;
  100. typedef typename NumTraits<Scalar>::Real RealScalar;
  101. MatrixType scaledA;
  102. RealScalar maxImagPartOfSpectrum = A.eigenvalues().imag().cwiseAbs().maxCoeff();
  103. if (maxImagPartOfSpectrum >= 0.9 * M_PI)
  104. scaledA = A * 0.9 * M_PI / maxImagPartOfSpectrum;
  105. else
  106. scaledA = A;
  107. // identity X.exp().log() = X only holds if Im(lambda) < pi for all eigenvalues of X
  108. MatrixType expA = scaledA.exp();
  109. MatrixType logExpA = expA.log();
  110. VERIFY_IS_APPROX(logExpA, scaledA);
  111. }
  112. template<typename MatrixType>
  113. void testHyperbolicFunctions(const MatrixType& A)
  114. {
  115. // Need to use absolute error because of possible cancellation when
  116. // adding/subtracting expA and expmA.
  117. VERIFY_IS_APPROX_ABS(A.sinh(), (A.exp() - (-A).exp()) / 2);
  118. VERIFY_IS_APPROX_ABS(A.cosh(), (A.exp() + (-A).exp()) / 2);
  119. }
  120. template<typename MatrixType>
  121. void testGonioFunctions(const MatrixType& A)
  122. {
  123. typedef typename MatrixType::Scalar Scalar;
  124. typedef typename NumTraits<Scalar>::Real RealScalar;
  125. typedef std::complex<RealScalar> ComplexScalar;
  126. typedef Matrix<ComplexScalar, MatrixType::RowsAtCompileTime,
  127. MatrixType::ColsAtCompileTime, MatrixType::Options> ComplexMatrix;
  128. ComplexScalar imagUnit(0,1);
  129. ComplexScalar two(2,0);
  130. ComplexMatrix Ac = A.template cast<ComplexScalar>();
  131. ComplexMatrix exp_iA = (imagUnit * Ac).exp();
  132. ComplexMatrix exp_miA = (-imagUnit * Ac).exp();
  133. ComplexMatrix sinAc = A.sin().template cast<ComplexScalar>();
  134. VERIFY_IS_APPROX_ABS(sinAc, (exp_iA - exp_miA) / (two*imagUnit));
  135. ComplexMatrix cosAc = A.cos().template cast<ComplexScalar>();
  136. VERIFY_IS_APPROX_ABS(cosAc, (exp_iA + exp_miA) / 2);
  137. }
  138. template<typename MatrixType>
  139. void testMatrix(const MatrixType& A)
  140. {
  141. testMatrixExponential(A);
  142. testMatrixLogarithm(A);
  143. testHyperbolicFunctions(A);
  144. testGonioFunctions(A);
  145. }
  146. template<typename MatrixType>
  147. void testMatrixType(const MatrixType& m)
  148. {
  149. // Matrices with clustered eigenvalue lead to different code paths
  150. // in MatrixFunction.h and are thus useful for testing.
  151. typedef typename MatrixType::Index Index;
  152. const Index size = m.rows();
  153. for (int i = 0; i < g_repeat; i++) {
  154. testMatrix(MatrixType::Random(size, size).eval());
  155. testMatrix(randomMatrixWithRealEivals<MatrixType>(size));
  156. testMatrix(randomMatrixWithImagEivals<MatrixType>::run(size));
  157. }
  158. }
  159. void test_matrix_function()
  160. {
  161. CALL_SUBTEST_1(testMatrixType(Matrix<float,1,1>()));
  162. CALL_SUBTEST_2(testMatrixType(Matrix3cf()));
  163. CALL_SUBTEST_3(testMatrixType(MatrixXf(8,8)));
  164. CALL_SUBTEST_4(testMatrixType(Matrix2d()));
  165. CALL_SUBTEST_5(testMatrixType(Matrix<double,5,5,RowMajor>()));
  166. CALL_SUBTEST_6(testMatrixType(Matrix4cd()));
  167. CALL_SUBTEST_7(testMatrixType(MatrixXd(13,13)));
  168. }