From a53f5e443341eaf6e8809cc30041128a50749583 Mon Sep 17 00:00:00 2001
From: Tobias Leibner <tobias.leibner@uni-muenster.de>
Date: Wed, 13 Dec 2017 13:38:52 +0100
Subject: [PATCH] [eigen-solver.internal.base] try to compute real eigenvectors
from complex ones
---
dune/xt/la/eigen-solver/internal/base.hh | 117 ++++++++++++++++++++---
1 file changed, 104 insertions(+), 13 deletions(-)
diff --git a/dune/xt/la/eigen-solver/internal/base.hh b/dune/xt/la/eigen-solver/internal/base.hh
index a41a1e107..23767b7b7 100644
--- a/dune/xt/la/eigen-solver/internal/base.hh
+++ b/dune/xt/la/eigen-solver/internal/base.hh
@@ -441,22 +441,113 @@ protected:
const size_t rows = CM::rows(*self.eigenvectors_);
const size_t cols = CM::cols(*self.eigenvectors_);
self.real_eigenvectors_ = std::make_unique<RealMatrixType>(RM::create(rows, cols));
- for (size_t ii = 0; ii < rows; ++ii)
+ bool is_complex = false;
+ for (size_t ii = 0; ii < rows; ++ii) {
for (size_t jj = 0; jj < cols; ++jj) {
const auto complex_value = CM::get_entry(*self.eigenvectors_, ii, jj);
- if (std::abs(complex_value.imag()) > tolerance)
- DUNE_THROW(Exceptions::eigen_solver_failed_bc_eigenvectors_are_not_real_as_requested,
- "These were the given options:\n\n"
- << self.options_
- << "\n\nThis was the given matrix: "
- << std::setprecision(17)
- << self.matrix_
- << "\nThese are the computed eigenvectors:\n\n"
- << std::setprecision(17)
- << *self.eigenvectors_);
+ if (std::abs(complex_value.imag()) > tolerance) {
+ is_complex = true;
+ ii = rows;
+ break;
+ }
RM::set_entry(*self.real_eigenvectors_, ii, jj, complex_value.real());
- }
- }
+ } // jj
+ } // ii
+
+ if (is_complex) {
+ try {
+ // try to get real eigenvectors from the complex ones. If both the matrix and the eigenvalues are real, the
+ // eigenvectors can also be chosen real. If there is a imaginary eigenvector, both real and imaginary part are
+ // eigenvectors (if non-zero) to the same eigenvalue. So to get real eigenvectors, sort the eigenvectors by
+ // eigenvalues and, separately for each eigenvalue, perform a Gram-Schmidt process with all real and imaginary
+ // parts of the eigenvectors
+ auto eigenvalues = self.real_eigenvalues();
+
+ // form groups of equal eigenvalues
+ struct Cmp
+ {
+ bool operator()(const RealType& a, const RealType& b) const
+ {
+ return XT::Common::FloatCmp::lt(a, b);
+ }
+ };
+ std::vector<std::vector<size_t>> eigenvalue_groups;
+ std::vector<size_t> eigenvalue_multiplicity;
+ std::set<RealType, Cmp> eigenvalues_done;
+ for (size_t jj = 0; jj < rows; ++jj) {
+ const auto curr_eigenvalue = eigenvalues[jj];
+ if (!eigenvalues_done.count(curr_eigenvalue)) {
+ std::vector<size_t> curr_group;
+ curr_group.push_back(jj);
+ eigenvalue_multiplicity.push_back(1);
+ for (size_t kk = jj + 1; kk < rows; ++kk) {
+ if (XT::Common::FloatCmp::eq(curr_eigenvalue, eigenvalues[kk])) {
+ curr_group.push_back(kk);
+ ++(eigenvalue_multiplicity.back());
+ }
+ } // kk
+ eigenvalue_groups.push_back(curr_group);
+ eigenvalues_done.insert(curr_eigenvalue);
+ }
+ } // jj
+
+ // For each eigenvalue, calculate a orthogonal basis of the n-dim real eigenspace from the 2n real
+ // and imaginary parts of the complex eigenvectors
+ for (size_t kk = 0; kk < eigenvalue_groups.size(); ++kk) {
+ const auto& group = eigenvalue_groups[kk];
+ typedef typename XT::LA::CommonDenseVector<RealType> RealVectorType;
+ std::vector<RealVectorType> input_vectors(2 * eigenvalue_multiplicity[kk], RealVectorType(rows, 0.));
+ size_t index = 0;
+ for (const auto& jj : group) {
+ for (size_t ll = 0; ll < cols; ++ll) {
+ input_vectors[index][ll] = CM::get_entry(*self.eigenvectors_, ll, jj).real();
+ input_vectors[index + 1][ll] = CM::get_entry(*self.eigenvectors_, ll, jj).imag();
+ }
+ index += 2;
+ } // jj
+
+ // orthonormalize
+ for (size_t ii = 1; ii < input_vectors.size(); ++ii) {
+ auto& v_i = input_vectors[ii];
+ for (size_t jj = 0; jj < ii; ++jj) {
+ const auto& v_j = input_vectors[jj];
+ const auto vj_vj = v_j.dot(v_j);
+ if (XT::Common::FloatCmp::eq(vj_vj, 0.))
+ continue;
+ const auto vj_vi = v_j.dot(v_i);
+ for (size_t rr = 0; rr < rows; ++rr)
+ v_i[rr] -= vj_vi / vj_vj * v_j[rr];
+ } // jj
+ auto l2_norm = std::accumulate(v_i.begin(), v_i.end(), 0.);
+ if (XT::Common::FloatCmp::ne(l2_norm, 0.))
+ v_i *= 1. / std::sqrt(l2_norm);
+ } // ii
+ // copy eigenvectors back to eigenvectors matrix
+ index = 0;
+ for (size_t ii = 0; ii < input_vectors.size(); ++ii) {
+ if (XT::Common::FloatCmp::ne(input_vectors[ii], RealVectorType(rows, 0.))) {
+ index++;
+ if (index >= eigenvalue_multiplicity[ii])
+ DUNE_THROW(Exceptions::eigen_solver_failed_bc_eigenvectors_are_not_real_as_requested,
+ "Eigenvectors are complex and calculating real eigenvectors failed!");
+ for (size_t rr = 0; rr < rows; ++rr)
+ RM::set_entry(*self.real_eigenvectors_, rr, group[index], input_vectors[ii].get_entry(rr));
+ } // if (input_vectors[ii] != 0)
+ } // ii
+ } // kk
+ } catch (const Dune::Exception& e) {
+ DUNE_THROW(Exceptions::eigen_solver_failed_bc_eigenvectors_are_not_real_as_requested,
+ "These were the given options:\n\n"
+ << self.options_
+ << "\n\nThis was the given matrix: "
+ << std::setprecision(17)
+ << self.matrix_
+ << "\nThese are the computed eigenvectors:\n\n"
+ << std::setprecision(17)
+ << *self.eigenvectors_);
+ } // try/catch
+ } // if(is_complex)
+ } // static void compute(...)
}; // real_eigenvectors_helper<true, ...>
template <class T>
--
GitLab