1#ifndef DUNE_FEM_ISTLMATRIXADAPTER_HH
2#define DUNE_FEM_ISTLMATRIXADAPTER_HH
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>
23 template <
class Matrix>
24 class PreconditionerWrapper ;
27 template <
class MatrixImp>
28 class ISTLParallelMatrixAdapterInterface
29 :
public AssembledLinearOperator< MatrixImp,
30 typename MatrixImp :: RowBlockVectorType,
31 typename MatrixImp :: ColBlockVectorType>
33 typedef ISTLParallelMatrixAdapterInterface< MatrixImp > ThisType;
35 typedef MatrixImp MatrixType;
36 typedef Fem::PreconditionerWrapper< MatrixType > PreconditionAdapterType;
38 typedef typename MatrixType :: RowDiscreteFunctionType RowDiscreteFunctionType;
39 typedef typename MatrixType :: ColDiscreteFunctionType ColumnDiscreteFunctionType;
41 typedef typename RowDiscreteFunctionType :: DiscreteFunctionSpaceType RowSpaceType;
43 typedef typename ColumnDiscreteFunctionType :: DiscreteFunctionSpaceType ColSpaceType;
44 typedef Fem::ParallelScalarProduct<ColumnDiscreteFunctionType> ParallelScalarProductType;
46 typedef typename RowDiscreteFunctionType :: DofStorageType X;
47 typedef typename ColumnDiscreteFunctionType :: DofStorageType Y;
50 typedef MatrixType matrix_type;
51 typedef X domain_type;
53 typedef typename X::field_type field_type;
57 ISTLParallelMatrixAdapterInterface (
const ISTLParallelMatrixAdapterInterface &org )
58 : matrix_( org.matrix_ ),
59 rowSpace_( org.rowSpace_ ),
60 colSpace_( org.colSpace_ ),
62 preconditioner_( org.preconditioner_ ),
63 averageCommTime_( org.averageCommTime_ ),
64 threading_( org.threading_ )
69 ISTLParallelMatrixAdapterInterface ( MatrixType &matrix,
70 const RowSpaceType &rowSpace,
71 const ColSpaceType &colSpace,
72 const PreconditionAdapterType& precon,
73 const bool threading =
true )
75 rowSpace_( rowSpace ),
76 colSpace_( colSpace ),
78 preconditioner_( precon ),
79 averageCommTime_( 0 ),
80 threading_( threading )
84 virtual double averageCommTime()
const {
return averageCommTime_ ; }
87 virtual PreconditionAdapterType &preconditionAdapter() {
return preconditioner_; }
90 virtual const PreconditionAdapterType &preconditionAdapter()
const {
return preconditioner_; }
93 virtual ParallelScalarProductType &scp() {
return scp_; }
96 void apply (
const X &x, Y &y )
const override
98 DUNE_THROW(NotImplemented,
"interface method apply not overloaded!");
102 void applyscaleadd ( field_type alpha,
const X &x, Y &y)
const override
104 DUNE_THROW(NotImplemented,
"interface method applyscaleadd not overloaded!");
108 virtual double residuum(
const Y& rhs, X& x)
const
110 DUNE_THROW(NotImplemented,
"interface method residuum not overloaded!");
115 const MatrixType& getmat ()
const override {
return matrix_; }
118 bool threading ()
const {
return threading_; }
121 void communicate( Y &y )
const
123 if( rowSpace_.grid().comm().size() <= 1 )
127 ColumnDiscreteFunctionType tmp(
"LagrangeParallelMatrixAdapter::communicate", colSpace_, y );
128 colSpace_.communicate( tmp );
129 averageCommTime_ += commTime.
elapsed();
133 const RowSpaceType &rowSpace_;
134 const ColSpaceType &colSpace_;
135 mutable ParallelScalarProductType scp_;
136 PreconditionAdapterType preconditioner_;
137 mutable double averageCommTime_;
138 const bool threading_;
141 template <
class MatrixImp>
142 class LagrangeParallelMatrixAdapter
143 :
public ISTLParallelMatrixAdapterInterface< MatrixImp >
145 typedef LagrangeParallelMatrixAdapter< MatrixImp > ThisType;
146 typedef ISTLParallelMatrixAdapterInterface< MatrixImp > BaseType;
148 typedef MatrixImp MatrixType;
149 typedef Fem::PreconditionerWrapper<MatrixType> PreconditionAdapterType;
151 typedef typename MatrixType :: RowDiscreteFunctionType RowDiscreteFunctionType;
152 typedef typename MatrixType :: ColDiscreteFunctionType ColumnDiscreteFunctionType;
154 typedef typename RowDiscreteFunctionType :: DiscreteFunctionSpaceType RowSpaceType;
156 typedef typename ColumnDiscreteFunctionType :: DiscreteFunctionSpaceType ColSpaceType;
157 typedef Fem::ParallelScalarProduct<ColumnDiscreteFunctionType> ParallelScalarProductType;
159 typedef typename RowDiscreteFunctionType :: DofStorageType X;
160 typedef typename ColumnDiscreteFunctionType :: DofStorageType Y;
163 typedef MatrixType matrix_type;
164 typedef X domain_type;
165 typedef Y range_type;
166 typedef typename X::field_type field_type;
168 using BaseType :: threading;
171 using BaseType :: matrix_;
172 using BaseType :: scp_ ;
173 using BaseType :: colSpace_ ;
174 using BaseType :: rowSpace_ ;
175 using BaseType :: averageCommTime_;
179 LagrangeParallelMatrixAdapter (
const LagrangeParallelMatrixAdapter &org )
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 )
194 void apply (
const X &x, Y &y )
const override
197 matrix_.mvThreaded(x, y );
205 void applyscaleadd ( field_type alpha,
const X &x, Y &y)
const override
207 if( rowSpace_.grid().comm().size() <= 1 )
210 matrix_.usmvThreaded(alpha, x, y );
212 matrix_.usmv(alpha,x,y);
220 y.axpy( alpha, tmp );
224 double residuum(
const Y& rhs, X& x)
const override
232 return scp_.norm(tmp);
238 void communicate( Y &y )
const
240 if( rowSpace_.grid().comm().size() <= 1 )
244 ColumnDiscreteFunctionType tmp(
"LagrangeParallelMatrixAdapter::communicate", colSpace_, y );
245 colSpace_.communicate( tmp );
246 averageCommTime_ += commTime.
elapsed();
250 template <
class MatrixImp>
251 class DGParallelMatrixAdapter
252 :
public ISTLParallelMatrixAdapterInterface< MatrixImp >
254 typedef DGParallelMatrixAdapter< MatrixImp > ThisType ;
255 typedef ISTLParallelMatrixAdapterInterface< MatrixImp > BaseType;
257 typedef MatrixImp MatrixType;
258 typedef Fem::PreconditionerWrapper<MatrixType> PreconditionAdapterType;
260 typedef typename MatrixType :: RowDiscreteFunctionType RowDiscreteFunctionType;
261 typedef typename MatrixType :: ColDiscreteFunctionType ColumnDiscreteFunctionType;
263 typedef typename RowDiscreteFunctionType :: DiscreteFunctionSpaceType RowSpaceType;
265 typedef typename ColumnDiscreteFunctionType :: DiscreteFunctionSpaceType ColSpaceType;
266 typedef Fem::ParallelScalarProduct<ColumnDiscreteFunctionType> ParallelScalarProductType;
268 typedef typename RowDiscreteFunctionType :: DofStorageType X;
269 typedef typename ColumnDiscreteFunctionType :: DofStorageType Y;
272 typedef MatrixType matrix_type;
273 typedef X domain_type;
274 typedef Y range_type;
275 typedef typename X::field_type field_type;
277 using BaseType :: threading;
280 using BaseType :: matrix_;
281 using BaseType :: scp_;
282 using BaseType :: rowSpace_;
283 using BaseType :: colSpace_;
284 using BaseType :: averageCommTime_;
288 DGParallelMatrixAdapter (
const DGParallelMatrixAdapter& org)
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 )
302 void apply (
const X& x, Y& y)
const override
309 matrix_.mvThreaded( x, y );
314 scp_.deleteNonInterior( y );
318 void applyscaleadd (field_type alpha,
const X& x, Y& y)
const override
325 matrix_.usmvThreaded(alpha, x, y );
327 matrix_.usmv(alpha,x,y);
330 scp_.deleteNonInterior( y );
333 double residuum(
const Y& rhs, X& x)
const override
338 typedef typename ParallelScalarProductType :: AuxiliaryDofsType AuxiliaryDofsType;
339 const AuxiliaryDofsType& auxiliaryDofs = scp_.auxiliaryDofs();
341 typedef typename Y :: block_type LittleBlockVectorType;
342 LittleBlockVectorType tmp;
347 const int auxiliarySize = auxiliaryDofs.size();
348 for(
int auxiliary = 0; auxiliary<auxiliarySize; ++auxiliary)
350 const int nextAuxiliary = auxiliaryDofs[auxiliary];
351 for(; i<nextAuxiliary; ++i)
355 typedef typename MatrixType :: row_type row_type;
357 const row_type& row = matrix_[i];
359 typedef typename MatrixType :: ConstColIterator ConstColIterator;
360 ConstColIterator endj = row.end();
361 for (ConstColIterator j = row.begin(); j!=endj; ++j)
363 (*j).umv(x[j.index()], tmp);
370 res += tmp.two_norm2();
375 res = rowSpace_.grid().comm().sum( res );
377 return std::sqrt( res );
383 void communicate(
const X& x)
const
385 if( rowSpace_.grid().comm().size() <= 1 ) return ;
390 RowDiscreteFunctionType tmp (
"DGParallelMatrixAdapter::communicate",
394 rowSpace_.communicate( tmp );
397 averageCommTime_ += commTime.
elapsed();
A simple stop watch.
Definition: timer.hh:43
double elapsed() const noexcept
Get elapsed user-time from last reset until now/last stop in seconds.
Definition: timer.hh:77
Various macros to work with Dune module version numbers.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
Dune namespace.
Definition: alignedallocator.hh:13
Define general, extensible interface for operators. The available implementation wraps a matrix.
Category
Definition: solvercategory.hh:23
@ sequential
Category for sequential solvers.
Definition: solvercategory.hh:25