dune-fem  2.8-git
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 
18 namespace 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 
Definition: bindguard.hh:11
static double sqrt(const Double &v)
Definition: double.hh:886