DUNE-FEM (unstable)

istlmatrixadapter.hh
1#ifndef DUNE_FEM_ISTLMATRIXADAPTER_HH
2#define DUNE_FEM_ISTLMATRIXADAPTER_HH
3
4#if HAVE_DUNE_ISTL
5
6#include <dune/common/timer.hh>
7#include <dune/common/version.hh>
8
9#include <dune/istl/operators.hh>
10
11#include <dune/fem/space/lagrange.hh>
12#include <dune/fem/space/discontinuousgalerkin.hh>
13#include <dune/fem/space/padaptivespace/lagrange.hh>
14#include <dune/fem/space/padaptivespace/discontinuousgalerkin.hh>
15#include <dune/fem/space/combinedspace.hh>
16#include <dune/fem/operator/common/operator.hh>
17
18namespace Dune
19{
20 namespace Fem
21 {
22
23 template <class Matrix>
24 class PreconditionerWrapper ;
25
26
27 template <class MatrixImp>
28 class ISTLParallelMatrixAdapterInterface
29 : public AssembledLinearOperator< MatrixImp,
30 typename MatrixImp :: RowBlockVectorType,
31 typename MatrixImp :: ColBlockVectorType>
32 {
33 typedef ISTLParallelMatrixAdapterInterface< MatrixImp > ThisType;
34 public:
35 typedef MatrixImp MatrixType;
36 typedef Fem::PreconditionerWrapper< MatrixType > PreconditionAdapterType;
37
38 typedef typename MatrixType :: RowDiscreteFunctionType RowDiscreteFunctionType;
39 typedef typename MatrixType :: ColDiscreteFunctionType ColumnDiscreteFunctionType;
40
41 typedef typename RowDiscreteFunctionType :: DiscreteFunctionSpaceType RowSpaceType;
42
43 typedef typename ColumnDiscreteFunctionType :: DiscreteFunctionSpaceType ColSpaceType;
44 typedef Fem::ParallelScalarProduct<ColumnDiscreteFunctionType> ParallelScalarProductType;
45
46 typedef typename RowDiscreteFunctionType :: DofStorageType X;
47 typedef typename ColumnDiscreteFunctionType :: DofStorageType Y;
48
50 typedef MatrixType matrix_type;
51 typedef X domain_type;
52 typedef Y range_type;
53 typedef typename X::field_type field_type;
54
55 public:
57 ISTLParallelMatrixAdapterInterface ( const ISTLParallelMatrixAdapterInterface &org )
58 : matrix_( org.matrix_ ),
59 rowSpace_( org.rowSpace_ ),
60 colSpace_( org.colSpace_ ),
61 scp_( colSpace_ ),
62 preconditioner_( org.preconditioner_ ),
63 averageCommTime_( org.averageCommTime_ ),
64 threading_( org.threading_ )
65 {}
66
69 ISTLParallelMatrixAdapterInterface ( MatrixType &matrix,
70 const RowSpaceType &rowSpace,
71 const ColSpaceType &colSpace,
72 const PreconditionAdapterType& precon,
73 const bool threading = true )
74 : matrix_( matrix ),
75 rowSpace_( rowSpace ),
76 colSpace_( colSpace ),
77 scp_( colSpace_ ),
78 preconditioner_( precon ),
79 averageCommTime_( 0 ),
80 threading_( threading )
81 {}
82
84 virtual double averageCommTime() const { return averageCommTime_ ; }
85
87 virtual PreconditionAdapterType &preconditionAdapter() { return preconditioner_; }
88
90 virtual const PreconditionAdapterType &preconditionAdapter() const { return preconditioner_; }
91
93 virtual ParallelScalarProductType &scp() { return scp_; }
94
96 void apply ( const X &x, Y &y ) const override
97 {
98 DUNE_THROW(NotImplemented,"interface method apply not overloaded!");
99 }
100
102 void applyscaleadd ( field_type alpha, const X &x, Y &y) const override
103 {
104 DUNE_THROW(NotImplemented,"interface method applyscaleadd not overloaded!");
105 }
106
108 virtual double residuum(const Y& rhs, X& x) const
109 {
110 DUNE_THROW(NotImplemented,"interface method residuum not overloaded!");
111 return 0.0;
112 }
113
115 const MatrixType& getmat () const override { return matrix_; }
116
118 bool threading () const { return threading_; }
119
120 protected:
121 void communicate( Y &y ) const
122 {
123 if( rowSpace_.grid().comm().size() <= 1 )
124 return;
125
126 Dune::Timer commTime;
127 ColumnDiscreteFunctionType tmp( "LagrangeParallelMatrixAdapter::communicate", colSpace_, y );
128 colSpace_.communicate( tmp );
129 averageCommTime_ += commTime.elapsed();
130 }
131
132 MatrixType &matrix_;
133 const RowSpaceType &rowSpace_;
134 const ColSpaceType &colSpace_;
135 mutable ParallelScalarProductType scp_;
136 PreconditionAdapterType preconditioner_;
137 mutable double averageCommTime_;
138 const bool threading_;
139 };
140
141 template <class MatrixImp>
142 class LagrangeParallelMatrixAdapter
143 : public ISTLParallelMatrixAdapterInterface< MatrixImp >
144 {
145 typedef LagrangeParallelMatrixAdapter< MatrixImp > ThisType;
146 typedef ISTLParallelMatrixAdapterInterface< MatrixImp > BaseType;
147 public:
148 typedef MatrixImp MatrixType;
149 typedef Fem::PreconditionerWrapper<MatrixType> PreconditionAdapterType;
150
151 typedef typename MatrixType :: RowDiscreteFunctionType RowDiscreteFunctionType;
152 typedef typename MatrixType :: ColDiscreteFunctionType ColumnDiscreteFunctionType;
153
154 typedef typename RowDiscreteFunctionType :: DiscreteFunctionSpaceType RowSpaceType;
155
156 typedef typename ColumnDiscreteFunctionType :: DiscreteFunctionSpaceType ColSpaceType;
157 typedef Fem::ParallelScalarProduct<ColumnDiscreteFunctionType> ParallelScalarProductType;
158
159 typedef typename RowDiscreteFunctionType :: DofStorageType X;
160 typedef typename ColumnDiscreteFunctionType :: DofStorageType Y;
161
163 typedef MatrixType matrix_type;
164 typedef X domain_type;
165 typedef Y range_type;
166 typedef typename X::field_type field_type;
167
168 using BaseType :: threading;
169
170 protected:
171 using BaseType :: matrix_;
172 using BaseType :: scp_ ;
173 using BaseType :: colSpace_ ;
174 using BaseType :: rowSpace_ ;
175 using BaseType :: averageCommTime_;
176
177 public:
179 LagrangeParallelMatrixAdapter ( const LagrangeParallelMatrixAdapter &org )
180 : BaseType( org )
181 {}
182
185 LagrangeParallelMatrixAdapter ( MatrixType &matrix,
186 const RowSpaceType &rowSpace,
187 const ColSpaceType &colSpace,
188 const PreconditionAdapterType& precon,
189 const bool threading = true )
190 : BaseType( matrix, rowSpace, colSpace, precon, threading )
191 {}
192
194 void apply ( const X &x, Y &y ) const override
195 {
196 if( threading() )
197 matrix_.mvThreaded(x, y );
198 else
199 matrix_.mv( x, y );
200
201 communicate( y );
202 }
203
205 void applyscaleadd ( field_type alpha, const X &x, Y &y) const override
206 {
207 if( rowSpace_.grid().comm().size() <= 1 )
208 {
209 if( threading() )
210 matrix_.usmvThreaded(alpha, x, y );
211 else
212 matrix_.usmv(alpha,x,y);
213
214 communicate( y );
215 }
216 else
217 {
218 Y tmp( y.size() );
219 apply( x, tmp );
220 y.axpy( alpha, tmp );
221 }
222 }
223
224 double residuum(const Y& rhs, X& x) const override
225 {
226 Y tmp( rhs );
227
228 this->apply(x,tmp);
229 tmp -= rhs;
230
231 // return global sum of residuum
232 return scp_.norm(tmp);
233 }
234
235 SolverCategory::Category category () const override { return SolverCategory::sequential; }
236
237 protected:
238 void communicate( Y &y ) const
239 {
240 if( rowSpace_.grid().comm().size() <= 1 )
241 return;
242
243 Dune::Timer commTime;
244 ColumnDiscreteFunctionType tmp( "LagrangeParallelMatrixAdapter::communicate", colSpace_, y );
245 colSpace_.communicate( tmp );
246 averageCommTime_ += commTime.elapsed();
247 }
248 };
249
250 template <class MatrixImp>
251 class DGParallelMatrixAdapter
252 : public ISTLParallelMatrixAdapterInterface< MatrixImp >
253 {
254 typedef DGParallelMatrixAdapter< MatrixImp > ThisType ;
255 typedef ISTLParallelMatrixAdapterInterface< MatrixImp > BaseType;
256 public:
257 typedef MatrixImp MatrixType;
258 typedef Fem::PreconditionerWrapper<MatrixType> PreconditionAdapterType;
259
260 typedef typename MatrixType :: RowDiscreteFunctionType RowDiscreteFunctionType;
261 typedef typename MatrixType :: ColDiscreteFunctionType ColumnDiscreteFunctionType;
262
263 typedef typename RowDiscreteFunctionType :: DiscreteFunctionSpaceType RowSpaceType;
264
265 typedef typename ColumnDiscreteFunctionType :: DiscreteFunctionSpaceType ColSpaceType;
266 typedef Fem::ParallelScalarProduct<ColumnDiscreteFunctionType> ParallelScalarProductType;
267
268 typedef typename RowDiscreteFunctionType :: DofStorageType X;
269 typedef typename ColumnDiscreteFunctionType :: DofStorageType Y;
270
272 typedef MatrixType matrix_type;
273 typedef X domain_type;
274 typedef Y range_type;
275 typedef typename X::field_type field_type;
276
277 using BaseType :: threading;
278
279 protected:
280 using BaseType :: matrix_;
281 using BaseType :: scp_;
282 using BaseType :: rowSpace_;
283 using BaseType :: colSpace_;
284 using BaseType :: averageCommTime_;
285
286 public:
288 DGParallelMatrixAdapter (const DGParallelMatrixAdapter& org)
289 : BaseType( org )
290 {}
291
293 DGParallelMatrixAdapter (MatrixType& A,
294 const RowSpaceType& rowSpace,
295 const ColSpaceType& colSpace,
296 const PreconditionAdapterType& precon,
297 const bool threading = true)
298 : BaseType( A, rowSpace, colSpace, precon, threading )
299 {}
300
302 void apply (const X& x, Y& y) const override
303 {
304 // exchange data first
305 communicate( x );
306
307 // apply vector to matrix
308 if( threading() )
309 matrix_.mvThreaded( x, y );
310 else
311 matrix_.mv(x,y);
312
313 // delete non-interior
314 scp_.deleteNonInterior( y );
315 }
316
318 void applyscaleadd (field_type alpha, const X& x, Y& y) const override
319 {
320 // exchange data first
321 communicate( x );
322
323 // apply matrix
324 if( threading() )
325 matrix_.usmvThreaded(alpha, x, y );
326 else
327 matrix_.usmv(alpha,x,y);
328
329 // delete non-interior
330 scp_.deleteNonInterior( y );
331 }
332
333 double residuum(const Y& rhs, X& x) const override
334 {
335 // exchange data
336 communicate( x );
337
338 typedef typename ParallelScalarProductType :: AuxiliaryDofsType AuxiliaryDofsType;
339 const AuxiliaryDofsType& auxiliaryDofs = scp_.auxiliaryDofs();
340
341 typedef typename Y :: block_type LittleBlockVectorType;
342 LittleBlockVectorType tmp;
343 double res = 0.0;
344
345 // counter for rows
346 int i = 0;
347 const int auxiliarySize = auxiliaryDofs.size();
348 for(int auxiliary = 0; auxiliary<auxiliarySize; ++auxiliary)
349 {
350 const int nextAuxiliary = auxiliaryDofs[auxiliary];
351 for(; i<nextAuxiliary; ++i)
352 {
353 tmp = 0;
354 // get row
355 typedef typename MatrixType :: row_type row_type;
356
357 const row_type& row = matrix_[i];
358 // multiply with row
359 typedef typename MatrixType :: ConstColIterator ConstColIterator;
360 ConstColIterator endj = row.end();
361 for (ConstColIterator j = row.begin(); j!=endj; ++j)
362 {
363 (*j).umv(x[j.index()], tmp);
364 }
365
366 // substract right hand side
367 tmp -= rhs[i];
368
369 // add scalar product
370 res += tmp.two_norm2();
371 }
372 ++i;
373 }
374
375 res = rowSpace_.grid().comm().sum( res );
376 // return global sum of residuum
377 return std::sqrt( res );
378 }
379
380 SolverCategory::Category category () const override { return SolverCategory::sequential; }
381
382 protected:
383 void communicate(const X& x) const
384 {
385 if( rowSpace_.grid().comm().size() <= 1 ) return ;
386
387 Dune::Timer commTime;
388
389 // create temporary discrete function object
390 RowDiscreteFunctionType tmp ("DGParallelMatrixAdapter::communicate",
391 rowSpace_, x );
392
393 // exchange data by copying
394 rowSpace_.communicate( tmp );
395
396 // accumulate communication time
397 averageCommTime_ += commTime.elapsed();
398 }
399 };
400
401 } // end namespace Fem
402} // end namespace Dune
403
404#endif // #if HAVE_DUNE_ISTL
405
406#endif // #ifndef DUNE_FEM_ISTLMATRIXADAPTER_HH
407
Dune::Timer
A simple stop watch.
Definition: timer.hh:43
Dune::Timer::elapsed
double elapsed() const noexcept
Get elapsed user-time from last reset until now/last stop in seconds.
Definition: timer.hh:77
version.hh
Various macros to work with Dune module version numbers.
DUNE_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
Dune
Dune namespace.
Definition: alignedallocator.hh:13
operators.hh
Define general, extensible interface for operators. The available implementation wraps a matrix.
Dune::SolverCategory::Category
Category
Definition: solvercategory.hh:23
Dune::SolverCategory::sequential
@ sequential
Category for sequential solvers.
Definition: solvercategory.hh:25
timer.hh
A simple timing class.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 27, 22:29, 2024)