dune-pdelab  2.0.0
petscnestedvectorbackend.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_PETSCNESTEDVECTORBACKEND_HH
4 #define DUNE_PETSCNESTEDVECTORBACKEND_HH
5 
6 #if HAVE_PETSC
7 
8 #include<vector>
9 
10 #include<dune/common/fvector.hh>
11 
14 
15 #include <petscvec.h>
16 
17 #include "backendselector.hh"
18 
19 namespace Dune {
20  namespace PDELab {
21 
22  class PetscNestedVectorBackend;
23 
24  class PetscNestedVectorContainer
25  {
26 
27  typedef std::vector<shared_ptr<PetscVectorContainer> > ChildContainer;
28 
29  public:
30  typedef PetscScalar ElementType;
31  typedef ElementType E;
32  typedef Vec ContainerType;
33 
34  typedef E block_type;
35  typedef block_type field_type;
36 
37  typedef std::size_t size_type;
38  typedef PetscVectorBackend Backend;
39 
40  struct create_child_vectors
41  : public TypeTree::DirectChildrenVisitor
42  , public TypeTree::DynamicTraversal
43  {
44  create_child_vectors(ChildContainer& children)
45  : _children(children)
46  {}
47 
48  template<typename T, typename Child, typename TreePath, typename ChildIndex>
49  void beforeChild(T&& t, Child&& child, TreePath treePath, ChildIndex childIndex)
50  {
51  _children[childIndex] = make_shared<PetscVectorContainer>(child);
52  }
53 
54  ChildContainer& _children;
55  };
56 
57  template<typename GFS>
58  PetscNestedVectorContainer (const GFS& gfs)
59  : _data(gfs.globalSize())
60  , _sub_data(GFS::CHILDREN)
61  , _checkedIn(true)
62  {
63  const size_t _blockCount = GFS::CHILDREN;
64  ChildContainer children(_blockCount);
65  create_child_vectors ccv(children);
66  TypeTree::applyToTree(gfs,ccv);
67  Vec child_array[_blockCount];
68  for (size_t i = 0; i < _blockCount; ++i)
69  child_array[i] = children[i]->_v;
70  PETSC_CALL(VecCreateNest(PETSC_COMM_SELF,_blockCount,PETSC_NULL,child_array,&_v));
71  }
72 
73  template<typename GFS>
74  PetscNestedVectorContainer (const GFS& gfs, const E& e)
75  : _data(gfs.globalSize())
76  , _sub_data(GFS::CHILDREN)
77  , _checkedIn(true)
78  {
79  const size_t _blockCount = GFS::CHILDREN;
80  ChildContainer children(_blockCount);
81  create_child_vectors ccv(children);
82  TypeTree::applyToTree(gfs,ccv);
83  Vec child_array[_blockCount];
84  for (size_t i = 0; i < _blockCount; ++i)
85  child_array[i] = *(children[i]);
86  PETSC_CALL(VecCreateNest(PETSC_COMM_SELF,_blockCount,PETSC_NULL,child_array,&_v));
87  this->operator=(e);
88  }
89 
90  PetscNestedVectorContainer (const PetscNestedVectorContainer& rhs)
91  : _data(rhs._data.size())
92  , _sub_data(rhs._sub_data.size())
93  , _checkedIn(true)
94  {
95  rhs.checkin();
96  PETSC_CALL(VecDuplicate(rhs._v,&_v));
97  PETSC_CALL(VecCopy(rhs._v,_v));
98  }
99 
100  PetscNestedVectorContainer& operator= (const PetscNestedVectorContainer& rhs)
101  {
102  checkin();
103  rhs.checkin();
104  PETSC_CALL(VecCopy(rhs._v,_v));
105  _checkedIn = true;
106  }
107 
108  ~PetscNestedVectorContainer()
109  {
110  checkin();
111  PETSC_CALL(VecDestroy(&_v));
112  }
113 
114  size_type N() const
115  {
116  int size = 0;
117  PETSC_CALL(VecGetLocalSize(_v,&size));
118  return size;
119  }
120 
121 
122  PetscNestedVectorContainer& operator= (const E& e)
123  {
124  checkin();
125  PETSC_CALL(VecSet(_v,e));
126  return *this;
127  }
128 
129  PetscNestedVectorContainer& operator*= (const E& e)
130  {
131  checkin();
132  PETSC_CALL(VecScale(_v,e));
133  return *this;
134  }
135 
136 
137  PetscNestedVectorContainer& operator+= (const E& e)
138  {
139  checkin();
140  PETSC_CALL(VecShift(_v,e));
141  return *this;
142  }
143 
144  PetscNestedVectorContainer& operator+= (const PetscNestedVectorContainer& rhs)
145  {
146  return axpy(1.0,rhs);
147  }
148 
149  PetscNestedVectorContainer& operator-= (const PetscNestedVectorContainer& rhs)
150  {
151  return axpy(-1.0,rhs);
152  }
153 
154  field_type& operator[](std::size_t i)
155  {
156  checkout();
157  return _data[i];
158  }
159 
160  const field_type& operator[](std::size_t i) const
161  {
162  checkout();
163  return _data[i];
164  }
165 
166  E two_norm() const
167  {
168  checkin();
169  double norm = 0;
170  PETSC_CALL(VecNorm(_v,NORM_2,&norm));
171  return norm;
172  }
173 
174  E one_norm() const
175  {
176  checkin();
177  double norm = 0;
178  PETSC_CALL(VecNorm(_v,NORM_1,&norm));
179  return norm;
180  }
181 
182  E infinity_norm() const
183  {
184  checkin();
185  double norm = 0;
186  PETSC_CALL(VecNorm(_v,NORM_INFINITY,&norm));
187  return norm;
188  }
189 
190  E operator*(const PetscNestedVectorContainer& y) const
191  {
192  checkin();
193  PetscScalar dotproduct = 0;
194  PETSC_CALL(VecTDot(_v,y._v,&dotproduct));
195  return dotproduct;
196  }
197 
198  E dot(const PetscNestedVectorContainer& y) const
199  {
200  checkin();
201  PetscScalar dotproduct = 0;
202  PETSC_CALL(VecDot(_v,y._v,&dotproduct));
203  return dotproduct;
204  }
205 
206 
207  PetscNestedVectorContainer& axpy(const E& a, const PetscNestedVectorContainer& y)
208  {
209  checkin();
210  y.checkin();
211  PETSC_CALL(VecAXPY(_v,a,y._v));
212  return *this;
213  }
214 
215  // for debugging and AMG access
216  ContainerType& base ()
217  {
218  checkin();
219  return _v;
220  }
221 
222  const ContainerType& base () const
223  {
224  checkin();
225  return _v;
226  }
227 
228  operator ContainerType&()
229  {
230  checkin();
231  return _v;
232  }
233 
234  operator const ContainerType&() const
235  {
236  checkin();
237  return _v;
238  }
239 
240  /*
241  iterator begin()
242  {
243  return container.begin();
244  }
245 
246 
247  const_iterator begin() const
248  {
249  return container.begin();
250  }
251 
252  iterator end()
253  {
254  return container.end();
255  }
256 
257 
258  const_iterator end() const
259  {
260  return container.end();
261  }
262  */
263 
264  size_t flatsize() const
265  {
266  return N();
267  }
268 
269  template<typename X>
270  void std_copy_to (std::vector<X>& x) const
271  {
272  checkout();
273  std::copy(_data.begin(),_data.end(),x.begin());
274  }
275 
276  template<typename X>
277  void std_copy_from (const std::vector<X>& x)
278  {
279  checkout();
280  std::copy(x.begin(),x.end(),_data.begin());
281  }
282 
283  PetscVectorContainer subVector(size_type i)
284  {
285  checkin();
286  Vec sub_vec;
287  PETSC_CALL(VecNestGetSubVec(_v,i,&sub_vec));
288  return PetscVectorContainer(sub_vec,false);
289  }
290 
291  private:
292  Vec _v;
293  mutable std::vector<PetscScalar> _data;
294  mutable std::vector<PetscScalar*> _sub_data;
295  mutable bool _checkedIn;
296 
297  void checkin() const
298  {
299  if (!_checkedIn)
300  {
301  int N;
302  Vec* sv;
303  PETSC_CALL(VecNestGetSubVecs(_v,&N,&sv));
304  auto it = _data.begin();
305  for (int n = 0; n < N; ++n)
306  {
307  int CN;
308  PETSC_CALL(VecGetLocalSize(sv[n],&CN));
309  auto endit = it + CN;
310  for(auto childit = _sub_data[n]; it != endit; ++it, ++childit)
311  *childit = *it;
312  PETSC_CALL(VecRestoreArray(sv[n],&(_sub_data[n])));
313  }
314  _checkedIn = true;
315  }
316  }
317 
318  void checkout() const
319  {
320  if (_checkedIn)
321  {
322  int N;
323  Vec* sv;
324  PETSC_CALL(VecNestGetSubVecs(_v,&N,&sv));
325  auto it = _data.begin();
326  for (int n = 0; n < N; ++n)
327  {
328  int CN;
329  PETSC_CALL(VecGetArray(sv[n],&(_sub_data[n])));
330  PETSC_CALL(VecGetLocalSize(sv[n],&CN));
331  auto endit = it + CN;
332  for(auto childit = _sub_data[n]; it != endit; ++it, ++childit)
333  *it = *childit;
334  }
335  _checkedIn = false;
336  }
337  }
338 
339 
340  };
341 
342 
344  class PetscNestedVectorBackend
345  {
346  public:
347  /* enum{
349  BlockSize = BLOCKSIZE
350  };*/
351 
352  //export Matrix Backend Type
353  //typedef ISTLBCRSMatrixBackend<BLOCKSIZE,BLOCKSIZE> MatrixBackend;
354 
356 
357  // extract type of container element
358  template<class C>
359  struct Value
360  {
361  typedef typename C::ElementType Type;
362  };
363 
365  typedef typename std::size_t size_type;
366 
367  // get const_reference to container element
368  // note: this method does not depend on T!
369  template<typename C>
370  static const typename C::field_type& access (const C& c, size_type i)
371  {
372  return c[i];
373  }
374 
375  // get non const_reference to container element
376  // note: this method does not depend on T!
377  template<typename C>
378  static typename C::field_type& access (C& c, size_type i)
379  {
380  return c[i];
381  }
382  };
383 
384  template<typename T, typename E>
385  struct BackendVectorSelectorHelper<PetscNestedVectorBackend,T,E>
386  {
387  typedef PetscNestedVectorContainer Type;
388  };
389 
390 
391 
392  } // namespace PDELab
393 } // namespace Dune
394 
395 #endif // HAVE_PETSC
396 
397 #endif // DUNE_PETSCNESTEDVECTORBACKEND_HH
const E & e
Definition: interpolate.hh:172