6#ifndef DUNE_PDELAB_BACKEND_ISTL_GENEO_ARPACK_GENEO_HH
7#define DUNE_PDELAB_BACKEND_ISTL_GENEO_ARPACK_GENEO_HH
23#if DUNE_VERSION_GTE(ISTL,2,8)
27#include <dune/istl/istlexception.hh>
30#include <dune/pdelab/backend/interface.hh>
59 template <
class BCRSMatrix>
60 class ArPackPlusPlus_BCRSMatrixWrapperGen
64 typedef typename BCRSMatrix::field_type Real;
71 ArPackPlusPlus_BCRSMatrixWrapperGen (
const BCRSMatrix& A,
const BCRSMatrix& B)
74 a_(A_.M() * mBlock), n_(A_.N() * nBlock)
77#if DUNE_VERSION_LT(DUNE_ISTL,2,8)
79 (BCRSMatrix::blocklevel == 2,
80 "Only BCRSMatrices with blocklevel 2 are supported.");
83 (Dune::blockLevel<BCRSMatrix>() == 2,
84 "Only BCRSMatrices with blocklevel 2 are supported.");
88 domainBlockVector.resize(A_.N());
89 rangeBlockVector.resize(A_.M());
93 inline void multMv (Real* v, Real* w)
96 arrayToDomainBlockVector(v,domainBlockVector);
99 B_.mv(domainBlockVector,rangeBlockVector);
102 auto rangeBlockVector2 = rangeBlockVector;
103 A_solver.apply(rangeBlockVector2, rangeBlockVector, result);
106 rangeBlockVectorToArray(rangeBlockVector2,w);
109 inline void multMvB (Real* v, Real* w)
112 arrayToDomainBlockVector(v,domainBlockVector);
115 B_.mv(domainBlockVector,rangeBlockVector);
118 rangeBlockVectorToArray(rangeBlockVector,w);
123 inline int nrows ()
const {
return a_; }
126 inline int ncols ()
const {
return n_; }
130 constexpr static int mBlock = BCRSMatrix::block_type::rows;
131 constexpr static int nBlock = BCRSMatrix::block_type::cols;
135 constexpr static int dbvBlockSize = nBlock;
141 constexpr static int rbvBlockSize = mBlock;
146 typedef typename DomainBlockVector::size_type dbv_size_type;
147 typedef typename RangeBlockVector::size_type rbv_size_type;
148 typedef typename DomainBlockVectorBlock::size_type dbvb_size_type;
149 typedef typename RangeBlockVectorBlock::size_type rbvb_size_type;
154 domainBlockVectorToArray (
const DomainBlockVector& dbv, Real* v)
156 for (dbv_size_type block = 0; block < dbv.N(); ++block)
157 for (dbvb_size_type iBlock = 0; iBlock < dbvBlockSize; ++iBlock)
158 v[block*dbvBlockSize + iBlock] = dbv[block][iBlock];
164 rangeBlockVectorToArray (
const RangeBlockVector& rbv, Real* v)
166 for (rbv_size_type block = 0; block < rbv.N(); ++block)
167 for (rbvb_size_type iBlock = 0; iBlock < rbvBlockSize; ++iBlock)
168 v[block*rbvBlockSize + iBlock] = rbv[block][iBlock];
174 static inline void arrayToDomainBlockVector (
const Real* v,
175 DomainBlockVector& dbv)
177 for (dbv_size_type block = 0; block < dbv.N(); ++block)
178 for (dbvb_size_type iBlock = 0; iBlock < dbvBlockSize; ++iBlock)
179 dbv[block][iBlock] = v[block*dbvBlockSize + iBlock];
182 template<
typename Vector>
183 static inline void arrayToVector(
const Real* data, Vector& v)
185 std::copy(data,data + v.flatsize(),v.begin());
190 static inline void arrayToRangeBlockVector (
const Real* v,
191 RangeBlockVector& rbv)
193 for (rbv_size_type block = 0; block < rbv.N(); ++block)
194 for (rbvb_size_type iBlock = 0; iBlock < rbvBlockSize; ++iBlock)
195 rbv[block][iBlock] = v[block*rbvBlockSize + iBlock];
200 const BCRSMatrix& A_;
201 const BCRSMatrix& B_;
210 mutable DomainBlockVector domainBlockVector;
211 mutable RangeBlockVector rangeBlockVector;
233 template <
typename BCRSMatrix,
typename BlockVectorWrapper>
234 class ArPackPlusPlus_Algorithms
238 using BlockVector = Dune::PDELab::Backend::Native<BlockVectorWrapper>;
239 typedef typename BlockVector::field_type Real;
260 ArPackPlusPlus_Algorithms (
const BCRSMatrix& m,
261 const unsigned int nIterationsMax = 100000,
262 const unsigned int verbosity_level = 0)
263 : a_(m), nIterationsMax_(nIterationsMax),
264 verbosity_level_(verbosity_level),
266 title_(
" ArPackPlusPlus_Algorithms: "),
267 blank_(title_.length(),
' ')
271 inline void computeGenNonSymMinMagnitude (
const BCRSMatrix& b_,
const Real& epsilon,
272 std::vector<BlockVectorWrapper>& x, std::vector<Real>& lambda, Real sigma)
const
275 if (verbosity_level_ > 0)
276 std::cout << title_ <<
"Computing an approximation of the "
277 <<
"least dominant eigenvalue of a matrix which "
278 <<
"is assumed to be symmetric." << std::endl;
283 const int nev = x.size();
285 const Real tol = epsilon;
286 const int maxit = nIterationsMax_*nev;
287 auto ev = std::vector<Real>(nev);
288 auto ev_imag = std::vector<Real>(nev);
289 const bool ivec =
true;
291 BCRSMatrix ashiftb(a_);
292 ashiftb.axpy(-sigma,b_);
296 typedef ArPackPlusPlus_BCRSMatrixWrapperGen<BCRSMatrix> WrappedMatrix;
297 WrappedMatrix A(ashiftb,b_);
300 const int nrows = A.nrows();
301 const int ncols = A.ncols();
306 << nrows <<
"x" << ncols <<
").");
314 ARNonSymStdEig<Real,WrappedMatrix>
315 dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
319 if (verbosity_level_ > 3) dprob.Trace();
323 auto ev_data = ev.data();
324 auto ev_imag_data = ev_imag.data();
325 dprob.Eigenvalues(ev_data,ev_imag_data,ivec);
328 std::vector<int> index(nev, 0);
329 std::iota(index.begin(),index.end(),0);
331 std::sort(index.begin(), index.end(),
332 [&](
const int& a,
const int& b) {
333 return (sigma+1./ev[a] < sigma+1./ev[b]);
338 for (
int i = 0; i < nev; i++) {
339 lambda[i] = sigma+1./ev[index[i]];
340 Real* x_raw = dprob.RawEigenvector(index[i]);
341 WrappedMatrix::arrayToVector(x_raw,x[i]);
345 nIterations_ = dprob.GetIter();
353 inline unsigned int getIterationCount ()
const
355 if (nIterations_ == 0)
363 const BCRSMatrix& a_;
364 const unsigned int nIterationsMax_;
367 const unsigned int verbosity_level_;
372 mutable unsigned int nIterations_;
375 const std::string title_;
376 const std::string blank_;
Helper functions for determining the vector/matrix block level.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
A vector of blocks with memory management.
Definition: bvector.hh:392
vector space out of a tensor product of fields.
Definition: fvector.hh:91
derive error class from the base class in common
Definition: istlexception.hh:19
The UMFPack direct sparse solver.
Definition: umfpack.hh:258
A few common exception classes.
Some generic functions for pretty printing vectors and matrices.
Implements a vector constructed from a given type representing a field and a compile-time given size.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
Statistics about the application of an inverse operator.
Definition: solver.hh:50