1#ifndef DUNE_PDELAB_BACKEND_ISTL_GENEO_GENEOBASIS_HH
2#define DUNE_PDELAB_BACKEND_ISTL_GENEO_GENEOBASIS_HH
7#include <dune/pdelab/backend/istl/geneo/subdomainbasis.hh>
8#include <dune/pdelab/backend/istl/geneo/arpackpp_geneo.hh>
21 template<
class GFS,
class M,
class X,
int dim>
22 class GenEOBasis :
public SubdomainBasis<X>
24 typedef Dune::PDELab::Backend::Native<M> ISTLM;
25 typedef Dune::PDELab::Backend::Native<X> ISTLX;
42 GenEOBasis(
const GFS& gfs,
const M& AF_exterior,
const M& AF_ovlp,
const double eigenvalue_threshold, X& part_unity,
43 int& nev,
int nev_arpack = -1,
double shift = 0.001,
bool add_part_unity =
false,
int verbose = 0) {
44 using Dune::PDELab::Backend::native;
51 const int block_size = X::block_type::dimension;
55 for (
auto row_iter = native(ovlp_mat).begin(); row_iter != native(ovlp_mat).end(); row_iter++) {
56 for (
auto col_iter = row_iter->begin(); col_iter != row_iter->end(); col_iter++) {
57 for (
int block_row = 0; block_row < block_size; block_row++) {
58 for (
int block_col = 0; block_col < block_size; block_col++) {
59 (*col_iter)[block_row][block_col] *= native(part_unity)[row_iter.index()][block_row] * native(part_unity)[col_iter.index()][block_col];
66 ArpackGeneo::ArPackPlusPlus_Algorithms<ISTLM, X> arpack(native(AF_exterior));
69 std::vector<double> eigenvalues(nev_arpack,0.0);
70 std::vector<X> eigenvectors(nev_arpack,X(gfs,0.0));
72 arpack.computeGenNonSymMinMagnitude(native(ovlp_mat), eps, eigenvectors, eigenvalues, shift);
76 if (eigenvalue_threshold >= 0) {
77 for (
int i = 0; i < nev; i++) {
78 if (eigenvalues[i] > eigenvalue_threshold) {
84 std::cout <<
"Process " << gfs.gridView().comm().rank() <<
" picked " << cnt <<
" eigenvectors" << std::endl;
92 this->local_basis.resize(cnt);
93 for (
int base_id = 0; base_id < cnt; base_id++) {
94 this->local_basis[base_id] = std::make_shared<X>(part_unity);
97 this->local_basis[base_id]->begin(),this->local_basis[base_id]->end(),
98 eigenvectors[base_id].begin(),
99 this->local_basis[base_id]->begin(),
105 for (
auto& v : this->local_basis) {
106 *v *= 1.0 / v->two_norm2();
111 if (add_part_unity && eigenvalues[0] > 1E-10) {
112 this->local_basis.insert (this->local_basis.begin(), std::make_shared<X>(part_unity));
113 this->local_basis.pop_back();
Base class for Dune-Exceptions.
Definition: exceptions.hh:96
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
Dune namespace.
Definition: alignedallocator.hh:13