From 8c8ab9ae103930e21da76ed6dcbaa585e8af7ff4 Mon Sep 17 00:00:00 2001
From: Jitse Niesen <jitse@maths.leeds.ac.uk>
Date: Tue, 7 Jun 2011 14:44:43 +0100
Subject: Implement matrix logarithm + test + docs. Currently, test
 matrix_function_1 fails due to bug #288.

---
 unsupported/test/matrix_function.cpp | 21 +++++++++++++++++++++
 1 file changed, 21 insertions(+)

(limited to 'unsupported/test/matrix_function.cpp')

diff --git a/unsupported/test/matrix_function.cpp b/unsupported/test/matrix_function.cpp
index 04167abfb..c2ca5d5f1 100644
--- a/unsupported/test/matrix_function.cpp
+++ b/unsupported/test/matrix_function.cpp
@@ -120,6 +120,26 @@ void testMatrixExponential(const MatrixType& A)
   VERIFY_IS_APPROX(A.exp(), A.matrixFunction(StdStemFunctions<ComplexScalar>::exp));
 }
 
+template<typename MatrixType>
+void testMatrixLogarithm(const MatrixType& A)
+{
+  typedef typename internal::traits<MatrixType>::Scalar Scalar;
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+  typedef std::complex<RealScalar> ComplexScalar;
+
+  MatrixType scaledA;
+  RealScalar maxImagPartOfSpectrum = A.eigenvalues().imag().cwiseAbs().maxCoeff();
+  if (maxImagPartOfSpectrum >= 0.9 * M_PI)
+    scaledA = A * 0.9 * M_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<typename MatrixType>
 void testHyperbolicFunctions(const MatrixType& A)
 {
@@ -157,6 +177,7 @@ template<typename MatrixType>
 void testMatrix(const MatrixType& A)
 {
   testMatrixExponential(A);
+  testMatrixLogarithm(A);
   testHyperbolicFunctions(A);
   testGonioFunctions(A);
 }
-- 
cgit v1.2.3