DUNE PDELab (git)

arpackpp_geneo.hh
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)
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
30#include <dune/pdelab/backend/interface.hh>
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
44namespace 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
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
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
238 using BlockVector = Dune::PDELab::Backend::Native<BlockVectorWrapper>;
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
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
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)