Add test case for complex eigenvalues and eigenvectors

diff --git a/tests/lapack/full_matrix_39.cc b/tests/lapack/full_matrix_39.cc
new file mode 100644 (file)
index 0000000..ff2a377
--- /dev/null
+++ b/tests/lapack/full_matrix_39.cc
@@ -0,0 +1,108 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Tests eigenvalues and eigenvectors of LAPACKFullMatrix with complex numbers
+
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+/*
+ * Eigenvalues and -vectors of this system are
+ * lambda = 3     v = (0.707,  0.707)
+ * lambda = 5     v = (0.707, -0.707)
+ */
+const double mat1[2][2] = {{4., -1.}, {-1., 4.}};
+
+/*
+ * Eigenvalues and -vectors of this system are
+ * lambda = 1+1i     v = (0.707, -0.707i)
+ * lambda = 1-1i     v = (0.707,  0.707i)
+ */
+const double mat2[2][2] = {{1., -1.}, {1., 1.}};
+
+
+void
+check_matrix(const double *matrix_pointer, const bool make_imaginary)
+{
+  LAPACKFullMatrix<std::complex<double>> LA(2, 2);
+  for (unsigned int i = 0; i < 2; ++i)
+    for (unsigned int j = 0; j < 2; ++j)
+      {
+        LA(i, j).real(matrix_pointer[i * 2 + j]);
+        if (make_imaginary)
+          LA(i, j).imag(1.0);
+      }
+
+  deallog << "Checking matrix " << std::endl;
+  LA.print_formatted(deallog.get_file_stream());
+
+  LA.compute_eigenvalues(true, true);
+  deallog << "Eigenvalues: ";
+  for (unsigned int i = 0; i < LA.m(); ++i)
+    {
+      std::complex<double> lambda = LA.eigenvalue(i);
+      deallog << lambda << " ";
+    }
+  deallog << std::endl;
+
+  // LAPACK might produce different signs in the eigenvectors depending on the
+  // implementation. To avoid these problems, print eigenvectors normalized by
+  // the sign of the first non-zero real part of the eigenvectors
+  const auto print_eigenvectors =
+    [](const FullMatrix<std::complex<double>> &eigenvectors) {
+      for (unsigned int col = 0; col < eigenvectors.n(); ++col)
+        {
+          double sign = 1.;
+          for (unsigned int row = 0; row < eigenvectors.n(); ++row)
+            if (std::abs(eigenvectors(row, col).real()) > 1e-12)
+              {
+                sign = (eigenvectors(row, col).real() > 0.) ? 1.0 : -1.0;
+                break;
+              }
+          for (unsigned int row = 0; row < eigenvectors.n(); ++row)
+            deallog << sign * eigenvectors(row, col) << "  ";
+          deallog << std::endl;
+        }
+    };
+  deallog << "Right eigenvectors" << std::endl;
+  print_eigenvectors(LA.get_right_eigenvectors());
+  deallog << "Left eigenvectors" << std::endl;
+  print_eigenvectors(LA.get_left_eigenvectors());
+  deallog << std::endl;
+}
+
+
+
+int
+main()
+{
+  initlog();
+
+  // Test symmetric system
+  check_matrix(&mat1[0][0], false);
+
+  // Test non-symmetric system with complex eigenvalues/eigenvectors
+  check_matrix(&mat2[0][0], false);
+
+  // Test with imaginary contribution in matrix
+  check_matrix(&mat1[0][0], true);
+  check_matrix(&mat2[0][0], true);
+}

diff --git a/tests/lapack/full_matrix_39.output b/tests/lapack/full_matrix_39.output
new file mode 100644 (file)
index 0000000..5773329
--- /dev/null
+++ b/tests/lapack/full_matrix_39.output
@@ -0,0 +1,45 @@
+
+DEAL::Checking matrix 
+(4.000e+00,0.000e+00) (-1.000e+00,0.000e+00) 
+(-1.000e+00,0.000e+00) (4.000e+00,0.000e+00) 
+DEAL::Eigenvalues: (5.00000,0.00000) (3.00000,0.00000) 
+DEAL::Right eigenvectors
+DEAL::(0.707107,0.00000)  (-0.707107,0.00000)  
+DEAL::(0.707107,0.00000)  (0.707107,0.00000)  
+DEAL::Left eigenvectors
+DEAL::(0.707107,0.00000)  (-0.707107,0.00000)  
+DEAL::(0.707107,0.00000)  (0.707107,0.00000)  
+DEAL::
+DEAL::Checking matrix 
+(1.000e+00,0.000e+00) (-1.000e+00,0.000e+00) 
+(1.000e+00,0.000e+00) (1.000e+00,0.000e+00) 
+DEAL::Eigenvalues: (1.00000,1.00000) (1.00000,-1.00000) 
+DEAL::Right eigenvectors
+DEAL::(0.707107,0.00000)  (0.00000,-0.707107)  
+DEAL::(0.707107,0.00000)  (-8.32667e-17,0.707107)  
+DEAL::Left eigenvectors
+DEAL::(-1.52519e-16,0.707107)  (0.707107,0.00000)  
+DEAL::(-1.94289e-16,-0.707107)  (0.707107,0.00000)  
+DEAL::
+DEAL::Checking matrix 
+(4.000e+00,1.000e+00) (-1.000e+00,1.000e+00) 
+(-1.000e+00,1.000e+00) (4.000e+00,1.000e+00) 
+DEAL::Eigenvalues: (5.00000,-4.44089e-16) (3.00000,2.00000) 
+DEAL::Right eigenvectors
+DEAL::(0.707107,0.00000)  (-0.707107,2.22045e-16)  
+DEAL::(0.707107,0.00000)  (0.707107,6.93889e-17)  
+DEAL::Left eigenvectors
+DEAL::(0.707107,-2.22045e-16)  (-0.707107,0.00000)  
+DEAL::(0.707107,0.00000)  (0.707107,-3.60822e-16)  
+DEAL::
+DEAL::Checking matrix 
+(1.000e+00,1.000e+00) (-1.000e+00,1.000e+00) 
+(1.000e+00,1.000e+00) (1.000e+00,1.000e+00) 
+DEAL::Eigenvalues: (1.00000,2.41421) (1.00000,-0.414214) 
+DEAL::Right eigenvectors
+DEAL::(0.707107,0.00000)  (0.500000,-0.500000)  
+DEAL::(0.500000,0.500000)  (-0.707107,0.00000)  
+DEAL::Left eigenvectors
+DEAL::(0.707107,0.00000)  (0.500000,-0.500000)  
+DEAL::(0.500000,0.500000)  (-0.707107,0.00000)  
+DEAL::