diff --git a/dune/xt/la/eigen-solver/internal/base.hh b/dune/xt/la/eigen-solver/internal/base.hh
index e55b8fca9c080ffdc827175aea7d6c05bfb499c2..59a15978745300cb985a0feaea0a0283c999c0ce 100644
--- a/dune/xt/la/eigen-solver/internal/base.hh
+++ b/dune/xt/la/eigen-solver/internal/base.hh
@@ -50,7 +50,7 @@ static inline void ensure_eigen_solver_type(const std::string& type, const std::
if (!contained)
DUNE_THROW(Exceptions::eigen_solver_failed_bc_it_was_not_set_up_correctly,
"Given type '" << type << "' is not one of the available types: " << available_types);
-} // ... ensure_type(...)
+} // ... ensure_eigen_solver_type(...)
static inline Common::Configuration default_eigen_solver_options()
@@ -245,7 +245,7 @@ public:
*this, check_eigendecomposition > 0 ? check_eigendecomposition : options_->get<double>("real_tolerance"));
}
return *eigenvectors_inverse_;
- } // ... eigenvectors(...)
+ } // ... eigenvectors_inverse(...)
const RealMatrixType& real_eigenvectors() const
{
@@ -289,7 +289,7 @@ public:
assert_real_eigendecomposition);
}
return *real_eigenvectors_inverse_;
- } // ... eigenvectors(...)
+ } // ... real_eigenvectors_inverse(...)
protected:
void compute_and_check() const
diff --git a/dune/xt/la/exceptions.hh b/dune/xt/la/exceptions.hh
index 7a9d2d5c7cfec50dda2ab8cb68752048b1f9afc2..368c63d249b165b38e12cbbc8cb0b60b0a9178e8 100644
--- a/dune/xt/la/exceptions.hh
+++ b/dune/xt/la/exceptions.hh
@@ -67,6 +67,36 @@ class eigen_solver_failed_bc_result_is_not_an_eigendecomposition : public eigen_
{};
+class generalized_eigen_solver_failed : public Dune::Exception
+{};
+
+class generalized_eigen_solver_failed_bc_data_did_not_fulfill_requirements : public generalized_eigen_solver_failed
+{};
+
+class generalized_eigen_solver_failed_bc_it_was_not_set_up_correctly : public generalized_eigen_solver_failed
+{};
+
+class generalized_eigen_solver_failed_bc_result_contained_inf_or_nan : public generalized_eigen_solver_failed
+{};
+
+class generalized_eigen_solver_failed_bc_eigenvalues_are_not_real_as_requested : public generalized_eigen_solver_failed
+{};
+
+class generalized_eigen_solver_failed_bc_eigenvalues_are_not_positive_as_requested
+ : public generalized_eigen_solver_failed
+{};
+
+class generalized_eigen_solver_failed_bc_eigenvalues_are_not_negative_as_requested
+ : public generalized_eigen_solver_failed
+{};
+
+class generalized_eigen_solver_failed_bc_eigenvectors_are_not_real_as_requested : public generalized_eigen_solver_failed
+{};
+
+// class generalized_eigen_solver_failed_bc_result_is_not_an_eigendecomposition : public generalized_eigen_solver_failed
+//{};
+
+
class matrix_invert_failed : public Dune::Exception
{};
diff --git a/dune/xt/la/generalized-eigen-solver.hh b/dune/xt/la/generalized-eigen-solver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8e049388ba8117edba54574e2c7171c062ec3f45
--- /dev/null
+++ b/dune/xt/la/generalized-eigen-solver.hh
@@ -0,0 +1,125 @@
+// This file is part of the dune-xt-la project:
+// https://github.com/dune-community/dune-xt-la
+// Copyright 2009-2018 dune-xt-la developers and contributors. All rights reserved.
+// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+// or GPL-2.0+ (http://opensource.org/licenses/gpl-license)
+// with "runtime exception" (http://www.dune-project.org/license.html)
+// Authors:
+// Felix Schindler (2019)
+
+#ifndef DUNE_XT_LA_GENERALIZED_EIGEN_SOLVER_HH
+#define DUNE_XT_LA_GENERALIZED_EIGEN_SOLVER_HH
+
+#include <string>
+#include <vector>
+
+#include <dune/common/typetraits.hh>
+
+#include <dune/xt/common/configuration.hh>
+#include <dune/xt/common/type_traits.hh>
+
+namespace Dune {
+namespace XT {
+namespace LA {
+
+
+/**
+ * \brief A means to obtain available options at compile time.
+ * \note This class needs to be specialized for each MatrixType, the purpose of this variant is merely to document the
+ * expected functionality.
+ */
+template <class MatrixType, bool is_matrix = XT::Common::is_matrix<MatrixType>::value>
+class GeneralizedEigenSolverOptions
+{
+ static_assert(AlwaysFalse<MatrixType>::value,
+ "Please implement for given MatrixType and add the respective include below!");
+
+ static std::vector<std::string> types();
+
+ static Common::Configuration options(const std::string /*type*/ = "");
+}; // class GeneralizedEigenSolverOptions
+
+
+template <class MatrixType>
+std::vector<std::string> generalized_eigen_solver_types(const MatrixType& /*matrix*/)
+{
+ return GeneralizedEigenSolverOptions<MatrixType>::types();
+}
+
+
+template <class MatrixType>
+Common::Configuration generalized_eigen_solver_options(const MatrixType& /*matrix*/, const std::string type = "")
+{
+ return GeneralizedEigenSolverOptions<MatrixType>::options(type);
+}
+
+
+template <class MatrixImp, bool is_matrix = XT::Common::is_matrix<MatrixImp>::value>
+class GeneralizedEigenSolver
+{
+ static_assert(AlwaysFalse<MatrixImp>::value,
+ "Please implement for given MatrixType and add the respective include below!");
+
+public:
+ typedef MatrixImp MatrixType;
+ typedef double FieldType;
+ typedef int RealMatrixType;
+ typedef int ComplexMatrixType;
+
+ GeneralizedEigenSolver(const MatrixType& /*matrix*/, const std::string& /*type*/ = "")
+ {
+ static_assert(AlwaysFalse<MatrixType>::value,
+ "Please implement for given MatrixType and add the respective include below!");
+ }
+
+ GeneralizedEigenSolver(const MatrixType& /*matrix*/, const Common::Configuration /*opts*/)
+ {
+ static_assert(AlwaysFalse<MatrixType>::value,
+ "Please implement for given MatrixType and add the respective include below!");
+ }
+
+ const Common::Configuration& options() const;
+
+ const MatrixType& lhs_matrix() const;
+
+ const MatrixType& rhs_matrix() const;
+
+ const std::vector<Common::complex_t<FieldType>>& eigenvalues() const;
+
+ const std::vector<Common::real_t<FieldType>>& real_eigenvalues() const;
+
+ const std::vector<Common::real_t<FieldType>>&
+ min_eigenvalues(const size_t /*num_evs*/ = std::numeric_limits<size_t>::max()) const;
+
+ const std::vector<Common::real_t<FieldType>>&
+ max_eigenvalues(const size_t /*num_evs*/ = std::numeric_limits<size_t>::max()) const;
+
+ const ComplexMatrixType& eigenvectors() const;
+
+ const RealMatrixType& real_eigenvectors() const;
+}; // class GeneralizedEigenSolver
+
+
+template <class M>
+GeneralizedEigenSolver<M>
+make_generalized_eigen_solver(const M& lhs_matrix, const M& rhs_matrix, const std::string& type = "")
+{
+ return GeneralizedEigenSolver<M>(lhs_matrix, rhs_matrix, type);
+}
+
+
+template <class M>
+GeneralizedEigenSolver<M>
+make_generalized_eigen_solver(const M& lhs_matrix, const M& rhs_matrix, const XT::Common::Configuration& options)
+{
+ return GeneralizedEigenSolver<M>(lhs_matrix, rhs_matrix, options);
+}
+
+
+} // namespace LA
+} // namespace XT
+} // namespace Dune
+
+#include "generalized-eigen-solver/default.hh"
+
+#endif // DUNE_XT_LA_GENERALIZED_EIGEN_SOLVER_HH
diff --git a/dune/xt/la/generalized-eigen-solver/default.hh b/dune/xt/la/generalized-eigen-solver/default.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3f03043cf5b159f1943662a1a8c54f5331d90627
--- /dev/null
+++ b/dune/xt/la/generalized-eigen-solver/default.hh
@@ -0,0 +1,124 @@
+// This file is part of the dune-xt-la project:
+// https://github.com/dune-community/dune-xt-la
+// Copyright 2009-2018 dune-xt-la developers and contributors. All rights reserved.
+// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+// or GPL-2.0+ (http://opensource.org/licenses/gpl-license)
+// with "runtime exception" (http://www.dune-project.org/license.html)
+// Authors:
+// Felix Schindler (2019)
+
+#ifndef DUNE_XT_LA_GENERALIZED_EIGEN_SOLVER_DEFAULT_HH
+#define DUNE_XT_LA_GENERALIZED_EIGEN_SOLVER_DEFAULT_HH
+
+#include <vector>
+
+#include <dune/xt/common/matrix.hh>
+#include <dune/xt/la/container/conversion.hh>
+#include <dune/xt/la/exceptions.hh>
+#include <dune/xt/la/generalized-eigen-solver.hh>
+
+#include "internal/base.hh"
+#include "internal/lapacke.hh"
+
+namespace Dune {
+namespace XT {
+namespace LA {
+
+
+template <class MatrixType>
+class GeneralizedEigenSolverOptions<MatrixType, true>
+{
+public:
+ static std::vector<std::string> types()
+ {
+ std::vector<std::string> tps;
+ if (Common::Lapacke::available())
+ tps.push_back("lapack");
+ DUNE_THROW_IF(tps.empty(),
+ Exceptions::generalized_eigen_solver_failed,
+ "No backend available for generalized eigenvalue problems!");
+ return tps;
+ }
+
+ static Common::Configuration options(const std::string type = "")
+ {
+ const std::string actual_type = type.empty() ? types()[0] : type;
+ internal::ensure_generalized_eigen_solver_type(actual_type, types());
+ Common::Configuration opts = internal::default_generalized_eigen_solver_options();
+ opts["type"] = actual_type;
+ return opts;
+ }
+}; // class GeneralizedEigenSolverOptions<>
+
+
+template <class MatrixImp>
+class GeneralizedEigenSolver<MatrixImp, true>
+ : public internal::GeneralizedEigenSolverBase<
+ MatrixImp,
+ typename Common::MatrixAbstraction<MatrixImp>::ScalarType,
+ typename Common::MatrixAbstraction<MatrixImp>::template MatrixTypeTemplate<
+ Common::MatrixAbstraction<MatrixImp>::static_rows,
+ Common::MatrixAbstraction<MatrixImp>::static_cols,
+ typename Common::MatrixAbstraction<MatrixImp>::RealType>,
+ typename Common::MatrixAbstraction<MatrixImp>::template MatrixTypeTemplate<
+ Common::MatrixAbstraction<MatrixImp>::static_rows,
+ Common::MatrixAbstraction<MatrixImp>::static_cols,
+ std::complex<typename Common::MatrixAbstraction<MatrixImp>::RealType>>>
+{
+ using M = Common::MatrixAbstraction<MatrixImp>;
+ using BaseType = internal::GeneralizedEigenSolverBase<
+ MatrixImp,
+ typename M::ScalarType,
+ typename M::template MatrixTypeTemplate<M::static_rows, M::static_cols, typename M::RealType>,
+ typename M::template MatrixTypeTemplate<M::static_rows, M::static_cols, std::complex<typename M::RealType>>>;
+
+public:
+ using typename BaseType::ComplexMatrixType;
+ using typename BaseType::ComplexType;
+ using typename BaseType::MatrixType;
+ using typename BaseType::RealMatrixType;
+ using typename BaseType::RealType;
+ using RealM = Common::MatrixAbstraction<RealMatrixType>;
+ using ComplexM = Common::MatrixAbstraction<ComplexMatrixType>;
+
+ template <class... Args>
+ explicit GeneralizedEigenSolver(Args&&... args)
+ : BaseType(std::forward<Args>(args)...)
+ {}
+
+protected:
+ void compute() const override final
+ {
+ const auto type = options_->template get<std::string>("type");
+ if (type == "lapack") {
+ DUNE_THROW_IF(!Common::Lapacke::available(),
+ Exceptions::generalized_eigen_solver_failed_bc_it_was_not_set_up_correctly,
+ "Lapacke backend not available!");
+ if (!options_->template get<bool>("compute_eigenvectors"))
+ eigenvalues_ = std::make_unique<std::vector<ComplexType>>(
+ internal::compute_generalized_eigenvalues_using_lapack(lhs_matrix_, rhs_matrix_));
+ else {
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed,
+ "Eigenvectors of generalized eigenvalue problems are not yet implemented using lapacke!");
+ }
+ } else
+ DUNE_THROW(Common::Exceptions::internal_error,
+ "Given type '" << type << "' is none of GeneralizedEigenSolverOptions<...>::types(),"
+ << " and internal::GeneralizedEigenSolverBase promised to check this!"
+ << "\n\nThese are the available types:\n\n"
+ << GeneralizedEigenSolverOptions<MatrixType>::types());
+ } //... compute(...)
+
+ using BaseType::eigenvalues_;
+ using BaseType::eigenvectors_;
+ using BaseType::lhs_matrix_;
+ using BaseType::options_;
+ using BaseType::rhs_matrix_;
+}; // class GeneralizedEigenSolver<MatrixType, true>
+
+
+} // namespace LA
+} // namespace XT
+} // namespace Dune
+
+#endif // DUNE_XT_LA_GENERALIZED_EIGEN_SOLVER_DEFAULT_HH
diff --git a/dune/xt/la/generalized-eigen-solver/internal/base.hh b/dune/xt/la/generalized-eigen-solver/internal/base.hh
new file mode 100644
index 0000000000000000000000000000000000000000..fdd5e5c2fed027a0e5af509272f7e6581356f75b
--- /dev/null
+++ b/dune/xt/la/generalized-eigen-solver/internal/base.hh
@@ -0,0 +1,819 @@
+// This file is part of the dune-xt-la project:
+// https://github.com/dune-community/dune-xt-la
+// Copyright 2009-2018 dune-xt-la developers and contributors. All rights reserved.
+// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+// or GPL-2.0+ (http://opensource.org/licenses/gpl-license)
+// with "runtime exception" (http://www.dune-project.org/license.html)
+// Authors:
+// Felix Schindler (2019)
+
+#ifndef DUNE_XT_LA_GENERALIZED_EIGEN_SOLVER_INTERNAL_BASE_HH
+#define DUNE_XT_LA_GENERALIZED_EIGEN_SOLVER_INTERNAL_BASE_HH
+
+#include <algorithm>
+#include <functional>
+#include <memory>
+#include <numeric>
+
+#include <dune/xt/common/configuration.hh>
+#include <dune/xt/common/type_traits.hh>
+#include <dune/xt/common/vector.hh>
+#include <dune/xt/common/matrix.hh>
+
+#include <dune/xt/la/container/common/vector/dense.hh>
+#include <dune/xt/la/container/conversion.hh>
+#include <dune/xt/la/container/matrix-interface.hh>
+#include <dune/xt/la/exceptions.hh>
+#include <dune/xt/la/matrix-inverter.hh>
+
+namespace Dune {
+namespace XT {
+namespace LA {
+
+
+// forward
+template <class MatrixType, bool is_matrix>
+class GeneralizedEigenSolverOptions;
+
+
+namespace internal {
+
+
+static inline void ensure_generalized_eigen_solver_type(const std::string& type,
+ const std::vector<std::string>& available_types)
+{
+ bool contained = false;
+ for (const auto& tp : available_types)
+ if (type == tp)
+ contained = true;
+ if (!contained)
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_it_was_not_set_up_correctly,
+ "Given type '" << type << "' is not one of the available types: " << available_types);
+} // ... ensure_generalized_eigen_solver_type(...)
+
+
+static inline Common::Configuration default_generalized_eigen_solver_options()
+{
+ Common::Configuration opts;
+ opts["compute_eigenvalues"] = "true";
+ opts["compute_eigenvectors"] = "false";
+ opts["check_for_inf_nan"] = "true";
+ opts["real_tolerance"] = "1e-15"; // is only used if the respective assert_... is negative
+ opts["assert_real_eigenvalues"] = "-1"; // if positive, this is the check tolerance
+ opts["assert_positive_eigenvalues"] = "-1"; // if positive, this is the check tolerance
+ opts["assert_negative_eigenvalues"] = "-1"; // if positive, this is the check tolerance
+ opts["assert_real_eigenvectors"] = "-1"; // if positive, this is the check tolerance
+ return opts;
+} // ... default_generalized_eigen_solver_options(...)
+
+
+/**
+ * \sa default_generalized_eigen_solver_options()
+ * \note If the provided options contain a subtree "matrix-inverter" that one is forwarded on eigenvector inversion.
+ */
+template <class MatrixImp, class FieldImp, class RealMatrixImp, class ComplexMatrixImp>
+class GeneralizedEigenSolverBase
+{
+ static_assert(is_matrix<MatrixImp>::value || XT::Common::is_matrix<MatrixImp>::value, "");
+ static_assert(is_matrix<RealMatrixImp>::value || XT::Common::is_matrix<RealMatrixImp>::value, "");
+ static_assert(is_matrix<ComplexMatrixImp>::value || XT::Common::is_matrix<ComplexMatrixImp>::value, "");
+ static_assert((is_matrix<MatrixImp>::value && is_matrix<RealMatrixImp>::value && is_matrix<ComplexMatrixImp>::value)
+ || (XT::Common::is_matrix<MatrixImp>::value && XT::Common::is_matrix<RealMatrixImp>::value
+ && XT::Common::is_matrix<ComplexMatrixImp>::value),
+ "");
+ using ThisType = GeneralizedEigenSolverBase<MatrixImp, FieldImp, RealMatrixImp, ComplexMatrixImp>;
+
+public:
+ using MatrixType = MatrixImp;
+ using RealType = XT::Common::field_t<FieldImp>;
+ using ComplexType = XT::Common::complex_t<RealType>;
+ using RealMatrixType = RealMatrixImp;
+ using ComplexMatrixType = ComplexMatrixImp;
+
+ GeneralizedEigenSolverBase(const MatrixType& lhs_matrix, const MatrixType& rhs_matrix, const std::string& type = "")
+ : GeneralizedEigenSolverBase(lhs_matrix, rhs_matrix, GeneralizedEigenSolverOptions<MatrixType, true>::options(type))
+ {}
+
+ GeneralizedEigenSolverBase(const MatrixType& lhs_matrix,
+ const MatrixType& rhs_matrix,
+ const Common::Configuration opts)
+ : lhs_matrix_(lhs_matrix)
+ , rhs_matrix_(rhs_matrix)
+ , stored_options_(opts)
+ , options_(&stored_options_)
+ , computed_(false)
+ , disable_checks_(options_->get<bool>("disable_checks", false))
+ {
+ pre_checks();
+ }
+
+ GeneralizedEigenSolverBase(const MatrixType& lhs_matrix, const MatrixType& rhs_matrix, Common::Configuration* opts)
+ : lhs_matrix_(lhs_matrix)
+ , rhs_matrix_(rhs_matrix)
+ , options_(opts)
+ , computed_(false)
+ , disable_checks_(options_->get<bool>("disable_checks", false))
+ {
+ pre_checks();
+ }
+
+ GeneralizedEigenSolverBase(const ThisType& other) = default;
+
+ GeneralizedEigenSolverBase(ThisType&& source) = default;
+
+ virtual ~GeneralizedEigenSolverBase() = default;
+
+ const Common::Configuration& options() const
+ {
+ return *options_;
+ }
+
+ const MatrixType& lhs_matrix() const
+ {
+ return lhs_matrix_;
+ }
+
+ const MatrixType& rhs_matrix() const
+ {
+ return rhs_matrix_;
+ }
+
+protected:
+ /**
+ * \brief Does the actual computation.
+ * \attention The implementor has to fill the appropriate members!
+ * \note The implementor can assume that the given options_ contain a valid 'type' and all default keys.
+ * \nte The implementor does not need to guard against multiple calls of this method.
+ */
+ virtual void compute() const = 0;
+
+public:
+ const std::vector<ComplexType>& eigenvalues() const
+ {
+ compute_and_check();
+ if (eigenvalues_)
+ return *eigenvalues_;
+ else if (options_->get<bool>("compute_eigenvalues"))
+ DUNE_THROW(Common::Exceptions::internal_error, "The eigenvalues_ member is not filled after calling compute()!");
+ else
+ DUNE_THROW(Common::Exceptions::you_are_using_this_wrong,
+ "Do not call eigenvalues() if 'compute_eigenvalues' is false!\n\nThese were the given options:\n\n"
+ << *options_);
+ } // ... eigenvalues(...)
+
+ const std::vector<RealType>& real_eigenvalues() const
+ {
+ compute_and_check();
+ if (eigenvalues_) {
+ if (!real_eigenvalues_)
+ compute_real_eigenvalues();
+ } else if (options_->get<bool>("compute_eigenvalues"))
+ DUNE_THROW(Common::Exceptions::internal_error, "The eigenvalues_ member is not filled after calling compute()!");
+ else
+ DUNE_THROW(
+ Common::Exceptions::you_are_using_this_wrong,
+ "Do not call real_eigenvalues() if 'compute_eigenvalues' is false!\n\nThese were the given options:\n\n"
+ << *options_);
+ assert(real_eigenvalues_ && "These have to exist after compute_real_eigenvalues()!");
+ return *real_eigenvalues_;
+ } // ... real_eigenvalues(...)
+
+ std::vector<RealType> min_eigenvalues(const size_t num_evs = std::numeric_limits<size_t>::max()) const
+ {
+ compute_and_check();
+ if (eigenvalues_) {
+ if (!real_eigenvalues_)
+ compute_real_eigenvalues();
+ } else if (options_->get<bool>("compute_eigenvalues"))
+ DUNE_THROW(Common::Exceptions::internal_error, "The eigenvalues_ member is not filled after calling compute()!");
+ else
+ DUNE_THROW(Common::Exceptions::you_are_using_this_wrong,
+ "Do not call min_eigenvalues() if 'compute_eigenvalues' is false!\n\nThese were the given options:\n\n"
+ << *options_);
+ assert(real_eigenvalues_ && "These have to exist after compute_real_eigenvalues()!");
+ std::vector<RealType> evs = *real_eigenvalues_;
+ std::sort(evs.begin(), evs.end(), [](const RealType& a, const RealType& b) { return a < b; });
+ evs.resize(std::min(evs.size(), num_evs));
+ return evs;
+ } // ... min_eigenvalues(...)
+
+ std::vector<RealType> max_eigenvalues(const size_t num_evs = std::numeric_limits<size_t>::max()) const
+ {
+ compute_and_check();
+ if (eigenvalues_) {
+ if (!real_eigenvalues_)
+ compute_real_eigenvalues();
+ } else if (options_->get<bool>("compute_eigenvalues"))
+ DUNE_THROW(Common::Exceptions::internal_error, "The eigenvalues_ member is not filled after calling compute()!");
+ else
+ DUNE_THROW(Common::Exceptions::you_are_using_this_wrong,
+ "Do not call max_eigenvalues() if 'compute_eigenvalues' is false!\n\nThese were the given options:\n\n"
+ << *options_);
+ assert(real_eigenvalues_ && "These have to exist after compute_real_eigenvalues()!");
+ std::vector<RealType> evs = *real_eigenvalues_;
+ std::sort(evs.begin(), evs.end(), [](const RealType& a, const RealType& b) { return a > b; });
+ evs.resize(std::min(evs.size(), num_evs));
+ return evs;
+ } // ... max_eigenvalues(...)
+
+ const ComplexMatrixType& eigenvectors() const
+ {
+ compute_and_check();
+ if (eigenvectors_)
+ return *eigenvectors_;
+ else if (options_->get<bool>("compute_eigenvectors"))
+ DUNE_THROW(Common::Exceptions::internal_error, "The eigenvectors_ member is not filled after calling compute()!");
+ else
+ DUNE_THROW(Common::Exceptions::you_are_using_this_wrong,
+ "Do not call eigenvectors() if 'compute_eigenvectors' is false!\n\nThese were the given options:\n\n"
+ << *options_);
+ } // ... eigenvectors(...)
+
+ // const ComplexMatrixType& eigenvectors_inverse() const
+ // {
+ // compute_and_check();
+ // if (!eigenvectors_) {
+ // if (options_->get<bool>("compute_eigenvectors"))
+ // DUNE_THROW(Common::Exceptions::internal_error,
+ // "The eigenvectors_ member is not filled after calling compute()!");
+ // else
+ // DUNE_THROW(Common::Exceptions::you_are_using_this_wrong,
+ // "Do not call eigenvectors_inverse() if 'compute_eigenvectors' is false!\n\nThese were the given "
+ // "options:\n\n"
+ // << *options_);
+ // }
+ // invert_eigenvectors();
+ // assert(eigenvectors_inverse_ && "This must not happen after calling invert_eigenvectors()!");
+ // if (!disable_checks_) {
+ // const double check_eigendecomposition = options_->get<double>("assert_eigendecomposition");
+ // if (check_eigendecomposition > 0)
+ // complex_eigendecomposition_helper<>::check(
+ // *this, check_eigendecomposition > 0 ? check_eigendecomposition :
+ // options_->get<double>("real_tolerance"));
+ // }
+ // return *eigenvectors_inverse_;
+ // } // ... eigenvectors(...)
+
+ const RealMatrixType& real_eigenvectors() const
+ {
+ compute_and_check();
+ if (eigenvectors_) {
+ if (!real_eigenvectors_)
+ compute_real_eigenvectors();
+ } else if (options_->get<bool>("compute_eigenvectors"))
+ DUNE_THROW(Common::Exceptions::internal_error, "The eigenvectors_ member is not filled after calling compute()!");
+ else
+ DUNE_THROW(
+ Common::Exceptions::you_are_using_this_wrong,
+ "Do not call real_eigenvectors() if 'compute_eigenvectors' is false!\n\nThese were the given options:\n\n"
+ << *options_);
+ assert(real_eigenvectors_ && "These have to exist after compute_real_eigenvectors()!");
+ return *real_eigenvectors_;
+ } // ... real_eigenvectors(...)
+
+ // const MatrixType& real_eigenvectors_inverse() const
+ // {
+ // compute_and_check();
+ // if (!real_eigenvectors_) {
+ // if (options_->get<double>("assert_real_eigendecomposition") > 0
+ // || options_->get<double>("assert_real_eigenvectors") > 0)
+ // DUNE_THROW(Common::Exceptions::internal_error,
+ // "The real_eigenvectors_ member is not filled after calling compute()!");
+ // else
+ // DUNE_THROW(Common::Exceptions::you_are_using_this_wrong,
+ // "Do not call real_eigenvectors_inverse() after providing these options:\n\n"
+ // << *options_);
+ // }
+ // invert_real_eigenvectors();
+ // assert(real_eigenvectors_inverse_ && "This must not happen after calling invert_real_eigenvectors()!");
+ // if (!disable_checks_) {
+ // const double assert_real_eigendecomposition = options_->get<double>("assert_real_eigendecomposition");
+ // if (assert_real_eigendecomposition > 0.)
+ // assert_eigendecomposition(matrix_,
+ // *real_eigenvalues_,
+ // *real_eigenvectors_,
+ // *real_eigenvectors_inverse_,
+ // assert_real_eigendecomposition);
+ // }
+ // return *real_eigenvectors_inverse_;
+ // } // ... eigenvectors(...)
+
+protected:
+ void compute_and_check() const
+ {
+ if (!computed_) {
+ compute();
+ post_checks();
+ }
+ computed_ = true;
+ }
+
+ void pre_checks()
+ {
+ if (!disable_checks_) {
+ // check options
+ if (!options_->has_key("type"))
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_it_was_not_set_up_correctly,
+ "Missing 'type' in given options:\n\n"
+ << *options_);
+ internal::ensure_generalized_eigen_solver_type(options_->get<std::string>("type"),
+ GeneralizedEigenSolverOptions<MatrixType, true>::types());
+ const Common::Configuration default_opts =
+ GeneralizedEigenSolverOptions<MatrixType, true>::options(options_->get<std::string>("type"));
+ for (const std::string& default_key : default_opts.getValueKeys()) {
+ if (!options_->has_key(default_key))
+ (*options_)[default_key] = default_opts.get<std::string>(default_key);
+ }
+ if (options_->get<double>("real_tolerance") <= 0)
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_it_was_not_set_up_correctly,
+ "It does not make sense to enforce a non-positive tolerance!");
+ if (options_->get<double>("assert_positive_eigenvalues") > 0
+ && options_->get<double>("assert_negative_eigenvalues") > 0)
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_it_was_not_set_up_correctly,
+ "It does not make sense to assert positive and negative eigenvalues!");
+ if (!options_->get<bool>("compute_eigenvalues")
+ && (options_->get<double>("assert_real_eigenvalues") > 0
+ || options_->get<double>("assert_positive_eigenvalues") > 0
+ || options_->get<double>("assert_negative_eigenvalues") > 0
+ /*|| options_->get<double>("assert_eigendecomposition") > 0
+ || options_->get<double>("assert_real_eigendecomposition") > 0*/))
+ (*options_)["compute_eigenvalues"] = "true";
+ if (options_->get<double>("assert_real_eigenvalues") <= 0
+ && (options_->get<double>("assert_positive_eigenvalues") > 0
+ || options_->get<double>("assert_negative_eigenvalues") > 0
+ /*|| options_->get<double>("assert_real_eigendecomposition") > 0*/))
+ (*options_)["assert_real_eigenvalues"] = (*options_)["real_tolerance"];
+ // check matrix
+ check_size(lhs_matrix_, rhs_matrix_);
+ if (options_->get<bool>("check_for_inf_nan") && contains_inf_or_nan(lhs_matrix_)) {
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_data_did_not_fulfill_requirements,
+ "Given LHS matrix contains inf or nan and you requested checking. To disable this check set "
+ "'check_for_inf_nan' to false in the options."
+ << "\n\nThese were the given options:\n\n"
+ << *options_ << "\nThis was the given LHS matrix:\n\n"
+ << lhs_matrix_);
+ }
+ if (options_->get<bool>("check_for_inf_nan") && contains_inf_or_nan(rhs_matrix_)) {
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_data_did_not_fulfill_requirements,
+ "Given RHS matrix contains inf or nan and you requested checking. To disable this check set "
+ "'check_for_inf_nan' to false in the options."
+ << "\n\nThese were the given options:\n\n"
+ << *options_ << "\nThis was the given RHS matrix:\n\n"
+ << rhs_matrix_);
+ }
+ }
+ } // ... pre_checks(...)
+
+ void post_checks() const
+ {
+ if (!disable_checks_) {
+ if (options_->get<bool>("compute_eigenvalues") && !eigenvalues_)
+ DUNE_THROW(Common::Exceptions::internal_error,
+ "The eigenvalues_ member is not filled after calling compute()!");
+ if (options_->get<bool>("compute_eigenvectors") && !eigenvectors_)
+ DUNE_THROW(Common::Exceptions::internal_error,
+ "The eigenvectors_ member is not filled after calling compute()!");
+ if (options_->get<bool>("check_for_inf_nan")) {
+ if (eigenvalues_ && contains_inf_or_nan(*eigenvalues_))
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_result_contained_inf_or_nan,
+ "Computed eigenvalues contain inf or nan and you requested checking. To disable this check set "
+ "'check_for_inf_nan' to false in the options."
+ << "\n\nThese were the given options:\n\n"
+ << *options_ << "\nThese are the computed eigenvalues:\n\n"
+ << *eigenvalues_);
+ if (eigenvectors_ && contains_inf_or_nan(*eigenvectors_))
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_result_contained_inf_or_nan,
+ "Computed eigenvectors contain inf or nan and you requested checking. To disable this check set "
+ "'check_for_inf_nan' to false in the options."
+ << "\n\nThese were the given options:\n\n"
+ << *options_ << "\nThese are the computed eigenvectors:\n\n"
+ << *eigenvectors_);
+ }
+ const double assert_real_eigenvalues = options_->get<double>("assert_real_eigenvalues");
+ const double assert_positive_eigenvalues = options_->get<double>("assert_positive_eigenvalues");
+ const double assert_negative_eigenvalues = options_->get<double>("assert_negative_eigenvalues");
+ // const double check_real_eigendecomposition = options_->get<double>("assert_real_eigendecomposition");
+ if (assert_real_eigenvalues > 0 || assert_positive_eigenvalues > 0 || assert_negative_eigenvalues > 0
+ /*|| check_real_eigendecomposition > 0*/)
+ compute_real_eigenvalues();
+ if (assert_positive_eigenvalues > 0) {
+ assert(real_eigenvalues_ && "This must not happen after compute_real_eigenvalues()!");
+ for (const auto& ev : *real_eigenvalues_) {
+ if (ev < assert_positive_eigenvalues)
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_eigenvalues_are_not_positive_as_requested,
+ "These were the given options:\n\n"
+ << *options_ << "\nThis were the given matrices:\n\nLHS\n"
+ << lhs_matrix_ << "\n\nRHS\n"
+ << rhs_matrix_ << "\nThese are the computed eigenvalues:\n\n"
+ << *real_eigenvalues_);
+ }
+ }
+ if (assert_negative_eigenvalues > 0) {
+ assert(real_eigenvalues_ && "This must not happen after compute_real_eigenvalues()!");
+ for (const auto& ev : *real_eigenvalues_) {
+ if (ev > -1 * assert_negative_eigenvalues)
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_eigenvalues_are_not_negative_as_requested,
+ "These were the given options:\n\n"
+ << *options_ << "\nThis were the given matrices:\n\nLHS\n"
+ << lhs_matrix_ << "\n\nRHS\n"
+ << rhs_matrix_ << "\nThese are the computed eigenvalues:\n\n"
+ << *real_eigenvalues_);
+ }
+ }
+ if (options_->get<double>("assert_real_eigenvectors") > 0 /*|| check_real_eigendecomposition > 0*/)
+ compute_real_eigenvectors();
+ // const double check_eigendecomposition = options_->get<double>("assert_eigendecomposition");
+ // if (check_eigendecomposition > 0)
+ // complex_eigendecomposition_helper<>::check(*this, check_eigendecomposition);
+ // if (check_real_eigendecomposition > 0) {
+ // invert_real_eigenvectors();
+ // assert_eigendecomposition(matrix_,
+ // *real_eigenvalues_,
+ // *real_eigenvectors_,
+ // *real_eigenvectors_inverse_,
+ // check_real_eigendecomposition);
+ // }
+ }
+ } // ... post_checks(...)
+
+ void compute_real_eigenvalues() const
+ {
+ assert(eigenvalues_ && "This should not happen!");
+ if (!real_eigenvalues_) {
+ real_eigenvalues_ = std::make_unique<std::vector<RealType>>(eigenvalues_->size());
+ for (size_t ii = 0; ii < eigenvalues_->size(); ++ii)
+ (*real_eigenvalues_)[ii] = (*eigenvalues_)[ii].real();
+
+ if (!disable_checks_) {
+ const double assert_real_eigenvalues = options_->get<double>("assert_real_eigenvalues");
+ const double tolerance =
+ (assert_real_eigenvalues > 0) ? assert_real_eigenvalues : options_->get<double>("real_tolerance");
+ for (size_t ii = 0; ii < eigenvalues_->size(); ++ii) {
+ if (std::abs((*eigenvalues_)[ii].imag()) > tolerance)
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_eigenvalues_are_not_real_as_requested,
+ "These were the given options:\n\n"
+ << *options_ << "\nThese are the computed eigenvalues:\n\n"
+ << *eigenvalues_);
+ } // ii
+ } // if (!disable_checks_)
+ } // if (!real_eigenvalues_)
+ } // ... compute_real_eigenvalues(...)
+
+ template <bool is_common_matrix = XT::Common::is_matrix<MatrixType>::value, class T = MatrixType>
+ struct real_eigenvectors_helper
+ {};
+
+ template <class T>
+ struct real_eigenvectors_helper<true, T>
+ {
+ static void compute(const ThisType& self, const double& tolerance)
+ {
+ using RM = XT::Common::MatrixAbstraction<RealMatrixType>;
+ using CM = XT::Common::MatrixAbstraction<ComplexMatrixType>;
+ 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));
+ 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) {
+ is_complex = true;
+ ii = rows;
+ break;
+ }
+ RM::set_entry(*self.real_eigenvectors_, ii, jj, complex_value.real());
+ } // jj
+ } // ii
+
+ if (is_complex) {
+ // 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
+ self.compute_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 = (*self.real_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, (*self.real_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 = 0; 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
+ RealType l2_norm = std::sqrt(std::accumulate(
+ v_i.begin(), v_i.end(), 0., [](const RealType& a, const RealType& b) { return a + b * b; }));
+ if (XT::Common::FloatCmp::ne(l2_norm, 0.))
+ v_i *= 1. / 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.))) {
+ if (index >= eigenvalue_multiplicity[kk]) {
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_eigenvectors_are_not_real_as_requested,
+ "Eigenvectors are complex and calculating real eigenvectors failed!"
+ << "These were the given options:\n\n"
+ << *self.options_ << "\n\nThis were the given matrices:\n\nLHS\n"
+ << self.lhs_matrix_ << "\n\nRHS\n"
+ << self.rhs_matrix_ << "\nThese are the computed eigenvectors:\n\n"
+ << std::setprecision(17) << *self.eigenvectors_);
+ }
+ for (size_t rr = 0; rr < rows; ++rr)
+ RM::set_entry(*self.real_eigenvectors_, rr, group[index], input_vectors[ii].get_entry(rr));
+ index++;
+ } // if (input_vectors[ii] != 0)
+ } // ii
+ if (index < eigenvalue_multiplicity[kk]) {
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_eigenvectors_are_not_real_as_requested,
+ "Eigenvectors are complex and calculating real eigenvectors failed!"
+ << "These were the given options:\n\n"
+ << *self.options_ << "\n\nThis were the given matrices:\n\nLHS\n"
+ << self.lhs_matrix_ << "\n\nRHS\n"
+ << self.rhs_matrix_ << "\nThese are the computed eigenvectors:\n\n"
+ << std::setprecision(17) << *self.eigenvectors_);
+ }
+ } // kk
+ } // if(is_complex)
+ } // static void compute(...)
+ }; // real_eigenvectors_helper<true, ...>
+
+ template <class T>
+ struct real_eigenvectors_helper<false, T>
+ {
+ static void compute(ThisType& self, const double& tolerance)
+ {
+ if (RealMatrixType::sparse) {
+ const size_t rows = self.eigenvectors_->rows();
+ const size_t cols = self.eigenvectors_->cols();
+ const auto pattern = self.eigenvectors_->pattern();
+ self.real_eigenvectors_ = std::make_unique<RealMatrixType>(rows, cols, pattern);
+ for (size_t ii = 0; ii < rows; ++ii)
+ for (size_t jj : pattern.inner(ii)) {
+ const auto complex_value = self.eigenvectors_->get_entry(ii, jj);
+ if (std::abs(complex_value.imag()) > tolerance)
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_eigenvectors_are_not_real_as_requested,
+ "These were the given options:\n\n"
+ << *self.options_ << "\nThese are the computed eigenvectors:\n\n"
+ << std::setprecision(17) << *self.eigenvectors_);
+ self.real_eigenvectors_->set_entry(ii, jj, complex_value.real());
+ }
+ } else {
+ const size_t rows = self.eigenvectors_->rows();
+ const size_t cols = self.eigenvectors_->cols();
+ self.real_eigenvectors_ = std::make_unique<RealMatrixType>(rows, cols);
+ for (size_t ii = 0; ii < rows; ++ii)
+ for (size_t jj = 0; jj < cols; ++jj) {
+ const auto complex_value = self.eigenvectors_->get_entry(ii, jj);
+ if (std::abs(complex_value.imag()) > tolerance)
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_eigenvectors_are_not_real_as_requested,
+ "These were the given options:\n\n"
+ << *self.options_ << "\nThese are the computed eigenvectors:\n\n"
+ << std::setprecision(17) << *self.eigenvectors_);
+ self.real_eigenvectors_->set_entry(ii, jj, complex_value.real());
+ }
+ }
+ }
+ }; // real_eigenvectors_helper<false, ...>
+
+ void compute_real_eigenvectors() const
+ {
+ assert(eigenvectors_ && "This should not happen!");
+ if (!real_eigenvectors_) {
+ const double assert_real_eigenvectors = options_->get<double>("assert_real_eigenvectors");
+ const double tolerance =
+ (assert_real_eigenvectors > 0) ? assert_real_eigenvectors : options_->get<double>("real_tolerance");
+ real_eigenvectors_helper<>::compute(*this, tolerance);
+ }
+ }
+
+ // void invert_eigenvectors() const
+ // {
+ // assert(eigenvectors_ && "This must not happen when you call this function!");
+ // try {
+ // if (options_->has_sub("matrix-inverter")) {
+ // eigenvectors_inverse_ =
+ // std::make_unique<ComplexMatrixType>(invert_matrix(*eigenvectors_, options_->sub("matrix-inverter")));
+ // } else
+ // eigenvectors_inverse_ = std::make_unique<ComplexMatrixType>(invert_matrix(*eigenvectors_));
+ // } catch (const Exceptions::matrix_invert_failed& ee) {
+ // DUNE_THROW(Exceptions::generalized_eigen_solver_failed,
+ // "The computed matrix of eigenvectors is not invertible!"
+ // << "\n\nmatrix = " << std::setprecision(17) << matrix_ << "\n\noptions: " << *options_
+ // << "\n\neigenvectors = " << std::setprecision(17) << *eigenvectors_
+ // << "\n\nThis was the original error: " << ee.what());
+ // }
+ // } // ... invert_eigenvectors(...)
+
+ // void invert_real_eigenvectors() const
+ // {
+ // assert(real_eigenvectors_ && "This must not happen when you call this function!");
+ // try {
+ // if (options_->has_sub("matrix-inverter")) {
+ // real_eigenvectors_inverse_ =
+ // std::make_unique<MatrixType>(invert_matrix(*real_eigenvectors_, options_->sub("matrix-inverter")));
+ // } else
+ // real_eigenvectors_inverse_ = std::make_unique<MatrixType>(invert_matrix(*real_eigenvectors_));
+ // } catch (const Exceptions::matrix_invert_failed& ee) {
+ // DUNE_THROW(Exceptions::generalized_eigen_solver_failed,
+ // "The computed matrix of real eigenvectors is not invertible!"
+ // << "\n\nmatrix = " << std::setprecision(17) << matrix_ << "\n\noptions: " << *options_
+ // << "\n\nreal_eigenvectors = " << std::setprecision(17) << *real_eigenvectors_
+ // << "\n\nThis was the original error: " << ee.what());
+ // }
+ // } // ... invert_real_eigenvectors(...)
+
+ // template <class A, class B, class C, class D>
+ // void assert_eigendecomposition(const std::unique_ptr<A>& mat,
+ // const B& eigenvalues,
+ // const C& eigenvectors,
+ // const D& eigenvectors_inverse,
+ // const double& tolerance) const
+ // {
+ // assert_eigendecomposition(*mat, eigenvalues, eigenvectors, eigenvectors_inverse, tolerance);
+ // }
+
+ // template <class A, class B, class C, class D>
+ // void assert_eigendecomposition(const A& mat,
+ // const B& eigenvalues,
+ // const C& eigenvectors,
+ // const D& eigenvectors_inverse,
+ // const double& tolerance) const
+ // {
+ // using M = Common::MatrixAbstraction<C>;
+ // const size_t rows = Common::get_matrix_rows(mat);
+ // const size_t cols = Common::get_matrix_cols(mat);
+ // auto eigenvalue_matrix = M::create(rows, cols, typename M::ScalarType(0.));
+ // for (size_t ii = 0; ii < rows; ++ii)
+ // Common::set_matrix_entry(eigenvalue_matrix, ii, ii, eigenvalues[ii]);
+ // const auto decomposition_error = eigenvectors * eigenvalue_matrix * eigenvectors_inverse - mat;
+ // for (size_t ii = 0; ii < rows; ++ii)
+ // for (size_t jj = 0; jj < cols; ++jj)
+ // if (std::abs(Common::get_matrix_entry(decomposition_error, ii, jj)) > tolerance)
+ // DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_result_is_not_an_eigendecomposition,
+ // "\n\nmatrix = " << std::setprecision(17) << matrix_ << "\n\noptions: " << *options_
+ // << "\n\neigenvalues (lambda)= " << std::setprecision(17) << eigenvalues
+ // << "\n\neigenvectors (T) = " << std::setprecision(17) << eigenvectors
+ // << "\n\n(T * (lambda * T^-1)) - matrix = " << decomposition_error);
+ // } // ... assert_eigendecomposition(...)
+
+ // template <bool upcast_required = !std::is_same<MatrixType, ComplexMatrixType>::value, bool anything = true>
+ // struct complex_eigendecomposition_helper;
+
+ // template <bool anything>
+ // struct complex_eigendecomposition_helper<true, anything>
+ // {
+ // static void check(const ThisType& self, const double& tolerance)
+ // {
+ // self.invert_eigenvectors();
+ // self.assert_eigendecomposition(Dune::XT::LA::convert_to<ComplexMatrixType>(self.matrix_),
+ // *self.eigenvalues_,
+ // *self.eigenvectors_,
+ // *self.eigenvectors_inverse_,
+ // tolerance);
+ // }
+ // };
+
+ // template <bool anything>
+ // struct complex_eigendecomposition_helper<false, anything>
+ // {
+ // static void check(const ThisType& self, const double& tolerance)
+ // {
+ // self.invert_eigenvectors();
+ // self.assert_eigendecomposition(
+ // self.matrix_, *self.eigenvalues_, *self.eigenvectors_, *self.eigenvectors_inverse_, tolerance);
+ // }
+ // };
+
+ template <class M>
+ void check_size(const MatrixInterface<M>& lhs_mat, const MatrixInterface<M>& rhs_mat) const
+ {
+ const auto sz = lhs_mat.rows();
+ if (lhs_mat.cols() != sz)
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_data_did_not_fulfill_requirements,
+ "LHS Matrix has to be square, is " << lhs_mat.rows() << "x" << lhs_mat.cols() << "!");
+ if (rhs_mat.rows() != sz)
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_data_did_not_fulfill_requirements,
+ "Matrices have to be of same size, are " << lhs_mat.rows() << "x" << lhs_mat.cols() << " and "
+ << rhs_mat.rows() << "x" << rhs_mat.cols() << "!");
+ if (rhs_mat.cols() != sz)
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_data_did_not_fulfill_requirements,
+ "RHS Matrix has to be square, is " << rhs_mat.rows() << "x" << rhs_mat.cols() << "!");
+ } // ... check_size(...)
+
+ template <class M>
+ typename std::enable_if<XT::Common::is_matrix<M>::value && !is_matrix<M>::value, void>::type
+ check_size(const M& lhs_mat, const M& rhs_mat) const
+ {
+ using Mat = XT::Common::MatrixAbstraction<M>;
+ const auto sz = Mat::rows(lhs_mat);
+ if (Mat::cols(lhs_mat) != sz)
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_data_did_not_fulfill_requirements,
+ "LHS Matrix has to be square, is " << Mat::rows(lhs_mat) << "x" << Mat::cols(lhs_mat) << "!");
+ if (Mat::rows(rhs_mat) != sz)
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_data_did_not_fulfill_requirements,
+ "Matrices have to be of same size, are " << Mat::rows(lhs_mat) << "x" << Mat::cols(lhs_mat) << " and "
+ << Mat::rows(rhs_mat) << "x" << Mat::cols(rhs_mat) << "!");
+ if (Mat::cols(rhs_mat) != sz)
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_data_did_not_fulfill_requirements,
+ "RHS Matrix has to be square, is " << Mat::rows(rhs_mat) << "x" << Mat::cols(rhs_mat) << "!");
+ } // ... check_size(...)
+
+ template <class T>
+ bool contains_inf_or_nan(const std::vector<T>& vec) const
+ {
+ for (const auto& element : vec)
+ if (XT::Common::isinf(element) || XT::Common::isnan(element))
+ return true;
+ return false;
+ }
+
+ template <class M>
+ bool contains_inf_or_nan(const MatrixInterface<M>& mat) const
+ {
+ return !mat.valid();
+ }
+
+ template <class M>
+ typename std::enable_if<XT::Common::is_matrix<M>::value && !is_matrix<M>::value, bool>::type
+ contains_inf_or_nan(const M& mat) const
+ {
+ using Mat = XT::Common::MatrixAbstraction<M>;
+ for (size_t ii = 0; ii < Mat::rows(mat); ++ii)
+ for (size_t jj = 0; jj < Mat::cols(mat); ++jj) {
+ const auto value = Mat::get_entry(mat, ii, jj);
+ if (XT::Common::isinf(value) || XT::Common::isnan(value))
+ return true;
+ }
+ return false;
+ } // ... contains_inf_or_nan(...)
+
+ const MatrixType& lhs_matrix_;
+ const MatrixType& rhs_matrix_;
+ mutable Common::Configuration stored_options_;
+ mutable Common::Configuration* options_;
+ mutable bool computed_;
+ mutable std::unique_ptr<std::vector<ComplexType>> eigenvalues_;
+ mutable std::unique_ptr<std::vector<RealType>> real_eigenvalues_;
+ mutable std::unique_ptr<ComplexMatrixType> eigenvectors_;
+ mutable std::unique_ptr<RealMatrixType> real_eigenvectors_;
+ // mutable std::unique_ptr<ComplexMatrixType> eigenvectors_inverse_;
+ // mutable std::unique_ptr<RealMatrixType> real_eigenvectors_inverse_;
+ const bool disable_checks_;
+}; // class GeneralizedEigenSolverBase
+
+
+} // namespace internal
+} // namespace LA
+} // namespace XT
+} // namespace Dune
+
+#endif // DUNE_XT_LA_GENERALIZED_EIGEN_SOLVER_INTERNAL_BASE_HH
diff --git a/dune/xt/la/generalized-eigen-solver/internal/lapacke.hh b/dune/xt/la/generalized-eigen-solver/internal/lapacke.hh
new file mode 100644
index 0000000000000000000000000000000000000000..36bf88463a26c89e9ae29993d01775b56c23d12a
--- /dev/null
+++ b/dune/xt/la/generalized-eigen-solver/internal/lapacke.hh
@@ -0,0 +1,170 @@
+// This file is part of the dune-xt-la project:
+// https://github.com/dune-community/dune-xt-la
+// Copyright 2009-2018 dune-xt-la developers and contributors. All rights reserved.
+// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+// or GPL-2.0+ (http://opensource.org/licenses/gpl-license)
+// with "runtime exception" (http://www.dune-project.org/license.html)
+// Authors:
+// Felix Schindler (2019)
+
+#ifndef DUNE_XT_LA_GENERALIZED_EIGEN_SOLVER_INTERNAL_LAPACKE_HH
+#define DUNE_XT_LA_GENERALIZED_EIGEN_SOLVER_INTERNAL_LAPACKE_HH
+
+#include <complex>
+#include <vector>
+#include <string>
+#include <numeric>
+#include <type_traits>
+
+#include <dune/common/typetraits.hh>
+
+#include <dune/xt/common/matrix.hh>
+#include <dune/xt/common/type_traits.hh>
+#include <dune/xt/common/lapacke.hh>
+
+#include <dune/xt/la/container/common/matrix/dense.hh>
+#include <dune/xt/la/exceptions.hh>
+#include <dune/xt/la/type_traits.hh>
+
+#include <dune/xt/la/eigen-solver/internal/lapacke.hh>
+
+namespace Dune {
+namespace XT {
+namespace LA {
+namespace internal {
+
+
+/**
+ * \sa https://software.intel.com/en-us/mkl-developer-reference-c-sygv
+ * \note Most likely, you do not want to use this function directly, but compute_generalized_eigenvalues_using_lapack.
+ */
+template <class RealMatrixType>
+typename std::enable_if<Common::is_matrix<std::decay_t<RealMatrixType>>::value, std::vector<std::complex<double>>>::type
+compute_generalized_eigenvalues_of_real_matrices_using_lapack(RealMatrixType&& lhs_matrix, RealMatrixType&& rhs_matrix)
+{
+ MatrixDataProvider<std::decay_t<RealMatrixType>, is_contiguous_and_mutable<RealMatrixType>::value>
+ lhs_matrix_data_provider(lhs_matrix);
+ MatrixDataProvider<std::decay_t<RealMatrixType>, is_contiguous_and_mutable<RealMatrixType>::value>
+ rhs_matrix_data_provider(rhs_matrix);
+ return compute_generalized_eigenvalues_of_real_matrices_using_lapack_impl(
+ lhs_matrix, lhs_matrix_data_provider, rhs_matrix, rhs_matrix_data_provider);
+}
+
+
+template <class RealMatrixType, bool contiguous_and_mutable>
+typename std::enable_if<Common::is_matrix<RealMatrixType>::value, std::vector<std::complex<double>>>::type
+compute_generalized_eigenvalues_of_real_matrices_using_lapack_impl(
+ const RealMatrixType& lhs_matrix,
+ MatrixDataProvider<RealMatrixType, contiguous_and_mutable>& lhs_matrix_data_provider,
+ const RealMatrixType& rhs_matrix,
+ MatrixDataProvider<RealMatrixType, contiguous_and_mutable>& rhs_matrix_data_provider)
+{
+ if (!Common::Lapacke::available())
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_it_was_not_set_up_correctly,
+ "Do not call any lapack related method if Common::Lapacke::available() is false!");
+ using real_type = typename Dune::XT::Common::MatrixAbstraction<RealMatrixType>::S;
+ static_assert(Dune::XT::Common::is_arithmetic<real_type>::value && !Dune::XT::Common::is_complex<real_type>::value,
+ "Not implemented for complex matrices (yet)!");
+ const size_t size = Dune::XT::Common::get_matrix_rows(lhs_matrix);
+ const auto sz = XT::Common::numeric_cast<int>(size);
+#ifdef DUNE_XT_LA_DISABLE_ALL_CHECKS
+ assert(Dune::XT::Common::get_matrix_cols(lhs_matrix) == size);
+ assert(Dune::XT::Common::get_matrix_rows(rhs_matrix) == size);
+ assert(Dune::XT::Common::get_matrix_cols(rhs_matrix) == size);
+#else
+ if (Dune::XT::Common::get_matrix_cols(lhs_matrix) != size)
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_data_did_not_fulfill_requirements,
+ "Given LHS matrix has to be square, is " << size << "x" << Dune::XT::Common::get_matrix_cols(lhs_matrix)
+ << "!");
+ if (Dune::XT::Common::get_matrix_rows(rhs_matrix) != size)
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_data_did_not_fulfill_requirements,
+ "Given matrices have to be of same size, are " << size << "x"
+ << Dune::XT::Common::get_matrix_cols(lhs_matrix) << "and "
+ << Dune::XT::Common::get_matrix_rows(rhs_matrix) << "x"
+ << Dune::XT::Common::get_matrix_cols(rhs_matrix) << "!");
+ if (Dune::XT::Common::get_matrix_cols(rhs_matrix) != size)
+ DUNE_THROW(Exceptions::generalized_eigen_solver_failed_bc_data_did_not_fulfill_requirements,
+ "Given RHS matrix has to be square, is " << Dune::XT::Common::get_matrix_rows(rhs_matrix) << "x"
+ << Dune::XT::Common::get_matrix_cols(rhs_matrix) << "!");
+#endif // DUNE_XT_LA_DISABLE_ALL_CHECKS
+ thread_local std::vector<double> real_part_of_eigenvalues(size, 0.);
+ int storage_layout = MatrixDataProvider<RealMatrixType, contiguous_and_mutable>::storage_layout
+ == Common::StorageLayout::dense_row_major
+ ? Common::Lapacke::row_major()
+ : Common::Lapacke::col_major();
+ const int info = XT::Common::Lapacke::dsygv(storage_layout,
+ /*problem type=*/1,
+ /*only eigenvalues*/ 'N',
+ /*upper triangles*/ 'U',
+ sz,
+ lhs_matrix_data_provider.data(),
+ sz,
+ rhs_matrix_data_provider.data(),
+ sz,
+ real_part_of_eigenvalues.data());
+ if (info != 0)
+ DUNE_THROW(Dune::XT::LA::Exceptions::generalized_eigen_solver_failed,
+ "The lapack backend reported '"
+ << info << "', see https://software.intel.com/en-us/mkl-developer-reference-c-sygv!");
+ std::vector<std::complex<double>> eigenvalues(size);
+ for (size_t ii = 0; ii < size; ++ii)
+ eigenvalues[ii] = {real_part_of_eigenvalues[ii], 0.};
+ return eigenvalues;
+} // ... compute_generalized_eigenvalues_of_real_matrices_using_lapack_impl(...)
+
+
+/**
+ * \note Most likely, you do not want to use this function directly, but compute_generalized_eigenvalues_using_lapack.
+ */
+template <class MatrixType>
+struct generalized_eigenvalues_lapack_helper
+{
+ static_assert(Common::is_matrix<MatrixType>::value, "");
+
+ template <bool is_complex = Common::is_complex<typename Common::MatrixAbstraction<MatrixType>::S>::value,
+ bool anything = true>
+ struct dtype_switch;
+
+ template <bool anything>
+ struct dtype_switch<true, anything>
+ {
+ template <class MatrixImp>
+ static inline std::vector<std::complex<double>> eigenvalues(MatrixImp&& /*lhs_matrix*/, MatrixImp&& /*rhs_matrix*/)
+ {
+ static_assert(AlwaysFalse<MatrixImp>::value,
+ "Not yet implemented for complex matrices, take a look at "
+ "https://software.intel.com/en-us/mkl-developer-reference-c-sygv "
+ "and add a corresponding free function like "
+ "compute_generalized_eigenvalues_of_real_matrices_using_lapack(...)!");
+ return std::vector<std::complex<double>>();
+ }
+ };
+
+ template <bool anything>
+ struct dtype_switch<false, anything>
+ {
+ template <class MatrixImp>
+ static inline std::vector<std::complex<double>> eigenvalues(MatrixImp&& lhs_matrix, MatrixImp&& rhs_matrix)
+ {
+ return compute_generalized_eigenvalues_of_real_matrices_using_lapack(std::forward<MatrixImp>(lhs_matrix),
+ std::forward<MatrixImp>(rhs_matrix));
+ }
+ };
+}; // class generalized_eigenvalues_lapack_helper
+
+
+template <class MatrixType>
+typename std::enable_if<Common::is_matrix<std::decay_t<MatrixType>>::value, std::vector<std::complex<double>>>::type
+compute_generalized_eigenvalues_using_lapack(MatrixType&& lhs_matrix, MatrixType&& rhs_matrix)
+{
+ return generalized_eigenvalues_lapack_helper<std::decay_t<MatrixType>>::template dtype_switch<>::eigenvalues(
+ std::forward<MatrixType>(lhs_matrix), std::forward<MatrixType>(rhs_matrix));
+}
+
+
+} // namespace internal
+} // namespace LA
+} // namespace XT
+} // namespace Dune
+
+#endif // DUNE_XT_LA_GENERALIZED_EIGEN_SOLVER_INTERNAL_LAPACKE_HH
diff --git a/dune/xt/la/test/generalized-eigensolver.hh b/dune/xt/la/test/generalized-eigensolver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d25eaed99724aecf619debc76b7c46cc93d1f6a7
--- /dev/null
+++ b/dune/xt/la/test/generalized-eigensolver.hh
@@ -0,0 +1,331 @@
+// This file is part of the dune-xt-la project:
+// https://github.com/dune-community/dune-xt-la
+// Copyright 2009-2018 dune-xt-la developers and contributors. All rights reserved.
+// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+// or GPL-2.0+ (http://opensource.org/licenses/gpl-license)
+// with "runtime exception" (http://www.dune-project.org/license.html)
+// Authors:
+// Felix Schindler (2019)
+
+#ifndef DUNE_XT_LA_TEST_GENERALIZED_EIGENSOLVER_HH
+#define DUNE_XT_LA_TEST_GENERALIZED_EIGENSOLVER_HH
+
+#include <dune/xt/common/exceptions.hh>
+#include <dune/xt/common/unused.hh>
+#include <dune/xt/common/logging.hh>
+#include <dune/xt/common/test/gtest/gtest.h>
+
+#include <dune/xt/la/container.hh>
+#include <dune/xt/la/container/conversion.hh>
+#include <dune/xt/la/container/eye-matrix.hh>
+#include <dune/xt/la/generalized-eigen-solver.hh>
+
+using namespace Dune;
+using namespace Dune::XT;
+using namespace Dune::XT::LA;
+
+
+template <class MatrixImp, class F, class C, class R>
+struct GeneralizedEigenSolverTest : public ::testing::Test
+{
+ using MatrixType = MatrixImp;
+ using FieldType = Common::real_t<F>;
+ using ComplexMatrixType = C;
+ using RealMatrixType = R;
+ using RealType = Common::real_t<FieldType>;
+ using EigenValuesType = std::vector<Common::complex_t<FieldType>>;
+ using RealEigenValuesType = std::vector<Common::real_t<FieldType>>;
+ typedef GeneralizedEigenSolver<MatrixType> GeneralizedEigenSolverType;
+ typedef GeneralizedEigenSolverOptions<MatrixType> GeneralizedEigenSolverOpts;
+
+ using M = Common::MatrixAbstraction<MatrixType>;
+
+ GeneralizedEigenSolverTest()
+ : all_matrices_and_expected_eigenvalues_and_vectors_are_computed_(false)
+ , broken_matrix_(M::create(M::has_static_size ? M::static_rows : 1, M::has_static_size ? M::static_cols : 1))
+ , unit_matrix_(
+ eye_matrix<MatrixType>(M::has_static_size ? M::static_rows : 1, M::has_static_size ? M::static_cols : 1))
+ {
+ M::set_entry(broken_matrix_, 0, 0, std::numeric_limits<typename M::S>::infinity());
+ }
+
+ void make_unit_matrix()
+ {
+ ASSERT_TRUE(all_matrices_and_expected_eigenvalues_and_vectors_are_computed_);
+ unit_matrix_ = eye_matrix<MatrixType>(M::has_static_size ? M::static_rows : M::rows(matrix_),
+ M::has_static_size ? M::static_cols : M::cols(matrix_));
+ }
+
+ static void exports_correct_types()
+ {
+ const bool MatrixType_is_correct = std::is_same<typename GeneralizedEigenSolverType::MatrixType, MatrixType>::value;
+ EXPECT_TRUE(MatrixType_is_correct);
+ const bool RealType_is_correct = std::is_same<typename GeneralizedEigenSolverType::RealType, RealType>::value;
+ EXPECT_TRUE(RealType_is_correct);
+ const bool RealMatrixType_is_correct =
+ std::is_same<typename GeneralizedEigenSolverType::RealMatrixType, RealMatrixType>::value;
+ EXPECT_TRUE(RealMatrixType_is_correct);
+ const bool ComplexMatrixType_is_correct =
+ std::is_same<typename GeneralizedEigenSolverType::ComplexMatrixType, ComplexMatrixType>::value;
+ EXPECT_TRUE(ComplexMatrixType_is_correct);
+ } // ... exports_correct_types()
+
+ static void has_types_and_options()
+ {
+ std::vector<std::string> types = GeneralizedEigenSolverOpts::types();
+ EXPECT_GT(types.size(), 0);
+ Common::Configuration opts = GeneralizedEigenSolverOpts::options();
+ EXPECT_TRUE(opts.has_key("type"));
+ EXPECT_EQ(types[0], opts.get<std::string>("type"));
+ for (const auto& tp : types) {
+ Common::Configuration tp_opts = GeneralizedEigenSolverOpts::options(tp);
+ EXPECT_TRUE(tp_opts.has_key("type"));
+ EXPECT_EQ(tp, tp_opts.get<std::string>("type"));
+ }
+ } // ... has_types_and_options(...)
+
+ void throws_on_broken_matrix_construction()
+ {
+ try {
+ GeneralizedEigenSolverType DXTC_UNUSED(default_solver){broken_matrix_, unit_matrix_};
+ FAIL() << "Expected LA::Exceptions::generalized_eigen_solver_failed_bc_data_did_not_fulfill_requirements";
+ } catch (const LA::Exceptions::generalized_eigen_solver_failed_bc_data_did_not_fulfill_requirements& /*ee*/) {
+ } catch (...) {
+ FAIL() << "Expected LA::Exceptions::generalized_eigen_solver_failed_bc_data_did_not_fulfill_requirements";
+ }
+ for (const auto& tp : GeneralizedEigenSolverOpts::types()) {
+ try {
+ GeneralizedEigenSolverType DXTC_UNUSED(default_opts_solver)(broken_matrix_, unit_matrix_, tp);
+ FAIL() << "Expected LA::Exceptions::generalized_eigen_solver_failed_bc_data_did_not_fulfill_requirements";
+ } catch (const LA::Exceptions::generalized_eigen_solver_failed_bc_data_did_not_fulfill_requirements& /*ee*/) {
+ } catch (...) {
+ FAIL() << "Expected LA::Exceptions::generalized_eigen_solver_failed_bc_data_did_not_fulfill_requirements";
+ }
+ try {
+ GeneralizedEigenSolverType DXTC_UNUSED(solver)(
+ broken_matrix_, unit_matrix_, GeneralizedEigenSolverOpts::options(tp));
+ FAIL() << "Expected LA::Exceptions::generalized_eigen_solver_failed_bc_data_did_not_fulfill_requirements";
+ } catch (const LA::Exceptions::generalized_eigen_solver_failed_bc_data_did_not_fulfill_requirements& /*ee*/) {
+ } catch (...) {
+ FAIL() << "Expected LA::Exceptions::generalized_eigen_solver_failed_bc_data_did_not_fulfill_requirements";
+ }
+ }
+ } // ... throws_on_broken_matrix_construction(...)
+
+ void allows_broken_matrix_construction_when_checks_disabled()
+ {
+ for (const auto& tp : GeneralizedEigenSolverOpts::types()) {
+ Common::Configuration opts_with_disabled_check = GeneralizedEigenSolverOpts::options(tp);
+ opts_with_disabled_check["check_for_inf_nan"] = "false";
+ GeneralizedEigenSolverType DXTC_UNUSED(solver)(broken_matrix_, unit_matrix_, opts_with_disabled_check);
+ }
+ }
+
+ void throws_on_inconsistent_given_options()
+ {
+ for (const auto& tp : GeneralizedEigenSolverOpts::types()) {
+ Common::Configuration inconsistent_opts = GeneralizedEigenSolverOpts::options(tp);
+ inconsistent_opts["assert_positive_eigenvalues"] = "1e-15";
+ inconsistent_opts["assert_negative_eigenvalues"] = "1e-15";
+ try {
+ GeneralizedEigenSolverType DXTC_UNUSED(solver)(unit_matrix_, unit_matrix_, inconsistent_opts);
+ FAIL() << "Expected LA::Exceptions::generalized_eigen_solver_failed_bc_it_was_not_set_up_correctly";
+ } catch (const LA::Exceptions::generalized_eigen_solver_failed_bc_it_was_not_set_up_correctly& /*ee*/) {
+ } catch (...) {
+ FAIL() << "Expected LA::Exceptions::generalized_eigen_solver_failed_bc_it_was_not_set_up_correctly";
+ }
+ }
+ } // ... throws_on_inconsistent_given_options(...)
+
+ void is_constructible()
+ {
+ ASSERT_TRUE(all_matrices_and_expected_eigenvalues_and_vectors_are_computed_);
+ this->make_unit_matrix();
+ GeneralizedEigenSolverType DXTC_UNUSED(default_solver){matrix_, unit_matrix_};
+ for (const auto& tp : GeneralizedEigenSolverOpts::types()) {
+ GeneralizedEigenSolverType DXTC_UNUSED(default_opts_solver)(matrix_, unit_matrix_, tp);
+ GeneralizedEigenSolverType DXTC_UNUSED(solver)(matrix_, unit_matrix_, GeneralizedEigenSolverOpts::options(tp));
+ }
+ } // ... is_constructible(...)
+
+ static bool find_ev(const EigenValuesType& expected_evs,
+ const XT::Common::complex_t<FieldType>& actual_ev,
+ const double& tolerance)
+ {
+ for (const auto& expected_ev : expected_evs) {
+ if (Common::FloatCmp::eq(actual_ev, expected_ev, {tolerance, tolerance}))
+ return true;
+ }
+ return false;
+ }
+
+ static bool find_ev(const RealEigenValuesType& expected_evs, const RealType& actual_ev, const double& tolerance)
+ {
+ for (const auto& expected_ev : expected_evs) {
+ if (Common::FloatCmp::eq(actual_ev, expected_ev, tolerance))
+ return true;
+ }
+ return false;
+ }
+
+ void gives_correct_eigenvalues(const Common::Configuration& tolerances = {})
+ {
+ ASSERT_TRUE(all_matrices_and_expected_eigenvalues_and_vectors_are_computed_);
+ this->make_unit_matrix();
+ for (const auto& tp : GeneralizedEigenSolverOpts::types()) {
+ const double tolerance = tolerances.get(tp, 1e-15);
+ GeneralizedEigenSolverType solver(matrix_, unit_matrix_, tp);
+ try {
+ const auto& eigenvalues = solver.eigenvalues();
+ EXPECT_EQ(Common::get_matrix_rows(matrix_), eigenvalues.size());
+ for (const auto& ev : eigenvalues) {
+ EXPECT_TRUE(find_ev(expected_eigenvalues_, ev, tolerance))
+ << "\n\nactual eigenvalue: " << ev << "\n\nexpected eigenvalues: " << expected_eigenvalues_
+ << "\n\ntype: " << tp << "\n\ntolerance: " << tolerance;
+ }
+ } catch (const Dune::MathError&) {
+ if (tolerance > 0) {
+ FAIL() << "Dune::MathError thrown when trying to get eigenvalues!"
+ << "\n\ntype: " << tp << "\n\ntolerance: " << tolerance;
+ }
+ }
+ }
+ } // ... gives_correct_eigenvalues(...)
+
+ bool all_matrices_and_expected_eigenvalues_and_vectors_are_computed_;
+ MatrixType broken_matrix_;
+ MatrixType unit_matrix_;
+ MatrixType matrix_;
+ EigenValuesType expected_eigenvalues_;
+}; // ... struct GeneralizedEigenSolverTest
+
+
+template <class M, class F, class C, class R>
+struct GeneralizedEigenSolverTestForMatricesWithRealEigenvaluesAndVectors
+ : public GeneralizedEigenSolverTest<M, F, C, R>
+{
+ using BaseType = GeneralizedEigenSolverTest<M, F, C, R>;
+ using BaseType::find_ev;
+ using typename BaseType::ComplexMatrixType;
+ using typename BaseType::EigenValuesType;
+ using typename BaseType::FieldType;
+ using typename BaseType::GeneralizedEigenSolverOpts;
+ using typename BaseType::GeneralizedEigenSolverType;
+ using typename BaseType::MatrixType;
+ using typename BaseType::RealEigenValuesType;
+ using typename BaseType::RealMatrixType;
+ using typename BaseType::RealType;
+
+ void gives_correct_real_eigenvalues(const Common::Configuration& tolerances = {})
+ {
+ ASSERT_TRUE(all_matrices_and_expected_eigenvalues_and_vectors_are_computed_);
+ this->make_unit_matrix();
+ for (const auto& tp : GeneralizedEigenSolverOpts::types()) {
+ const double tolerance = tolerances.get(tp, 1e-15);
+ GeneralizedEigenSolverType solver(matrix_, unit_matrix_, tp);
+ const auto& eigenvalues = solver.eigenvalues();
+ EXPECT_EQ(Common::get_matrix_rows(matrix_), eigenvalues.size());
+ for (const auto& complex_ev : eigenvalues) {
+ if (tolerance > 0)
+ EXPECT_TRUE(Common::FloatCmp::eq(0., complex_ev.imag(), tolerance))
+ << "\n type: " << tp << "\n tolerance: " << tolerance;
+ else {
+ // negative tolerance: we expect a failure
+ EXPECT_FALSE(Common::FloatCmp::eq(0., complex_ev.imag()))
+ << "\n\nTHIS IS A GOOD THING! UPDATE THE EXPECTATIONS IN tolerances!\n\n"
+ << "\n type: " << tp;
+ }
+ const auto real_ev = complex_ev.real();
+ if (tolerance > 0)
+ EXPECT_TRUE(find_ev(expected_real_eigenvalues_, real_ev, tolerance))
+ << "\n\nactual eigenvalue: " << real_ev << "\n\nexpected eigenvalues: " << expected_real_eigenvalues_
+ << "\n\ntype: " << tp << "\n\ntolerance: " << tolerance;
+ else {
+ // negative tolerance: we expect a failure
+ EXPECT_FALSE(find_ev(expected_real_eigenvalues_, real_ev, 1e-15))
+ << "\n\nTHIS IS A GOOD THING! UPDATE THE EXPECTATIONS IN tolerances!\n\n"
+ << "\n\nactual eigenvalue: " << real_ev << "\n\nexpected eigenvalues: " << expected_real_eigenvalues_
+ << "\n\ntype: " << tp;
+ }
+ }
+ if (tolerance > 0) {
+ const auto& real_eigenvalues = solver.real_eigenvalues();
+ EXPECT_EQ(Common::get_matrix_rows(matrix_), real_eigenvalues.size());
+ for (const auto& real_ev : real_eigenvalues) {
+ EXPECT_TRUE(find_ev(expected_real_eigenvalues_, real_ev, tolerance))
+ << "\n\nactual eigenvalue: " << real_ev << "\n\nexpected eigenvalues: " << expected_real_eigenvalues_
+ << "\n\ntype: " << tp << "\n\ntolerance: " << tolerance;
+ }
+ } else {
+ // negative tolerance: we expect a failure
+ const auto& real_eigenvalues = solver.real_eigenvalues(); /// \todo: Add try/catch around this, too!
+ EXPECT_EQ(Common::get_matrix_rows(matrix_), real_eigenvalues.size());
+ for (const auto& real_ev : real_eigenvalues) {
+ EXPECT_FALSE(find_ev(expected_real_eigenvalues_, real_ev, 1e-15))
+ << "\n\nTHIS IS A GOOD THING! UPDATE THE EXPECTATIONS IN tolerances!\n\n"
+ << "\n\nactual eigenvalue: " << real_ev << "\n\nexpected eigenvalues: " << expected_real_eigenvalues_
+ << "\n\ntype: " << tp;
+ }
+ }
+ }
+ } // ... gives_correct_real_eigenvalues(...)
+
+ void gives_correct_max_eigenvalue(const Common::Configuration& tolerances = {})
+ {
+ ASSERT_TRUE(all_matrices_and_expected_eigenvalues_and_vectors_are_computed_);
+ this->make_unit_matrix();
+ for (const auto& tp : GeneralizedEigenSolverOpts::types()) {
+ const double tolerance = tolerances.get(tp, 1e-15);
+ GeneralizedEigenSolverType solver(matrix_, unit_matrix_, tp);
+ const auto actual_max_eigenvalues = solver.max_eigenvalues(1); /// \todo: Add try/catch around this, too!
+ ASSERT_GE(1, actual_max_eigenvalues.size());
+ if (tolerance > 0)
+ EXPECT_TRUE(Common::FloatCmp::eq(expected_max_ev_, actual_max_eigenvalues[0], tolerance))
+ << "\n\nactual max eigenvalue: " << actual_max_eigenvalues[0]
+ << "\n\nexpected max eigenvalue: " << expected_max_ev_ << "\n\ntolerance: " << tolerance
+ << "\n\ntype: " << tp;
+ else {
+ // negative tolerance: we expect a failure
+ EXPECT_NE(expected_max_ev_, actual_max_eigenvalues[0])
+ << "\n\nTHIS IS A GOOD THING! UPDATE THE EXPECTATIONS IN tolerances!\n\n"
+ << "\n\nactual max eigenvalue: " << actual_max_eigenvalues[0]
+ << "\n\nexpected max eigenvalue: " << expected_max_ev_ << "\n\ntype: " << tp;
+ }
+ }
+ }
+
+ void gives_correct_min_eigenvalue(const Common::Configuration& tolerances = {})
+ {
+ ASSERT_TRUE(all_matrices_and_expected_eigenvalues_and_vectors_are_computed_);
+ this->make_unit_matrix();
+ for (const auto& tp : GeneralizedEigenSolverOpts::types()) {
+ const double tolerance = tolerances.get(tp, 1e-15);
+ GeneralizedEigenSolverType solver(matrix_, unit_matrix_, tp);
+ const auto actual_min_eigenvalues = solver.min_eigenvalues(1); /// \todo: Add try/catch around this, too!
+ ASSERT_GE(1, actual_min_eigenvalues.size());
+ if (tolerance > 0)
+ EXPECT_TRUE(Common::FloatCmp::eq(expected_min_ev_, actual_min_eigenvalues[0], tolerance))
+ << "\n\nactual min eigenvalue: " << actual_min_eigenvalues[0]
+ << "\n\nexpected min eigenvalue: " << expected_min_ev_ << "\n\ntolerance: " << tolerance
+ << "\n\ntype: " << tp;
+ else {
+ // negative tolerance: we expect a failure
+ EXPECT_NE(expected_min_ev_, actual_min_eigenvalues[0])
+ << "\n\nTHIS IS A GOOD THING! UPDATE THE EXPECTATIONS IN tolerances!\n\n"
+ << "\n\nactual min eigenvalue: " << actual_min_eigenvalues[0]
+ << "\n\nexpected min eigenvalue: " << expected_min_ev_ << "\n\ntype: " << tp;
+ }
+ }
+ }
+
+ using BaseType::all_matrices_and_expected_eigenvalues_and_vectors_are_computed_;
+ using BaseType::expected_eigenvalues_;
+ using BaseType::matrix_;
+ using BaseType::unit_matrix_;
+ RealEigenValuesType expected_real_eigenvalues_;
+ RealType expected_max_ev_;
+ RealType expected_min_ev_;
+}; // struct GeneralizedEigenSolverTestForMatricesWithRealEigenvaluesAndVectors
+
+
+#endif // DUNE_XT_LA_TEST_GENERALIZED_EIGENSOLVER_HH
diff --git a/dune/xt/la/test/generalized-eigensolver_for_real_matrix_with_real_evs.py b/dune/xt/la/test/generalized-eigensolver_for_real_matrix_with_real_evs.py
new file mode 100644
index 0000000000000000000000000000000000000000..eeba057a279a3a08d5ea5cbfa89cc810857e5900
--- /dev/null
+++ b/dune/xt/la/test/generalized-eigensolver_for_real_matrix_with_real_evs.py
@@ -0,0 +1,35 @@
+# ~~~
+# This file is part of the dune-xt-la project:
+# https://github.com/dune-community/dune-xt-la
+# Copyright 2009-2018 dune-xt-la developers and contributors. All rights reserved.
+# License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+# or GPL-2.0+ (http://opensource.org/licenses/gpl-license)
+# with "runtime exception" (http://www.dune-project.org/license.html)
+# Authors:
+# Felix Schindler (2017)
+# René Fritze (2017 - 2019)
+# Tobias Leibner (2017 - 2018)
+# ~~~
+
+from dune.xt.codegen import typeid_to_typedef_name as safe_name, have_eigen
+
+matrix = [
+ 'EigenDenseMatrix<double>', 'FieldMatrix<double, 2, 2>', 'CommonDenseMatrix<double>', 'CommonSparseMatrix<double>'
+]
+field = ['double', 'double', 'double', 'double']
+complex_matrix = [
+ 'EigenDenseMatrix<std::complex<double>>', 'FieldMatrix<std::complex<double>, 2, 2>',
+ 'CommonDenseMatrix<std::complex<double>>', 'CommonSparseMatrix<std::complex<double>>'
+]
+real_matrix = [
+ 'EigenDenseMatrix<double>', 'FieldMatrix<double, 2, 2>', 'CommonDenseMatrix<double>', 'CommonSparseMatrix<double>'
+]
+
+
+def _ok(ft):
+ if 'Eigen' in ft[0]:
+ return have_eigen(cache)
+ return True
+
+
+testtypes = [(safe_name('_'.join(ft)), *ft) for ft in zip(matrix, field, complex_matrix, real_matrix) if _ok(ft)]
diff --git a/dune/xt/la/test/generalized-eigensolver_for_real_matrix_with_real_evs.tpl b/dune/xt/la/test/generalized-eigensolver_for_real_matrix_with_real_evs.tpl
new file mode 100644
index 0000000000000000000000000000000000000000..6b3acd02c12a6bc5fe9f171d3f8a983ef509a797
--- /dev/null
+++ b/dune/xt/la/test/generalized-eigensolver_for_real_matrix_with_real_evs.tpl
@@ -0,0 +1,97 @@
+// This file is part of the dune-xt-la project:
+// https://github.com/dune-community/dune-xt-la
+// Copyright 2009-2017 dune-xt-la developers and contributors. All rights reserved.
+// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+// or GPL-2.0+ (http://opensource.org/licenses/gpl-license)
+// with "runtime exception" (http://www.dune-project.org/license.html)
+// Authors:
+// Felix Schindler (2019)
+
+#include <dune/xt/common/test/main.hxx> // <- has to come first (includes the config.h)!
+
+#include <dune/xt/la/test/generalized-eigensolver.hh>
+
+{% for T_NAME, TESTMATRIXTYPE, TESTFIELDTYPE, TESTCOMPLEXMATRIXTYPE, TESTREALMATRIXTYPE in config.testtypes %}
+struct GeneralizedEigenSolverForMatrixFullOfOnes_{{T_NAME}}
+ : public GeneralizedEigenSolverTestForMatricesWithRealEigenvaluesAndVectors<{{TESTMATRIXTYPE}},
+ {{TESTFIELDTYPE}},
+ {{TESTCOMPLEXMATRIXTYPE}},
+ {{TESTREALMATRIXTYPE}}>
+{
+ using BaseType = GeneralizedEigenSolverTestForMatricesWithRealEigenvaluesAndVectors;
+ using typename BaseType::MatrixType;
+ using typename BaseType::ComplexMatrixType;
+ using typename BaseType::RealMatrixType;
+ using typename BaseType::EigenValuesType;
+ using typename BaseType::RealEigenValuesType;
+
+ GeneralizedEigenSolverForMatrixFullOfOnes_{{T_NAME}}()
+ {
+ matrix_ = XT::Common::from_string<MatrixType>("[1 1; 1 1]");
+ expected_eigenvalues_ = XT::Common::from_string<EigenValuesType>("[2 0]");
+ expected_real_eigenvalues_ = XT::Common::from_string<RealEigenValuesType>("[2 0]");
+ expected_max_ev_ = 2;
+ expected_min_ev_ = 0;
+ all_matrices_and_expected_eigenvalues_and_vectors_are_computed_ = true;
+ }
+
+ using BaseType::all_matrices_and_expected_eigenvalues_and_vectors_are_computed_;
+ using BaseType::matrix_;
+ using BaseType::expected_eigenvalues_;
+ using BaseType::expected_real_eigenvalues_;
+ using BaseType::expected_max_ev_;
+ using BaseType::expected_min_ev_;
+}; // struct GeneralizedEigenSolverForMatrixFullOfOnes_{{T_NAME}}
+
+
+TEST_F(GeneralizedEigenSolverForMatrixFullOfOnes_{{T_NAME}}, exports_correct_types)
+{
+ exports_correct_types();
+}
+
+TEST_F(GeneralizedEigenSolverForMatrixFullOfOnes_{{T_NAME}}, has_types_and_options)
+{
+ has_types_and_options();
+}
+
+TEST_F(GeneralizedEigenSolverForMatrixFullOfOnes_{{T_NAME}}, throws_on_broken_matrix_construction)
+{
+ throws_on_broken_matrix_construction();
+}
+
+TEST_F(GeneralizedEigenSolverForMatrixFullOfOnes_{{T_NAME}}, allows_broken_matrix_construction_when_checks_disabled)
+{
+ allows_broken_matrix_construction_when_checks_disabled();
+}
+
+TEST_F(GeneralizedEigenSolverForMatrixFullOfOnes_{{T_NAME}}, throws_on_inconsistent_given_options)
+{
+ throws_on_inconsistent_given_options();
+}
+
+TEST_F(GeneralizedEigenSolverForMatrixFullOfOnes_{{T_NAME}}, is_constructible)
+{
+ is_constructible();
+}
+
+TEST_F(GeneralizedEigenSolverForMatrixFullOfOnes_{{T_NAME}}, gives_correct_eigenvalues)
+{
+ gives_correct_eigenvalues();
+}
+
+TEST_F(GeneralizedEigenSolverForMatrixFullOfOnes_{{T_NAME}}, gives_correct_real_eigenvalues)
+{
+ gives_correct_real_eigenvalues();
+}
+
+TEST_F(GeneralizedEigenSolverForMatrixFullOfOnes_{{T_NAME}}, gives_correct_max_eigenvalue)
+{
+ gives_correct_max_eigenvalue();
+}
+
+TEST_F(GeneralizedEigenSolverForMatrixFullOfOnes_{{T_NAME}}, gives_correct_min_eigenvalue)
+{
+ gives_correct_min_eigenvalue();
+}
+
+{% endfor %}