6#ifndef DUNE_PDELAB_BACKEND_ISTL_GENEO_ARPACK_GENEO_HH
7#define DUNE_PDELAB_BACKEND_ISTL_GENEO_ARPACK_GENEO_HH
20#include <dune/common/fvector.hh>
21#include <dune/common/exceptions.hh>
23#if DUNE_VERSION_GTE(ISTL,2,8)
24#include <dune/istl/blocklevel.hh>
26#include <dune/istl/bvector.hh>
27#include <dune/istl/istlexception.hh>
28#include <dune/istl/io.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);
101 Dune::InverseOperatorResult result;
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;
136 typedef Dune::FieldVector<Real,dbvBlockSize> DomainBlockVectorBlock;
137 typedef Dune::BlockVector<DomainBlockVectorBlock> DomainBlockVector;
141 constexpr static int rbvBlockSize = mBlock;
142 typedef Dune::FieldVector<Real,rbvBlockSize> RangeBlockVectorBlock;
143 typedef Dune::BlockVector<RangeBlockVectorBlock> RangeBlockVector;
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_;
203 Dune::UMFPack<BCRSMatrix> A_solver;
210 mutable DomainBlockVector domainBlockVector;
211 mutable RangeBlockVector rangeBlockVector;
233 template <
typename BCRSMatrix,
typename BlockVectorWrapper>
234 class ArPackPlusPlus_Algorithms
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();
305 DUNE_THROW(Dune::ISTLError,
"Matrix is not square ("
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);
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)
356 DUNE_THROW(Dune::ISTLError,
"No algorithm applied, yet.");
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_;
VTKWriter & w
Definition: function.hh:842
std::size_t index
Definition: interpolate.hh:97
typename native_type< T >::type Native
Alias of the native container type associated with T or T itself if it is not a backend wrapper.
Definition: backend/interface.hh:176