From bbc02cecda1ba7efc60fc2ec603db912805b20e6 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 13 Mar 2023 18:42:24 +0100 Subject: [PATCH] Add test case for complex eigenvalues and eigenvectors --- tests/lapack/full_matrix_39.cc | 108 +++++++++++++++++++++++++++++ tests/lapack/full_matrix_39.output | 45 ++++++++++++ 2 files changed, 153 insertions(+) create mode 100644 tests/lapack/full_matrix_39.cc create mode 100644 tests/lapack/full_matrix_39.output diff --git a/tests/lapack/full_matrix_39.cc b/tests/lapack/full_matrix_39.cc new file mode 100644 index 0000000000..ff2a37783c --- /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 +#include +#include + +#include + +#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> 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 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> &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 index 0000000000..5773329d9a --- /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:: -- 2.39.5