Reference documentation for deal.II version 8.1.0
petsc_vector.h
1 // ---------------------------------------------------------------------
2 // @f$Id: petsc_vector.h 30036 2013-07-18 16:55:32Z maier @f$
3 //
4 // Copyright (C) 2004 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__petsc_vector_h
18 #define __deal2__petsc_vector_h
19 
20 
21 #include <deal.II/base/config.h>
22 
23 #ifdef DEAL_II_WITH_PETSC
24 
25 # include <deal.II/base/subscriptor.h>
26 # include <deal.II/lac/exceptions.h>
27 # include <deal.II/lac/petsc_vector_base.h>
28 # include <deal.II/lac/petsc_parallel_vector.h>
29 # include <deal.II/lac/vector.h>
30 
32 
33 
34 namespace PETScWrappers
35 {
52  class Vector : public VectorBase
53  {
54  public:
59 
73  static const bool supports_distributed_data = false;
74 
79  Vector ();
80 
95  explicit Vector (const size_type n);
96 
103  template <typename Number>
104  explicit Vector (const ::Vector<Number> &v);
105 
116  explicit Vector (const Vec &v);
117 
122  Vector (const Vector &v);
123 
138  explicit Vector (const MPI::Vector &v);
139 
144  Vector &operator = (const Vector &v);
145 
159  Vector &operator = (const MPI::Vector &v);
160 
180  Vector &operator = (const PetscScalar s);
181 
188  template <typename number>
189  Vector &operator = (const ::Vector<number> &v);
190 
208  void reinit (const size_type N,
209  const bool fast = false);
210 
221  void reinit (const Vector &v,
222  const bool fast = false);
223 
224  protected:
232  void create_vector (const size_type n);
233  };
234 
237 // ------------------ template and inline functions -------------
238 
239 
248  inline
249  void swap (Vector &u, Vector &v)
250  {
251  u.swap (v);
252  }
253 
254 
255 #ifndef DOXYGEN
256 
257  template <typename number>
258  Vector::Vector (const ::Vector<number> &v)
259  {
260  Vector::create_vector (v.size());
261 
262  *this = v;
263  }
264 
265 
266 
267  inline
268  Vector::Vector (const Vec &v)
269  :
270  VectorBase(v)
271  {}
272 
273 
274  inline
275  Vector &
276  Vector::operator = (const PetscScalar s)
277  {
279 
280  return *this;
281  }
282 
283 
284  inline
285  Vector &
286  Vector::operator = (const Vector &v)
287  {
288  // if the vectors have different sizes,
289  // then first resize the present one
290  if (size() != v.size())
291  reinit (v.size(), true);
292 
293  const int ierr = VecCopy (v.vector, vector);
294  AssertThrow (ierr == 0, ExcPETScError(ierr));
295 
296  return *this;
297  }
298 
299 
300 
301  inline
302  Vector &
304  {
305  int ierr;
306  if (attained_ownership)
307  {
308  // the petsc function we call wants to
309  // generate the vector itself, so destroy
310  // the old one first
311 #if DEAL_II_PETSC_VERSION_LT(3,2,0)
312  ierr = VecDestroy (vector);
313 #else
314  ierr = VecDestroy (&vector);
315 #endif
316  AssertThrow (ierr == 0, ExcPETScError(ierr));
317  }
318 
319  attained_ownership = true;
320 
321  // then do the gather
322  // operation. <rant>petsc has changed its
323  // interface again, and replaced a single
324  // function call by several calls that
325  // are hard to understand. gets me all
326  // annoyed at their development
327  // model</rant>
328 #if DEAL_II_PETSC_VERSION_LT(2,2,0)
329  ierr = VecConvertMPIToSeqAll (static_cast<const Vec &>(v),
330  &vector);
331  AssertThrow (ierr == 0, ExcPETScError(ierr));
332 
333 #else
334 
335  VecScatter ctx;
336 
337  ierr = VecScatterCreateToAll (static_cast<const Vec &>(v), &ctx, &vector);
338  AssertThrow (ierr == 0, ExcPETScError(ierr));
339 
340 #if ((PETSC_VERSION_MAJOR == 2) && \
341  ((PETSC_VERSION_MINOR < 3) || \
342  ((PETSC_VERSION_MINOR == 3) && \
343  (PETSC_VERSION_SUBMINOR < 3))))
344  ierr = VecScatterBegin (static_cast<const Vec &>(v), vector,
345  INSERT_VALUES, SCATTER_FORWARD, ctx);
346  AssertThrow (ierr == 0, ExcPETScError(ierr));
347 
348  ierr = VecScatterEnd (static_cast<const Vec &>(v), vector,
349  INSERT_VALUES, SCATTER_FORWARD, ctx);
350 
351 #else
352 
353  ierr = VecScatterBegin (ctx,static_cast<const Vec &>(v), vector,
354  INSERT_VALUES, SCATTER_FORWARD);
355  AssertThrow (ierr == 0, ExcPETScError(ierr));
356 
357  ierr = VecScatterEnd (ctx, static_cast<const Vec &>(v), vector,
358  INSERT_VALUES, SCATTER_FORWARD);
359 
360 #endif
361  AssertThrow (ierr == 0, ExcPETScError(ierr));
362 
363 #if DEAL_II_PETSC_VERSION_LT(3,2,0)
364  ierr = VecScatterDestroy (ctx);
365 #else
366  ierr = VecScatterDestroy (&ctx);
367 #endif
368  AssertThrow (ierr == 0, ExcPETScError(ierr));
369 #endif
370 
371  return *this;
372  }
373 
374 
375 
376  template <typename number>
377  inline
378  Vector &
379  Vector::operator = (const ::Vector<number> &v)
380  {
381  reinit (v.size(), true);
382  // the following isn't necessarily fast,
383  // but this is due to the fact that PETSc
384  // doesn't offer an inlined access
385  // operator.
386  //
387  // if someone wants to contribute some
388  // code: to make this code faster, one
389  // could either first convert all values
390  // to PetscScalar, and then set them all
391  // at once using VecSetValues. This has
392  // the drawback that it could take quite
393  // some memory, if the vector is large,
394  // and it would in addition allocate
395  // memory on the heap, which is
396  // expensive. an alternative would be to
397  // split the vector into chunks of, say,
398  // 128 elements, convert a chunk at a
399  // time and set it in the output vector
400  // using VecSetValues. since 128 elements
401  // is small enough, this could easily be
402  // allocated on the stack (as a local
403  // variable) which would make the whole
404  // thing much more efficient.
405  //
406  // a second way to make things faster is
407  // for the special case that
408  // number==PetscScalar. we could then
409  // declare a specialization of this
410  // template, and omit the conversion. the
411  // problem with this is that the best we
412  // can do is to use VecSetValues, but
413  // this isn't very efficient either: it
414  // wants to see an array of indices,
415  // which in this case a) again takes up a
416  // whole lot of memory on the heap, and
417  // b) is totally dumb since its content
418  // would simply be the sequence
419  // 0,1,2,3,...,n. the best of all worlds
420  // would probably be a function in Petsc
421  // that would take a pointer to an array
422  // of PetscScalar values and simply copy
423  // n elements verbatim into the vector...
424  for (size_type i=0; i<v.size(); ++i)
425  (*this)(i) = v(i);
426 
427  compress (::VectorOperation::insert);
428 
429  return *this;
430  }
431 #endif // DOXYGEN
432 
433 }
434 
435 DEAL_II_NAMESPACE_CLOSE
436 
437 #endif // DEAL_II_WITH_PETSC
438 
439 /*---------------------------- petsc_vector.h ---------------------------*/
440 
441 #endif
442 /*---------------------------- petsc_vector.h ---------------------------*/
types::global_dof_index size_type
Definition: petsc_vector.h:58
#define AssertThrow(cond, exc)
Definition: exceptions.h:362
Vector & operator=(const Vector &v)
VectorBase & operator=(const PetscScalar s)
static const bool supports_distributed_data
Definition: petsc_vector.h:73
void compress() DEAL_II_DEPRECATED
size_type size() const
unsigned int global_dof_index
Definition: types.h:100
void reinit(const size_type N, const bool fast=false)
void create_vector(const size_type n)
void swap(Vector &u, Vector &v)
Definition: petsc_vector.h:249