dune-pdelab  2.7-git
arpackpp_geneo.hh
Go to the documentation of this file.
1 /*
2  * Modified version of the ISTL Arpack++ wrapper for supporting
3  * generalized eigenproblems as required by the GenEO coarse space.
4  */
5 
6 #ifndef DUNE_PDELAB_BACKEND_ISTL_GENEO_ARPACK_GENEO_HH
7 #define DUNE_PDELAB_BACKEND_ISTL_GENEO_ARPACK_GENEO_HH
8 
9 #if HAVE_ARPACKPP
10 
11 #include <cmath> // provides std::abs, std::pow, std::sqrt
12 
13 #include <iostream> // provides std::cout, std::endl
14 #include <string> // provides std::string
15 
16 #include <algorithm>
17 #include <numeric>
18 #include <vector>
19 
20 #include <dune/common/fvector.hh> // provides Dune::FieldVector
21 #include <dune/common/exceptions.hh> // provides DUNE_THROW(...)
22 
23 #if DUNE_VERSION_GTE(ISTL,2,8)
24 #include <dune/istl/blocklevel.hh>
25 #endif
26 #include <dune/istl/bvector.hh> // provides Dune::BlockVector
27 #include <dune/istl/istlexception.hh> // provides Dune::ISTLError
28 #include <dune/istl/io.hh> // provides Dune::printvector(...)
29 
31 
32 #ifdef Status
33 #undef Status // prevent preprocessor from damaging the ARPACK++
34  // code when "X11/Xlib.h" is included (the latter
35  // defines Status as "#define Status int" and
36  // ARPACK++ provides a class with a method called
37  // Status)
38 #endif
39 #include "arssym.h" // provides ARSymStdEig
40 #include "argsym.h" // provides ARSymGenEig
41 #include "argnsym.h" // provides ARSymGenEig
42 
43 
44 namespace ArpackGeneo
45 {
46 
59  template <class BCRSMatrix>
60  class ArPackPlusPlus_BCRSMatrixWrapperGen
61  {
62  public:
64  typedef typename BCRSMatrix::field_type Real;
65 
69  public:
71  ArPackPlusPlus_BCRSMatrixWrapperGen (const BCRSMatrix& A, const BCRSMatrix& B)
72  : A_(A), B_(B),
73  A_solver(A_),
74  a_(A_.M() * mBlock), n_(A_.N() * nBlock)
75  {
76  // assert that BCRSMatrix type has blocklevel 2
77 #if DUNE_VERSION_LT(DUNE_ISTL,2,8)
78  static_assert
79  (BCRSMatrix::blocklevel == 2,
80  "Only BCRSMatrices with blocklevel 2 are supported.");
81 #else
82  static_assert
83  (Dune::blockLevel<BCRSMatrix>() == 2,
84  "Only BCRSMatrices with blocklevel 2 are supported.");
85 #endif
86  // allocate memory for auxiliary block vector objects
87  // which are compatible to matrix rows / columns
88  domainBlockVector.resize(A_.N());
89  rangeBlockVector.resize(A_.M());
90  }
91 
93  inline void multMv (Real* v, Real* w)
94  {
95  // get vector v as an object of appropriate type
96  arrayToDomainBlockVector(v,domainBlockVector);
97 
98  // perform matrix-vector product
99  B_.mv(domainBlockVector,rangeBlockVector);
100 
101  Dune::InverseOperatorResult result;
102  auto rangeBlockVector2 = rangeBlockVector;
103  A_solver.apply(rangeBlockVector2, rangeBlockVector, result);
104 
105  // get vector w from object of appropriate type
106  rangeBlockVectorToArray(rangeBlockVector2,w);
107  };
108 
109  inline void multMvB (Real* v, Real* w)
110  {
111  // get vector v as an object of appropriate type
112  arrayToDomainBlockVector(v,domainBlockVector);
113 
114  // perform matrix-vector product
115  B_.mv(domainBlockVector,rangeBlockVector);
116 
117  // get vector w from object of appropriate type
118  rangeBlockVectorToArray(rangeBlockVector,w);
119  };
120 
121 
123  inline int nrows () const { return a_; }
124 
126  inline int ncols () const { return n_; }
127 
128  protected:
129  // Number of rows and columns in each block of the matrix
130  constexpr static int mBlock = BCRSMatrix::block_type::rows;
131  constexpr static int nBlock = BCRSMatrix::block_type::cols;
132 
133  // Type of vectors in the domain of the linear map associated with
134  // the matrix, i.e. block vectors compatible to matrix rows
135  constexpr static int dbvBlockSize = nBlock;
136  typedef Dune::FieldVector<Real,dbvBlockSize> DomainBlockVectorBlock;
137  typedef Dune::BlockVector<DomainBlockVectorBlock> DomainBlockVector;
138 
139  // Type of vectors in the range of the linear map associated with
140  // the matrix, i.e. block vectors compatible to matrix columns
141  constexpr static int rbvBlockSize = mBlock;
142  typedef Dune::FieldVector<Real,rbvBlockSize> RangeBlockVectorBlock;
143  typedef Dune::BlockVector<RangeBlockVectorBlock> RangeBlockVector;
144 
145  // Types for vector index access
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;
150 
151  // Get vector v from a block vector object which is compatible to
152  // matrix rows
153  static inline void
154  domainBlockVectorToArray (const DomainBlockVector& dbv, Real* v)
155  {
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];
159  }
160 
161  // Get vector v from a block vector object which is compatible to
162  // matrix columns
163  static inline void
164  rangeBlockVectorToArray (const RangeBlockVector& rbv, Real* v)
165  {
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];
169  }
170 
171  public:
174  static inline void arrayToDomainBlockVector (const Real* v,
175  DomainBlockVector& dbv)
176  {
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];
180  }
181 
182  template<typename Vector>
183  static inline void arrayToVector(const Real* data, Vector& v)
184  {
185  std::copy(data,data + v.flatsize(),v.begin());
186  }
187 
190  static inline void arrayToRangeBlockVector (const Real* v,
191  RangeBlockVector& rbv)
192  {
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];
196  }
197 
198  protected:
199  // The DUNE-ISTL BCRSMatrix
200  const BCRSMatrix& A_;
201  const BCRSMatrix& B_;
202 
203  Dune::UMFPack<BCRSMatrix> A_solver;
204 
205  // Number of rows and columns in the matrix
206  const int a_, n_;
207 
208  // Auxiliary block vector objects which are
209  // compatible to matrix rows / columns
210  mutable DomainBlockVector domainBlockVector;
211  mutable RangeBlockVector rangeBlockVector;
212  };
213 
214 
233  template <typename BCRSMatrix, typename BlockVectorWrapper>
234  class ArPackPlusPlus_Algorithms
235  {
236  public:
237 
239  typedef typename BlockVector::field_type Real;
240 
241  public:
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),
265  nIterations_(0),
266  title_(" ArPackPlusPlus_Algorithms: "),
267  blank_(title_.length(),' ')
268  {}
269 
270 
271  inline void computeGenNonSymMinMagnitude (const BCRSMatrix& b_, const Real& epsilon,
272  std::vector<BlockVectorWrapper>& x, std::vector<Real>& lambda, Real sigma) const
273  {
274  // print verbosity information
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;
279 
280 
281 
282  // allocate memory for variables, set parameters
283  const int nev = x.size(); // Number of eigenvalues to compute
284  const int ncv = 0; // Number of Arnoldi vectors generated at each iteration (0 == auto)
285  const Real tol = epsilon; // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
286  const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed (0 == 100*nev)
287  auto ev = std::vector<Real>(nev); // Computed generalized eigenvalues
288  auto ev_imag = std::vector<Real>(nev); // Computed generalized eigenvalues
289  const bool ivec = true; // Flag deciding if eigenvectors shall be determined
290 
291  BCRSMatrix ashiftb(a_);
292  ashiftb.axpy(-sigma,b_);
293 
294  // use type ArPackPlusPlus_BCRSMatrixWrapperGen to store matrix information
295  // and to perform the product (A-sigma B)^-1 v (LU decomposition is not used)
296  typedef ArPackPlusPlus_BCRSMatrixWrapperGen<BCRSMatrix> WrappedMatrix;
297  WrappedMatrix A(ashiftb,b_);
298 
299  // get number of rows and columns in A
300  const int nrows = A.nrows();
301  const int ncols = A.ncols();
302 
303  // assert that A is square
304  if (nrows != ncols)
305  DUNE_THROW(Dune::ISTLError,"Matrix is not square ("
306  << nrows << "x" << ncols << ").");
307 
308  // define what we need: eigenvalues with smallest magnitude
309  char which[] = "LM";
310  //ARNonSymGenEig<Real,WrappedMatrix,WrappedMatrix>
311  // dprob(nrows, nev, &A, &WrappedMatrix::multMv, &A, &WrappedMatrix::multMvB, sigma, which, ncv, tol, maxit);
312 
313  //Non generalised problem
314  ARNonSymStdEig<Real,WrappedMatrix>
315  dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
316 
317 
318  // set ARPACK verbosity mode if requested
319  if (verbosity_level_ > 3) dprob.Trace();
320 
321  // find eigenvalues and eigenvectors of A
322  // The broken interface of ARPACK++ actually wants a reference to a pointer...
323  auto ev_data = ev.data();
324  auto ev_imag_data = ev_imag.data();
325  dprob.Eigenvalues(ev_data,ev_imag_data,ivec);
326 
327  // Get sorting permutation for un-shifted eigenvalues
328  std::vector<int> index(nev, 0);
329  std::iota(index.begin(),index.end(),0);
330 
331  std::sort(index.begin(), index.end(),
332  [&](const int& a, const int& b) {
333  return (sigma+1./ev[a] < sigma+1./ev[b]);
334  }
335  );
336 
337  // Unshift eigenpairs
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]);
342  }
343 
344  // obtain number of Arnoldi update iterations actually taken
345  nIterations_ = dprob.GetIter();
346 
347  }
348 
353  inline unsigned int getIterationCount () const
354  {
355  if (nIterations_ == 0)
356  DUNE_THROW(Dune::ISTLError,"No algorithm applied, yet.");
357 
358  return nIterations_;
359  }
360 
361  protected:
362  // parameters related to iterative eigenvalue algorithms
363  const BCRSMatrix& a_;
364  const unsigned int nIterationsMax_;
365 
366  // verbosity setting
367  const unsigned int verbosity_level_;
368 
369  // memory for storing temporary variables (mutable as they shall
370  // just be effectless auxiliary variables of the const apply*(...)
371  // methods)
372  mutable unsigned int nIterations_;
373 
374  // constants for printing verbosity information
375  const std::string title_;
376  const std::string blank_;
377  };
378 
379 } // namespace Dune
380 
381 
382 #endif // HAVE_ARPACKPP
383 
384 #endif // DUNE_PDELAB_BACKEND_ISTL_GENEO_ARPACK_GENEO_HH
VTKWriter & w
Definition: function.hh:842
std::size_t index
Definition: interpolate.hh:97
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:106
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