diff --git a/dune/xt/la/CMakeLists.txt b/dune/xt/la/CMakeLists.txt index 3fc8f62f814fcaa3286ef5bcfd98440af7136d4f..d83863657a50ea28e77eb575cde4a43e8cd58fc0 100644 --- a/dune/xt/la/CMakeLists.txt +++ b/dune/xt/la/CMakeLists.txt @@ -10,8 +10,9 @@ # Tobias Leibner (2016) set(lib_dune_xt_la_sources + algorithms/solve_sym_tridiag_posdef.cc container/pattern.cc - algorithms/solve_sym_tridiag_posdef.cc) + eigen-solver/internal/numpy.cc) if(DUNE_XT_WITH_PYTHON_BINDINGS) list(APPEND lib_dune_xt_la_sources diff --git a/dune/xt/la/eigen-solver/internal/numpy.cc b/dune/xt/la/eigen-solver/internal/numpy.cc new file mode 100644 index 0000000000000000000000000000000000000000..e861ab2a51b441418f467335863f96758a53998f --- /dev/null +++ b/dune/xt/la/eigen-solver/internal/numpy.cc @@ -0,0 +1,38 @@ +// 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 (2017) + +#include "config.h" + +#if HAVE_DUNE_PYBINDXI +#include <dune/pybindxi/interpreter.hh> +#endif + +namespace Dune { +namespace XT { +namespace LA { +namespace internal { + + +bool numpy_eigensolver_available() +{ +#if HAVE_DUNE_PYBINDXI + try { + PybindXI::GlobalInterpreter().import_module("numpy.linalg"); + } catch (...) { + return false; + } +#endif // HAVE_DUNE_PYBINDXI + return true; +} // ... numpy_eigensolver_available(...) + + +} // namespace internal +} // namespace LA +} // namespace XT +} // namespace Dune diff --git a/dune/xt/la/eigen-solver/internal/numpy.hh b/dune/xt/la/eigen-solver/internal/numpy.hh new file mode 100644 index 0000000000000000000000000000000000000000..0f38d4083728fba99005726ae014f6673025d251 --- /dev/null +++ b/dune/xt/la/eigen-solver/internal/numpy.hh @@ -0,0 +1,93 @@ +// 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 (2017) + +#ifndef DUNE_XT_LA_EIGEN_SOLVER_INTERNAL_NUMPY_HH +#define DUNE_XT_LA_EIGEN_SOLVER_INTERNAL_NUMPY_HH + +#include <vector> + +#include <dune/common/typetraits.hh> + +#if HAVE_DUNE_PYBINDXI +#include <dune/pybindxi/interpreter.hh> +#endif + +#include <dune/xt/common/vector.hh> +#include <dune/xt/common/fmatrix.pbh> +#include <dune/xt/common/type_traits.hh> +#include <dune/xt/la/exceptions.hh> + +namespace Dune { +namespace XT { +namespace LA { +namespace internal { + + +bool numpy_eigensolver_available(); + + +#if HAVE_DUNE_PYBINDXI + + +/** + * \todo Extend this for all matrix types M where is_matrix<M>::value or Common::is_matrix<M>::value is true, which + * would require pybind11 type casting for these (should happen in <dune/xt/common/matrix/bindings/hh>). + */ +template <class K, int SIZE> +void compute_eigenvalues_and_right_eigenvectors_of_a_fieldmatrix_using_numpy( + const FieldMatrix<K, SIZE, SIZE>& matrix, + std::vector<Common::complex_t<K>>& eigenvalues, + FieldMatrix<Common::complex_t<K>, SIZE, SIZE>& right_eigenvectors) +{ + if (!numpy_eigensolver_available()) + DUNE_THROW(Exceptions::eigen_solver_failed_bc_it_was_not_set_up_correctly, + "Do not call me if numpy_eigensolver_available() is false!"); + PybindXI::GlobalInterpreter().import_module("warnings").attr("filterwarnings")("error"); + PybindXI::GlobalInterpreter().import_module("numpy").attr("seterr")("raise"); + auto eig = PybindXI::GlobalInterpreter().import_module("numpy.linalg").attr("eig"); + try { + auto result = eig(matrix) + .template cast<std::tuple<FieldVector<Common::complex_t<K>, SIZE>, + FieldMatrix<Common::complex_t<K>, SIZE, SIZE>>>(); + eigenvalues = Common::convert_to<std::vector<Common::complex_t<K>>>(std::get<0>(result)); + right_eigenvectors = std::get<1>(result); + } catch (const pybind11::cast_error&) { + auto result = eig(matrix) + .template cast<std::tuple<FieldVector<Common::real_t<K>, SIZE>, + FieldMatrix<Common::real_t<K>, SIZE, SIZE>>>(); + eigenvalues = Common::convert_to<std::vector<Common::complex_t<K>>>(std::get<0>(result)); + right_eigenvectors = Common::convert_to<FieldMatrix<Common::complex_t<K>, SIZE, SIZE>>(std::get<1>(result)); + } catch (const std::runtime_error& ee) { + DUNE_THROW(Exceptions::eigen_solver_failed, + "Could not convert result!\n\nThis was the original error: " << ee.what()); + } +} // ... compute_eigenvalues_and_right_eigenvectors_of_a_fieldmatrix_using_numpy(...) + + +#else // HAVE_DUNE_PYBINDXI + + +template <class K, int SIZE> +void compute_eigenvalues_and_right_eigenvectors_of_a_fieldmatrix_using_numpy( + const FieldMatrix<K, SIZE, SIZE>& /*matrix*/, + std::vector<Common::complex_t<K>>& /*eigenvalues*/, + FieldMatrix<Common::complex_t<K>, SIZE, SIZE>& /*right_eigenvectors*/) +{ + static_assert(AlwaysFalse<K>::value, "You are missing dune-pybindxi!"); +} + + +#endif // HAVE_DUNE_PYBINDXI + +} // namespace internal +} // namespace LA +} // namespace XT +} // namespace Dune + +#endif // DUNE_XT_LA_EIGEN_SOLVER_INTERNAL_NUMPY_HH