From b18f737aa103c7589f589b354676f476d644c631 Mon Sep 17 00:00:00 2001
From: Jitse Niesen <jitse@maths.leeds.ac.uk>
Date: Sat, 13 Feb 2010 22:49:01 +0000
Subject: Test matrix functions with matrices with clustered imaginary eivals.
 The idea is that these test MatrixFunction::swapEntriesInSchur(), which is
 not covered by existing tests. This did not work out as expected, but
 nevertheless it is a good test so I left it in.

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

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

diff --git a/unsupported/test/matrix_function.cpp b/unsupported/test/matrix_function.cpp
index d1f5824d8..de63937ad 100644
--- a/unsupported/test/matrix_function.cpp
+++ b/unsupported/test/matrix_function.cpp
@@ -40,6 +40,58 @@ MatrixType randomMatrixWithRealEivals(const int size)
   return A.inverse() * diag * A;
 }
 
+template <typename MatrixType, int IsComplex = NumTraits<typename ei_traits<MatrixType>::Scalar>::IsComplex>
+struct randomMatrixWithImagEivals
+{
+  // Returns a matrix with eigenvalues clustered around 0 and +/- i.
+  static MatrixType run(const int size);
+};
+
+// Partial specialization for real matrices
+template<typename MatrixType>
+struct randomMatrixWithImagEivals<MatrixType, 0>
+{
+  static MatrixType run(const int size)
+  {
+    typedef typename MatrixType::Scalar Scalar;
+    MatrixType diag = MatrixType::Zero(size, size);
+    int i = 0;
+    while (i < size) {
+      int randomInt = ei_random<int>(-1, 1);
+      if (randomInt == 0 || i == size-1) {
+        diag(i, i) = ei_random<Scalar>() * Scalar(0.01);
+        ++i;
+      } else {
+        Scalar alpha = Scalar(randomInt) + ei_random<Scalar>() * Scalar(0.01);
+        diag(i, i+1) = alpha;
+        diag(i+1, i) = -alpha;
+        i += 2;
+      }
+    }
+    MatrixType A = MatrixType::Random(size, size);
+    return A.inverse() * diag * A;
+  }
+};
+
+// Partial specialization for complex matrices
+template<typename MatrixType>
+struct randomMatrixWithImagEivals<MatrixType, 1>
+{
+  static MatrixType run(const int size)
+  {
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename MatrixType::RealScalar RealScalar;
+    const Scalar imagUnit(0, 1);
+    MatrixType diag = MatrixType::Zero(size, size);
+    for (int i = 0; i < size; ++i) {
+      diag(i, i) = Scalar(RealScalar(ei_random<int>(-1, 1))) * imagUnit
+        + ei_random<Scalar>() * Scalar(RealScalar(0.01));
+    }
+    MatrixType A = MatrixType::Random(size, size);
+    return A.inverse() * diag * A;
+  }
+};
+
 template<typename MatrixType>
 void testMatrixExponential(const MatrixType& A)
 {
@@ -117,6 +169,7 @@ void testMatrixType(const MatrixType& m)
   for (int i = 0; i < g_repeat; i++) {
     testMatrix(MatrixType::Random(size, size).eval());
     testMatrix(randomMatrixWithRealEivals<MatrixType>(size));
+    testMatrix(randomMatrixWithImagEivals<MatrixType>::run(size));
   }
 }
 
-- 
cgit v1.2.3