From ea97ca0a6035419b862a110f69f71aa244abac7c Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 13 Mar 2023 18:42:03 +0100 Subject: [PATCH] Return eigenvectors as full matrix. Fix complex case --- include/deal.II/lac/lapack_full_matrix.h | 108 +++++++++-------- source/lac/lapack_full_matrix.cc | 140 ++++++++++++++++++++++- tests/lapack/full_matrix_38.cc | 11 +- tests/lapack/full_matrix_38.output | 18 ++- 4 files changed, 213 insertions(+), 64 deletions(-) diff --git a/include/deal.II/lac/lapack_full_matrix.h b/include/deal.II/lac/lapack_full_matrix.h index 6d91586177..78b1ad807e 100644 --- a/include/deal.II/lac/lapack_full_matrix.h +++ b/include/deal.II/lac/lapack_full_matrix.h @@ -75,14 +75,12 @@ public: */ explicit LAPACKFullMatrix(const size_type size = 0); - /** * Constructor. Initialize the matrix as a rectangular matrix $\rm{rows} * \times \rm{cols}$. */ LAPACKFullMatrix(const size_type rows, const size_type cols); - /** * Copy constructor. This constructor does a deep copy of the matrix. * Therefore, it poses a possible efficiency problem, if for example, @@ -481,8 +479,9 @@ public: * @note It is assumed that @p A, @p B and @p V have compatible sizes and that * @p C already has the right size. * - * @note This function is not provided by LAPACK. The function first forms $\rm{diag}(\mathbf V) \cdot \mathbf B$ product and - * then uses Xgemm function. + * @note This function is not provided by LAPACK. The function first forms + * $\rm{diag}(\mathbf V) \cdot \mathbf B$ product and then uses the Xgemm + * function. */ void Tmmult(LAPACKFullMatrix & C, @@ -830,30 +829,32 @@ public: compute_inverse_svd_with_kernel(const unsigned int kernel_size); /** - * Retrieve eigenvalue after compute_eigenvalues() was called. + * Retrieve eigenvalue after compute_eigenvalues() was called. The return + * type is always a complex number: In case the underlying matrix is + * real-valued, the type `std::complex` is returned, whereas + * `number` is returned for complex numbers. */ - std::complex + std::complex::real_type> eigenvalue(const size_type i) const; /** - * After a call to compute_eigenvalues(), this function returns all (right) - * eigenvectors as returned by LAPACK. This means that eigenvectors are - * contained in column-major ordering in a matrix associated to the given - * flat vector. Note that for real-valued matrices, there might appear - * complex eigenvalues with complex-conjugate values and eigenvectors. For - * those cases, LAPACK places a column of n() entries with the real part and - * the next column (corresponding to the complex conjugate eigenvalue) for - * the imaginary part of the eigenvector. + * After a call to compute_eigenvalues(), this function returns the $n\times + * n$ matrix of (right) eigenvectors in a decomposition of the form $A V = V + * \Lambda$. Note that this function constructs the associated matrix on the + * fly, since LAPACK packs complex-conjugate eigenvalue/eigenvector pairs of + * real-valued matrices into a real-valued return matrix. This call only + * succeeds in case the respective flag `right_eigenvectors` in + * compute_eigenvalues() has been set to true. */ - std::vector + FullMatrix::real_type>> get_right_eigenvectors() const; /** * Return the matrix of left eigenvectors after a call to - * compute_eigenvalues(), following the same convention as + * compute_eigenvalues(), following the same principles as * get_right_eigenvectors(). */ - std::vector + FullMatrix::real_type>> get_left_eigenvectors() const; /** @@ -1037,6 +1038,7 @@ LAPACKFullMatrix::set(const size_type i, } + template inline typename LAPACKFullMatrix::size_type LAPACKFullMatrix::m() const @@ -1044,6 +1046,8 @@ LAPACKFullMatrix::m() const return static_cast(this->n_rows()); } + + template inline typename LAPACKFullMatrix::size_type LAPACKFullMatrix::n() const @@ -1051,6 +1055,8 @@ LAPACKFullMatrix::n() const return static_cast(this->n_cols()); } + + template template inline void @@ -1109,6 +1115,7 @@ LAPACKFullMatrix::fill(const MatrixType &M, } + template template void @@ -1122,6 +1129,7 @@ LAPACKFullMatrix::vmult(Vector &, } + template template void @@ -1134,6 +1142,7 @@ LAPACKFullMatrix::vmult_add(Vector &, } + template template void @@ -1147,6 +1156,7 @@ LAPACKFullMatrix::Tvmult(Vector &, } + template template void @@ -1154,54 +1164,50 @@ LAPACKFullMatrix::Tvmult_add(Vector &, const Vector &) const { Assert(false, - ExcMessage( - "LAPACKFullMatrix::Tvmult_add must be called with a " - "matching Vector vector type.")); + ExcMessage("LAPACKFullMatrix::Tvmult_add must be called " + "with a matching Vector vector type.")); } -template -inline std::complex -LAPACKFullMatrix::eigenvalue(const size_type i) const -{ - Assert(state & LAPACKSupport::eigenvalues, ExcInvalidState()); - Assert(wr.size() == this->n_rows(), ExcInternalError()); - Assert(wi.size() == this->n_rows(), ExcInternalError()); - AssertIndexRange(i, this->n_rows()); - - if (numbers::NumberTraits::is_complex) - return std::complex(wi[i]); - else - return std::complex(wr[i], wi[i]); -} - -template -inline std::vector -LAPACKFullMatrix::get_right_eigenvectors() const +namespace internal { - Assert(state & LAPACKSupport::eigenvalues, ExcInvalidState()); - Assert(vr.size() == this->n_rows() * this->n_cols(), - ExcMessage("Right eigenvectors are not available! Did you " - "set the associated flag in compute_eigenvalues()")); + namespace LAPACKFullMatrixImplementation + { + template + std::complex + pack_complex(const RealNumber &real_part, const RealNumber &imaginary_part) + { + return std::complex(real_part, imaginary_part); + } + + // The eigenvalues in LAPACKFullMatrix with complex-valued matrices are + // contained in the 'wi' array, ignoring the 'wr' array. + template + std::complex + pack_complex(const Number &, const std::complex &complex_number) + { + return complex_number; + } + } // namespace LAPACKFullMatrixImplementation +} // namespace internal - return vr; -} template -inline std::vector -LAPACKFullMatrix::get_left_eigenvectors() const +inline std::complex::real_type> +LAPACKFullMatrix::eigenvalue(const size_type i) const { Assert(state & LAPACKSupport::eigenvalues, ExcInvalidState()); - Assert(vl.size() == this->n_rows() * this->n_cols(), - ExcMessage("Left eigenvectors are not available! Did you " - "set the associated flag in compute_eigenvalues()")); + Assert(wr.size() == this->n_rows(), ExcInternalError()); + Assert(wi.size() == this->n_rows(), ExcInternalError()); + AssertIndexRange(i, this->n_rows()); - return vl; + return internal::LAPACKFullMatrixImplementation::pack_complex(wr[i], wi[i]); } + template inline number LAPACKFullMatrix::singular_value(const size_type i) const @@ -1214,6 +1220,7 @@ LAPACKFullMatrix::singular_value(const size_type i) const } + template inline const LAPACKFullMatrix & LAPACKFullMatrix::get_svd_u() const @@ -1225,6 +1232,7 @@ LAPACKFullMatrix::get_svd_u() const } + template inline const LAPACKFullMatrix & LAPACKFullMatrix::get_svd_vt() const diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc index 765aea4529..1d61dea4cb 100644 --- a/source/lac/lapack_full_matrix.cc +++ b/source/lac/lapack_full_matrix.cc @@ -95,6 +95,8 @@ namespace internal &info); } + + template void geev_helper(const char vl, @@ -125,10 +127,10 @@ namespace internal Assert(static_cast(n_rows * n_rows) <= right_eigenvectors.size(), ExcInternalError()); - Assert(std::max(1, work_flag) <= real_work.size(), - ExcInternalError()); Assert(work_flag == -1 || - std::max(1, 2 * n_rows) <= (work_flag), + std::max(1, work_flag) <= real_work.size(), + ExcInternalError()); + Assert(work_flag == -1 || std::max(1, 2 * n_rows) <= work_flag, ExcInternalError()); geev(&vl, @@ -240,6 +242,8 @@ namespace internal } // namespace LAPACKFullMatrixImplementation } // namespace internal + + template LAPACKFullMatrix::LAPACKFullMatrix(const size_type n) : TransposeTable(n, n) @@ -248,6 +252,7 @@ LAPACKFullMatrix::LAPACKFullMatrix(const size_type n) {} + template LAPACKFullMatrix::LAPACKFullMatrix(const size_type m, const size_type n) : TransposeTable(m, n) @@ -256,6 +261,7 @@ LAPACKFullMatrix::LAPACKFullMatrix(const size_type m, const size_type n) {} + template LAPACKFullMatrix::LAPACKFullMatrix(const LAPACKFullMatrix &M) : TransposeTable(M) @@ -264,6 +270,7 @@ LAPACKFullMatrix::LAPACKFullMatrix(const LAPACKFullMatrix &M) {} + template LAPACKFullMatrix & LAPACKFullMatrix::operator=(const LAPACKFullMatrix &M) @@ -370,6 +377,7 @@ LAPACKFullMatrix::reinit(const size_type m, const size_type n) } + template template LAPACKFullMatrix & @@ -387,6 +395,7 @@ LAPACKFullMatrix::operator=(const FullMatrix &M) } + template template LAPACKFullMatrix & @@ -404,6 +413,7 @@ LAPACKFullMatrix::operator=(const SparseMatrix &M) } + template LAPACKFullMatrix & LAPACKFullMatrix::operator=(const number d) @@ -419,6 +429,7 @@ LAPACKFullMatrix::operator=(const number d) } + template LAPACKFullMatrix & LAPACKFullMatrix::operator*=(const number factor) @@ -447,6 +458,7 @@ LAPACKFullMatrix::operator*=(const number factor) } + template LAPACKFullMatrix & LAPACKFullMatrix::operator/=(const number factor) @@ -780,6 +792,7 @@ LAPACKFullMatrix::vmult(Vector & w, } + template void LAPACKFullMatrix::Tvmult(Vector & w, @@ -918,6 +931,7 @@ LAPACKFullMatrix::Tvmult(Vector & w, } + template void LAPACKFullMatrix::vmult_add(Vector & w, @@ -927,6 +941,7 @@ LAPACKFullMatrix::vmult_add(Vector & w, } + template void LAPACKFullMatrix::Tvmult_add(Vector & w, @@ -936,6 +951,7 @@ LAPACKFullMatrix::Tvmult_add(Vector & w, } + template void LAPACKFullMatrix::mmult(LAPACKFullMatrix & C, @@ -970,6 +986,7 @@ LAPACKFullMatrix::mmult(LAPACKFullMatrix & C, } + template void LAPACKFullMatrix::mmult(FullMatrix & C, @@ -1088,6 +1105,7 @@ LAPACKFullMatrix::transpose(LAPACKFullMatrix &B) const } + template void LAPACKFullMatrix::scale_rows(const Vector &V) @@ -1162,6 +1180,7 @@ LAPACKFullMatrix::Tmmult(LAPACKFullMatrix & C, } + template void LAPACKFullMatrix::Tmmult(FullMatrix & C, @@ -1197,6 +1216,7 @@ LAPACKFullMatrix::Tmmult(FullMatrix & C, } + template void LAPACKFullMatrix::mTmult(LAPACKFullMatrix & C, @@ -1290,6 +1310,7 @@ LAPACKFullMatrix::mTmult(FullMatrix & C, } + template void LAPACKFullMatrix::TmTmult(LAPACKFullMatrix & C, @@ -1324,6 +1345,7 @@ LAPACKFullMatrix::TmTmult(LAPACKFullMatrix & C, } + template void LAPACKFullMatrix::TmTmult(FullMatrix & C, @@ -1359,6 +1381,7 @@ LAPACKFullMatrix::TmTmult(FullMatrix & C, } + template void LAPACKFullMatrix::compute_lu_factorization() @@ -1647,6 +1670,7 @@ LAPACKFullMatrix::compute_svd() } + template void LAPACKFullMatrix::compute_inverse_svd(const double threshold) @@ -1876,6 +1900,7 @@ LAPACKFullMatrix::determinant() const } + template void LAPACKFullMatrix::compute_eigenvalues(const bool right, const bool left) @@ -1932,6 +1957,7 @@ LAPACKFullMatrix::compute_eigenvalues(const bool right, const bool left) // resize workspace array work.resize(lwork); + real_work.resize(lwork); // Finally compute the eigenvalues. internal::LAPACKFullMatrixImplementation::geev_helper(jobvl, @@ -1956,6 +1982,111 @@ LAPACKFullMatrix::compute_eigenvalues(const bool right, const bool left) } + +namespace +{ + // This function extracts complex eigenvectors from the underlying 'number' + // array 'vr' of the LAPACK eigenvalue routine. For real-valued matrices + // addressed by this function specialization, we might get complex + // eigenvalues, which come in complex-conjugate pairs. In LAPACK, a compact + // storage scheme is applied that stores the real and imaginary part of + // eigenvectors only once, putting the real parts in one column and the + // imaginary part in the next of a real-valued array. Here, we do the + // unpacking into the usual complex values. + template + void + unpack_lapack_eigenvector_and_increment_index( + const std::vector & vr, + const std::complex & eigenvalue, + FullMatrix> &result, + unsigned int & index) + { + const std::size_t n = result.n(); + if (eigenvalue.imag() != 0.) + { + for (std::size_t j = 0; j < n; ++j) + { + result(j, index).real(vr[index * n + j]); + result(j, index + 1).real(vr[index * n + j]); + result(j, index).imag(vr[(index + 1) * n + j]); + result(j, index + 1).imag(-vr[(index + 1) * n + j]); + } + + // we filled two columns with the complex-conjugate pair, so increment + // returned index by 2 + index += 2; + } + else + { + for (unsigned int j = 0; j < n; ++j) + result(j, index).real(vr[index * n + j]); + + // real-valued case, we only filled one column + ++index; + } + } + + // This specialization fills the eigenvectors for complex-valued matrices, + // in which case we simply read off the entry in the 'vr' array. + template + void + unpack_lapack_eigenvector_and_increment_index( + const std::vector &vr, + const ComplexNumber &, + FullMatrix &result, + unsigned int & index) + { + const std::size_t n = result.n(); + for (unsigned int j = 0; j < n; ++j) + result(j, index) = vr[index * n + j]; + + // complex-valued case always only fills a single column + ++index; + } +} // namespace + + + +template +FullMatrix::real_type>> +LAPACKFullMatrix::get_right_eigenvectors() const +{ + Assert(state & LAPACKSupport::eigenvalues, ExcInvalidState()); + Assert(vr.size() == this->n_rows() * this->n_cols(), + ExcMessage("Right eigenvectors are not available! Did you " + "set the associated flag in compute_eigenvalues()?")); + + FullMatrix::real_type>> + result(n(), n()); + + for (unsigned int i = 0; i < n();) + unpack_lapack_eigenvector_and_increment_index(vr, eigenvalue(i), result, i); + + return result; +} + + + +template +FullMatrix::real_type>> +LAPACKFullMatrix::get_left_eigenvectors() const +{ + Assert(state & LAPACKSupport::eigenvalues, ExcInvalidState()); + Assert(vl.size() == this->n_rows() * this->n_cols(), + ExcMessage("Left eigenvectors are not available! Did you " + "set the associated flag in compute_eigenvalues()?")); + + FullMatrix::real_type>> + result(n(), n()); + + for (unsigned int i = 0; i < n();) + unpack_lapack_eigenvector_and_increment_index(vl, eigenvalue(i), result, i); + + return result; +} + + + template void LAPACKFullMatrix::compute_eigenvalues_symmetric( @@ -2067,6 +2198,7 @@ LAPACKFullMatrix::compute_eigenvalues_symmetric( } + template void LAPACKFullMatrix::compute_generalized_eigenvalues_symmetric( @@ -2192,6 +2324,7 @@ LAPACKFullMatrix::compute_generalized_eigenvalues_symmetric( } + template void LAPACKFullMatrix::compute_generalized_eigenvalues_symmetric( @@ -2284,6 +2417,7 @@ LAPACKFullMatrix::compute_generalized_eigenvalues_symmetric( } + template void LAPACKFullMatrix::print_formatted(std::ostream & out, diff --git a/tests/lapack/full_matrix_38.cc b/tests/lapack/full_matrix_38.cc index 2475522972..f84e30fb55 100644 --- a/tests/lapack/full_matrix_38.cc +++ b/tests/lapack/full_matrix_38.cc @@ -57,13 +57,10 @@ check_matrix(const double *matrix_pointer) deallog << lambda << " "; } deallog << std::endl; - deallog << "Right eigenvectors "; - for (const auto d : LA.get_right_eigenvectors()) - deallog << d << " "; - deallog << std::endl; - deallog << "Left eigenvectors "; - for (const auto d : LA.get_left_eigenvectors()) - deallog << d << " "; + deallog << "Right eigenvectors" << std::endl; + LA.get_right_eigenvectors().print_formatted(deallog.get_file_stream()); + deallog << "Left eigenvectors" << std::endl; + LA.get_left_eigenvectors().print_formatted(deallog.get_file_stream()); deallog << std::endl; } diff --git a/tests/lapack/full_matrix_38.output b/tests/lapack/full_matrix_38.output index 1dd51a200a..28e5a4ca21 100644 --- a/tests/lapack/full_matrix_38.output +++ b/tests/lapack/full_matrix_38.output @@ -3,11 +3,21 @@ DEAL::Checking matrix 4.000e+00 -1.000e+00 -1.000e+00 4.000e+00 DEAL::Eigenvalues: (5.00000,0.00000) (3.00000,0.00000) -DEAL::Right eigenvectors 0.707107 -0.707107 0.707107 0.707107 -DEAL::Left eigenvectors 0.707107 -0.707107 0.707107 0.707107 +DEAL::Right eigenvectors +(7.071e-01,0.000e+00) (7.071e-01,0.000e+00) +(-7.071e-01,0.000e+00) (7.071e-01,0.000e+00) +DEAL::Left eigenvectors +(7.071e-01,0.000e+00) (7.071e-01,0.000e+00) +(-7.071e-01,0.000e+00) (7.071e-01,0.000e+00) +DEAL:: DEAL::Checking matrix 1.000e+00 -1.000e+00 1.000e+00 1.000e+00 DEAL::Eigenvalues: (1.00000,1.00000) (1.00000,-1.00000) -DEAL::Right eigenvectors 0.707107 0.00000 0.00000 -0.707107 -DEAL::Left eigenvectors -0.707107 0.00000 0.00000 0.707107 +DEAL::Right eigenvectors +(7.071e-01,0.000e+00) (7.071e-01,0.000e+00) +(0.000e+00,-7.071e-01) (0.000e+00,7.071e-01) +DEAL::Left eigenvectors +(-7.071e-01,0.000e+00) (-7.071e-01,0.000e+00) +(0.000e+00,7.071e-01) (0.000e+00,-7.071e-01) +DEAL:: -- 2.39.5