dune-fem 2.8.0
Loading...
Searching...
No Matches
istlmatrixadapter.hh
Go to the documentation of this file.
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
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#if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
57 enum { category=SolverCategory::sequential };
58#endif // #if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
59
60 public:
62 ISTLParallelMatrixAdapterInterface ( const ISTLParallelMatrixAdapterInterface &org )
63 : matrix_( org.matrix_ ),
64 rowSpace_( org.rowSpace_ ),
65 colSpace_( org.colSpace_ ),
66 scp_( colSpace_ ),
67 preconditioner_( org.preconditioner_ ),
68 averageCommTime_( org.averageCommTime_ )
69 {}
70
73 ISTLParallelMatrixAdapterInterface ( MatrixType &matrix,
74 const RowSpaceType &rowSpace,
75 const ColSpaceType &colSpace,
76 const PreconditionAdapterType& precon )
77 : matrix_( matrix ),
78 rowSpace_( rowSpace ),
79 colSpace_( colSpace ),
80 scp_( colSpace_ ),
81 preconditioner_( precon ),
82 averageCommTime_( 0 )
83 {}
84
86 virtual double averageCommTime() const { return averageCommTime_ ; }
87
89 virtual PreconditionAdapterType &preconditionAdapter() { return preconditioner_; }
90
92 virtual const PreconditionAdapterType &preconditionAdapter() const { return preconditioner_; }
93
95 virtual ParallelScalarProductType &scp() { return scp_; }
96
98 void apply ( const X &x, Y &y ) const override
99 {
100 DUNE_THROW(NotImplemented,"interface method apply not overloaded!");
101 }
102
104 void applyscaleadd ( field_type alpha, const X &x, Y &y) const override
105 {
106 DUNE_THROW(NotImplemented,"interface method applyscaleadd not overloaded!");
107 }
108
110 virtual double residuum(const Y& rhs, X& x) const
111 {
112 DUNE_THROW(NotImplemented,"interface method residuum not overloaded!");
113 return 0.0;
114 }
115
117 const MatrixType& getmat () const override { return matrix_; }
118
119#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
120 SolverCategory::Category category () const override { return SolverCategory::sequential; }
121#endif // #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
122
123 protected:
124 void communicate( Y &y ) const
125 {
126 if( rowSpace_.grid().comm().size() <= 1 )
127 return;
128
129 Dune::Timer commTime;
130 ColumnDiscreteFunctionType tmp( "LagrangeParallelMatrixAdapter::communicate", colSpace_, y );
131 colSpace_.communicate( tmp );
132 averageCommTime_ += commTime.elapsed();
133 }
134
135 MatrixType &matrix_;
136 const RowSpaceType &rowSpace_;
137 const ColSpaceType &colSpace_;
138 mutable ParallelScalarProductType scp_;
139 PreconditionAdapterType preconditioner_;
140 mutable double averageCommTime_;
141 };
142
143 template <class MatrixImp>
144 class LagrangeParallelMatrixAdapter
145 : public ISTLParallelMatrixAdapterInterface< MatrixImp >
146 {
147 typedef LagrangeParallelMatrixAdapter< MatrixImp > ThisType;
148 typedef ISTLParallelMatrixAdapterInterface< MatrixImp > BaseType;
149 public:
150 typedef MatrixImp MatrixType;
151 typedef Fem::PreconditionerWrapper<MatrixType> PreconditionAdapterType;
152
153 typedef typename MatrixType :: RowDiscreteFunctionType RowDiscreteFunctionType;
154 typedef typename MatrixType :: ColDiscreteFunctionType ColumnDiscreteFunctionType;
155
156 typedef typename RowDiscreteFunctionType :: DiscreteFunctionSpaceType RowSpaceType;
157
158 typedef typename ColumnDiscreteFunctionType :: DiscreteFunctionSpaceType ColSpaceType;
159 typedef Fem::ParallelScalarProduct<ColumnDiscreteFunctionType> ParallelScalarProductType;
160
161 typedef typename RowDiscreteFunctionType :: DofStorageType X;
162 typedef typename ColumnDiscreteFunctionType :: DofStorageType Y;
163
165 typedef MatrixType matrix_type;
166 typedef X domain_type;
167 typedef Y range_type;
168 typedef typename X::field_type field_type;
169
170#if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
172 enum { category=SolverCategory::sequential };
173#endif // #if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
174
175 protected:
176 using BaseType :: matrix_;
177 using BaseType :: scp_ ;
178 using BaseType :: colSpace_ ;
179 using BaseType :: rowSpace_ ;
180 using BaseType :: averageCommTime_;
181
182 public:
184 LagrangeParallelMatrixAdapter ( const LagrangeParallelMatrixAdapter &org )
185 : BaseType( org )
186 {}
187
190 LagrangeParallelMatrixAdapter ( MatrixType &matrix,
191 const RowSpaceType &rowSpace,
192 const ColSpaceType &colSpace,
193 const PreconditionAdapterType& precon )
194 : BaseType( matrix, rowSpace, colSpace, precon )
195 {}
196
198 void apply ( const X &x, Y &y ) const override
199 {
200 matrix_.mv( x, y );
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 matrix_.usmv(alpha,x,y);
210 communicate( y );
211 }
212 else
213 {
214 Y tmp( y.size() );
215 apply( x, tmp );
216 y.axpy( alpha, tmp );
217 }
218 }
219
220 double residuum(const Y& rhs, X& x) const override
221 {
222 Y tmp( rhs );
223
224 this->apply(x,tmp);
225 tmp -= rhs;
226
227 // return global sum of residuum
228 return scp_.norm(tmp);
229 }
230
231#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
232 SolverCategory::Category category () const override { return SolverCategory::sequential; }
233#endif // #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
234
235 protected:
236 void communicate( Y &y ) const
237 {
238 if( rowSpace_.grid().comm().size() <= 1 )
239 return;
240
241 Dune::Timer commTime;
242 ColumnDiscreteFunctionType tmp( "LagrangeParallelMatrixAdapter::communicate", colSpace_, y );
243 colSpace_.communicate( tmp );
244 averageCommTime_ += commTime.elapsed();
245 }
246 };
247
248 template <class MatrixImp>
249 class DGParallelMatrixAdapter
250 : public ISTLParallelMatrixAdapterInterface< MatrixImp >
251 {
252 typedef DGParallelMatrixAdapter< MatrixImp > ThisType ;
253 typedef ISTLParallelMatrixAdapterInterface< MatrixImp > BaseType;
254 public:
255 typedef MatrixImp MatrixType;
256 typedef Fem::PreconditionerWrapper<MatrixType> PreconditionAdapterType;
257
258 typedef typename MatrixType :: RowDiscreteFunctionType RowDiscreteFunctionType;
259 typedef typename MatrixType :: ColDiscreteFunctionType ColumnDiscreteFunctionType;
260
261 typedef typename RowDiscreteFunctionType :: DiscreteFunctionSpaceType RowSpaceType;
262
263 typedef typename ColumnDiscreteFunctionType :: DiscreteFunctionSpaceType ColSpaceType;
264 typedef Fem::ParallelScalarProduct<ColumnDiscreteFunctionType> ParallelScalarProductType;
265
266 typedef typename RowDiscreteFunctionType :: DofStorageType X;
267 typedef typename ColumnDiscreteFunctionType :: DofStorageType Y;
268
270 typedef MatrixType matrix_type;
271 typedef X domain_type;
272 typedef Y range_type;
273 typedef typename X::field_type field_type;
274
275#if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
277 enum { category=SolverCategory::sequential };
278#endif // #if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
279
280 protected:
281 using BaseType :: matrix_;
282 using BaseType :: scp_;
283 using BaseType :: rowSpace_;
284 using BaseType :: colSpace_;
285 using BaseType :: averageCommTime_;
286
287 public:
289 DGParallelMatrixAdapter (const DGParallelMatrixAdapter& org)
290 : BaseType( org )
291 {}
292
294 DGParallelMatrixAdapter (MatrixType& A,
295 const RowSpaceType& rowSpace,
296 const ColSpaceType& colSpace,
297 const PreconditionAdapterType& precon )
298 : BaseType( A, rowSpace, colSpace, precon )
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 matrix_.mv(x,y);
309
310 // delete non-interior
311 scp_.deleteNonInterior( y );
312 }
313
315 void applyscaleadd (field_type alpha, const X& x, Y& y) const override
316 {
317 // exchange data first
318 communicate( x );
319
320 // apply matrix
321 matrix_.usmv(alpha,x,y);
322
323 // delete non-interior
324 scp_.deleteNonInterior( y );
325 }
326
327 double residuum(const Y& rhs, X& x) const override
328 {
329 // exchange data
330 communicate( x );
331
332 typedef typename ParallelScalarProductType :: AuxiliaryDofsType AuxiliaryDofsType;
333 const AuxiliaryDofsType& auxiliaryDofs = scp_.auxiliaryDofs();
334
335 typedef typename Y :: block_type LittleBlockVectorType;
336 LittleBlockVectorType tmp;
337 double res = 0.0;
338
339 // counter for rows
340 int i = 0;
341 const int auxiliarySize = auxiliaryDofs.size();
342 for(int auxiliary = 0; auxiliary<auxiliarySize; ++auxiliary)
343 {
344 const int nextAuxiliary = auxiliaryDofs[auxiliary];
345 for(; i<nextAuxiliary; ++i)
346 {
347 tmp = 0;
348 // get row
349 typedef typename MatrixType :: row_type row_type;
350
351 const row_type& row = matrix_[i];
352 // multiply with row
353 typedef typename MatrixType :: ConstColIterator ConstColIterator;
354 ConstColIterator endj = row.end();
355 for (ConstColIterator j = row.begin(); j!=endj; ++j)
356 {
357 (*j).umv(x[j.index()], tmp);
358 }
359
360 // substract right hand side
361 tmp -= rhs[i];
362
363 // add scalar product
364 res += tmp.two_norm2();
365 }
366 ++i;
367 }
368
369 res = rowSpace_.grid().comm().sum( res );
370 // return global sum of residuum
371 return std::sqrt( res );
372 }
373
374#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
375 SolverCategory::Category category () const override { return SolverCategory::sequential; }
376#endif // #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
377
378 protected:
379 void communicate(const X& x) const
380 {
381 if( rowSpace_.grid().comm().size() <= 1 ) return ;
382
383 Dune::Timer commTime;
384
385 // create temporary discrete function object
386 RowDiscreteFunctionType tmp ("DGParallelMatrixAdapter::communicate",
387 rowSpace_, x );
388
389 // exchange data by copying
390 rowSpace_.communicate( tmp );
391
392 // accumulate communication time
393 averageCommTime_ += commTime.elapsed();
394 }
395 };
396
397 } // end namespace Fem
398} // end namespace Dune
399
400#endif // #if HAVE_DUNE_ISTL
401
402#endif // #ifndef DUNE_FEM_ISTLMATRIXADAPTER_HH
403
double sqrt(const Dune::Fem::Double &v)
Definition: double.hh:977
Definition: bindguard.hh:11