dune-pdelab  2.0.0
petscvectorbackend.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_PETSCVECTORBACKEND_HH
4 #define DUNE_PETSCVECTORBACKEND_HH
5 
6 #if HAVE_PETSC
7 
8 #include<vector>
9 
10 #include<dune/common/fvector.hh>
11 
13 #include <petscvec.h>
14 
15 #include "backendselector.hh"
16 
17 namespace Dune {
18  namespace PDELab {
19 
20  class PetscVectorBackend;
21 
22  class PetscVectorContainer
23  {
24  public:
25  typedef PetscScalar ElementType;
26  typedef ElementType E;
27  typedef Vec ContainerType;
28 
29  typedef E block_type;
30  typedef block_type field_type;
31 
32  typedef std::size_t size_type;
33  typedef PetscVectorBackend Backend;
34 
35  template<typename T>
36  PetscVectorContainer (const T& t_)
37  : _data(NULL)
38  , _managed(true)
39  {
40  PETSC_CALL(VecCreate(PETSC_COMM_SELF,&_v));
41  PETSC_CALL(VecSetSizes(_v,t_.globalSize(),PETSC_DECIDE));
42  PETSC_CALL(VecSetType(_v,VECSEQ));
43  }
44 
45  template<typename T>
46  PetscVectorContainer (const T& t_, const E& e)
47  : _data(NULL)
48  , _managed(true)
49  {
50  PETSC_CALL(VecCreate(PETSC_COMM_SELF,&_v));
51  PETSC_CALL(VecSetSizes(_v,t_.globalSize(),PETSC_DECIDE));
52  PETSC_CALL(VecSetType(_v,VECSEQ));
53  this->operator=(e);
54  }
55 
56  PetscVectorContainer (const PetscVectorContainer& rhs)
57  : _data(NULL)
58  , _managed(true)
59  {
60  rhs.checkin();
61  PETSC_CALL(VecDuplicate(rhs._v,&_v));
62  PETSC_CALL(VecCopy(rhs._v,_v));
63  }
64 
65  PetscVectorContainer(Vec vec, bool managed = true)
66  : _v(vec)
67  , _data(NULL)
68  , _managed(managed)
69  {
70  }
71 
72  PetscVectorContainer& operator= (const PetscVectorContainer& rhs)
73  {
74  checkin();
75  rhs.checkin();
76  PETSC_CALL(VecCopy(rhs._v,_v));
77  _data = NULL;
78  }
79 
80  ~PetscVectorContainer()
81  {
82  checkin();
83  if (_managed)
84  {
85 #if DUNE_PETSC_NEWER(3,2,0)
86  PETSC_CALL(VecDestroy(&_v));
87 #else
88  PETSC_CALL(VecDestroy(_v));
89 #endif
90  }
91  }
92 
93  size_type N() const
94  {
95  int size = 0;
96  PETSC_CALL(VecGetLocalSize(_v,&size));
97  return size;
98  }
99 
100 
101  PetscVectorContainer& operator= (const E& e)
102  {
103  checkin();
104  PETSC_CALL(VecSet(_v,e));
105  return *this;
106  }
107 
108  PetscVectorContainer& operator*= (const E& e)
109  {
110  checkin();
111  PETSC_CALL(VecScale(_v,e));
112  return *this;
113  }
114 
115 
116  PetscVectorContainer& operator+= (const E& e)
117  {
118  checkin();
119  PETSC_CALL(VecShift(_v,e));
120  return *this;
121  }
122 
123  PetscVectorContainer& operator+= (const PetscVectorContainer& rhs)
124  {
125  return axpy(1.0,rhs);
126  }
127 
128  PetscVectorContainer& operator-= (const PetscVectorContainer& rhs)
129  {
130  return axpy(-1.0,rhs);
131  }
132 
133  field_type& operator[](std::size_t i)
134  {
135  checkout();
136  return _data[i];
137  }
138 
139  const field_type& operator[](std::size_t i) const
140  {
141  checkout();
142  return _data[i];
143  }
144 
145  E two_norm() const
146  {
147  checkin();
148  double norm = 0;
149  PETSC_CALL(VecNorm(_v,NORM_2,&norm));
150  return norm;
151  }
152 
153  E one_norm() const
154  {
155  checkin();
156  double norm = 0;
157  PETSC_CALL(VecNorm(_v,NORM_1,&norm));
158  return norm;
159  }
160 
161  E infinity_norm() const
162  {
163  checkin();
164  double norm = 0;
165  PETSC_CALL(VecNorm(_v,NORM_INFINITY,&norm));
166  return norm;
167  }
168 
169  E dot(const PetscVectorContainer& y) const
170  {
171  checkin();
172  PetscScalar dotproduct = 0;
173  PETSC_CALL(VecDot(_v,y._v,&dotproduct));
174  return dotproduct;
175  }
176 
177  E operator*(const PetscVectorContainer& y) const
178  {
179  checkin();
180  PetscScalar dotproduct = 0;
181  PETSC_CALL(VecTDot(_v,y._v,&dotproduct));
182  return dotproduct;
183  }
184 
185 
186  PetscVectorContainer& axpy(const E& a, const PetscVectorContainer& y)
187  {
188  checkin();
189  y.checkin();
190  PETSC_CALL(VecAXPY(_v,a,y._v));
191  return *this;
192  }
193 
194  // for debugging and AMG access
195  ContainerType& base ()
196  {
197  checkin();
198  return _v;
199  }
200 
201  const ContainerType& base () const
202  {
203  checkin();
204  return _v;
205  }
206 
207  operator ContainerType&()
208  {
209  checkin();
210  return _v;
211  }
212 
213  operator const ContainerType&() const
214  {
215  checkin();
216  return _v;
217  }
218 
219  /*
220  iterator begin()
221  {
222  return container.begin();
223  }
224 
225 
226  const_iterator begin() const
227  {
228  return container.begin();
229  }
230 
231  iterator end()
232  {
233  return container.end();
234  }
235 
236 
237  const_iterator end() const
238  {
239  return container.end();
240  }
241  */
242 
243  size_t flatsize() const
244  {
245  return N();
246  }
247 
248  template<typename X>
249  void std_copy_to (std::vector<X>& x) const
250  {
251  size_t n = flatsize();
252  x.resize(n);
253  checkout();
254  for (size_t i=0; i<n; i++)
255  x[i] = _data[i];
256  }
257 
258  template<typename X>
259  void std_copy_from (const std::vector<X>& x)
260  {
261  //test if x has the same size as the container
262  assert (x.size() == flatsize());
263  checkout();
264  for (size_t i=0; i<flatsize(); i++)
265  _data[i] = x[i];
266  }
267 
268  private:
269  Vec _v;
270  mutable E* _data;
271  const bool _managed;
272 
273  void checkin() const
274  {
275  if (_data)
276  {
277  PETSC_CALL(VecRestoreArray(_v,&_data));
278  _data = NULL;
279  }
280  }
281 
282  void checkout() const
283  {
284  if (!_data)
285  {
286  PETSC_CALL(VecGetArray(_v,&_data));
287  }
288  }
289 
290 
291  };
292 
293 
295  class PetscVectorBackend
296  {
297  public:
298  // extract type of container element
299  template<class C>
300  struct Value
301  {
302  typedef typename C::ElementType Type;
303  };
304 
306  typedef typename std::size_t size_type;
307 
308  // get const_reference to container element
309  // note: this method does not depend on T!
310  template<typename C>
311  static const typename C::field_type& access (const C& c, size_type i)
312  {
313  return c[i];
314  }
315 
316  // get non const_reference to container element
317  // note: this method does not depend on T!
318  template<typename C>
319  static typename C::field_type& access (C& c, size_type i)
320  {
321  return c[i];
322  }
323  };
324 
325  template<typename T, typename E>
326  struct BackendVectorSelectorHelper<PetscVectorBackend,T,E>
327  {
328  typedef PetscVectorContainer Type;
329  };
330 
331 
332 
333  } // namespace PDELab
334 } // namespace Dune
335 
336 #endif // HAVE_PETSC
337 
338 #endif // DUNE_PETSCVECTORBACKEND_HH
const E & e
Definition: interpolate.hh:172