Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
petsc_parallel_vector.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/mpi.h>
17 
18 #include <deal.II/lac/petsc_vector.h>
19 
20 #ifdef DEAL_II_WITH_PETSC
21 
22 # include <algorithm>
23 # include <cmath>
24 
25 DEAL_II_NAMESPACE_OPEN
26 
27 namespace PETScWrappers
28 {
29  namespace MPI
30  {
32  : communicator(MPI_COMM_SELF)
33  {
34  // virtual functions called in constructors and destructors never use the
35  // override in a derived class
36  // for clarity be explicit on which function is called
38  }
39 
40 
41 
42  Vector::Vector(const MPI_Comm &communicator,
43  const size_type n,
44  const size_type local_size)
45  : communicator(communicator)
46  {
48  }
49 
50 
51 
52  Vector::Vector(const MPI_Comm & communicator,
53  const VectorBase &v,
54  const size_type local_size)
55  : VectorBase(v)
56  , communicator(communicator)
57  {
58  // In the past (before it was deprecated) this constructor did a
59  // byte-for-byte copy of v. This choice resulted in two problems:
60  // 1. The created vector will have the same size as v, not local size.
61  // 2. Since both the created vector and v maintain ownership of the same
62  // PETSc Vec, both will try to destroy it: this does not make sense.
63  //
64  // For the sake of backwards compatibility, preserve the behavior of the
65  // copy, but correct the ownership bug. Note that in both this (and the
66  // original) implementation local_size is ultimately unused.
67  (void)local_size;
68  }
69 
70 
71 
72  Vector::Vector(const IndexSet &local,
73  const IndexSet &ghost,
74  const MPI_Comm &communicator)
75  : communicator(communicator)
76  {
79 
80  IndexSet ghost_set = ghost;
81  ghost_set.subtract_set(local);
82 
83  Vector::create_vector(local.size(), local.n_elements(), ghost_set);
84  }
85 
86 
87 
89  : VectorBase()
90  , communicator(v.communicator)
91  {
92  if (v.has_ghost_elements())
94  else
96 
97  this->operator=(v);
98  }
99 
100 
101 
102  Vector::Vector(const IndexSet &local, const MPI_Comm &communicator)
103  : communicator(communicator)
104  {
107  Vector::create_vector(local.size(), local.n_elements());
108  }
109 
110 
111 
112  Vector &
114  {
115  // make sure left- and right-hand side of the assignment are
116  // compress()'ed:
118  internal::VectorReference::ExcWrongMode(VectorOperation::unknown,
119  v.last_action));
121  internal::VectorReference::ExcWrongMode(VectorOperation::unknown,
122  last_action));
123 
124  // if the vectors have different sizes,
125  // then first resize the present one
126  if (size() != v.size())
127  {
128  if (v.has_ghost_elements())
130  else
131  reinit(v.communicator, v.size(), v.local_size(), true);
132  }
133 
134  PetscErrorCode ierr = VecCopy(v.vector, vector);
135  AssertThrow(ierr == 0, ExcPETScError(ierr));
136 
137  if (has_ghost_elements())
138  {
139  ierr = VecGhostUpdateBegin(vector, INSERT_VALUES, SCATTER_FORWARD);
140  AssertThrow(ierr == 0, ExcPETScError(ierr));
141  ierr = VecGhostUpdateEnd(vector, INSERT_VALUES, SCATTER_FORWARD);
142  AssertThrow(ierr == 0, ExcPETScError(ierr));
143  }
144  return *this;
145  }
146 
147 
148 
149  void
151  {
152  obtained_ownership = true;
154 
155  create_vector(0, 0);
156  }
157 
158 
159 
160  void
161  Vector::reinit(const MPI_Comm &comm,
162  const size_type n,
163  const size_type local_sz,
164  const bool omit_zeroing_entries)
165  {
166  communicator = comm;
167 
168  // only do something if the sizes
169  // mismatch (may not be true for every proc)
170 
171  int k_global, k = ((size() != n) || (local_size() != local_sz));
172  {
173  const int ierr =
174  MPI_Allreduce(&k, &k_global, 1, MPI_INT, MPI_LOR, communicator);
175  AssertThrowMPI(ierr);
176  }
177 
178  if (k_global || has_ghost_elements())
179  {
180  // FIXME: I'd like to use this here,
181  // but somehow it leads to odd errors
182  // somewhere down the line in some of
183  // the tests:
184  // const PetscErrorCode ierr = VecSetSizes (vector, n, n);
185  // AssertThrow (ierr == 0, ExcPETScError(ierr));
186 
187  // so let's go the slow way:
188 
189  const PetscErrorCode ierr = VecDestroy(&vector);
190  AssertThrow(ierr == 0, ExcPETScError(ierr));
191 
192  create_vector(n, local_sz);
193  }
194 
195  // finally clear the new vector if so
196  // desired
197  if (omit_zeroing_entries == false)
198  *this = 0;
199  }
200 
201 
202 
203  void
204  Vector::reinit(const Vector &v, const bool omit_zeroing_entries)
205  {
206  if (v.has_ghost_elements())
207  {
209  if (!omit_zeroing_entries)
210  {
211  const PetscErrorCode ierr = VecSet(vector, 0.0);
212  AssertThrow(ierr == 0, ExcPETScError(ierr));
213  }
214  }
215  else
216  reinit(v.communicator, v.size(), v.local_size(), omit_zeroing_entries);
217  }
218 
219 
220 
221  void
222  Vector::reinit(const IndexSet &local,
223  const IndexSet &ghost,
224  const MPI_Comm &comm)
225  {
226  const PetscErrorCode ierr = VecDestroy(&vector);
227  AssertThrow(ierr == 0, ExcPETScError(ierr));
228 
229  communicator = comm;
230 
232 
233  IndexSet ghost_set = ghost;
234  ghost_set.subtract_set(local);
235 
236  create_vector(local.size(), local.n_elements(), ghost_set);
237  }
238 
239  void
240  Vector::reinit(const IndexSet &local, const MPI_Comm &comm)
241  {
242  const PetscErrorCode ierr = VecDestroy(&vector);
243  AssertThrow(ierr == 0, ExcPETScError(ierr));
244 
245  communicator = comm;
246 
248  Assert(local.size() > 0, ExcMessage("can not create vector of size 0."));
249  create_vector(local.size(), local.n_elements());
250  }
251 
252 
253  void
254  Vector::create_vector(const size_type n, const size_type local_size)
255  {
256  (void)n;
258  ghosted = false;
259 
260  const PetscErrorCode ierr =
261  VecCreateMPI(communicator, local_size, PETSC_DETERMINE, &vector);
262  AssertThrow(ierr == 0, ExcPETScError(ierr));
263 
264  Assert(size() == n, ExcDimensionMismatch(size(), n));
265  }
266 
267 
268 
269  void
271  const size_type local_size,
272  const IndexSet &ghostnodes)
273  {
274  (void)n;
276  ghosted = true;
277  ghost_indices = ghostnodes;
278 
279  std::vector<size_type> ghostindices;
280  ghostnodes.fill_index_vector(ghostindices);
281 
282  const PetscInt *ptr =
283  (ghostindices.size() > 0 ?
284  reinterpret_cast<const PetscInt *>(ghostindices.data()) :
285  nullptr);
286 
287  PetscErrorCode ierr = VecCreateGhost(communicator,
288  local_size,
289  PETSC_DETERMINE,
290  ghostindices.size(),
291  ptr,
292  &vector);
293  AssertThrow(ierr == 0, ExcPETScError(ierr));
294 
295  Assert(size() == n, ExcDimensionMismatch(size(), n));
296 
297 # if DEBUG
298  {
299  // test ghost allocation in debug mode
300  PetscInt begin, end;
301 
302  ierr = VecGetOwnershipRange(vector, &begin, &end);
303  AssertThrow(ierr == 0, ExcPETScError(ierr));
304 
305  AssertDimension(local_size, static_cast<size_type>(end - begin));
306 
307  Vec l;
308  ierr = VecGhostGetLocalForm(vector, &l);
309  AssertThrow(ierr == 0, ExcPETScError(ierr));
310 
311  PetscInt lsize;
312  ierr = VecGetSize(l, &lsize);
313  AssertThrow(ierr == 0, ExcPETScError(ierr));
314 
315  ierr = VecGhostRestoreLocalForm(vector, &l);
316  AssertThrow(ierr == 0, ExcPETScError(ierr));
317 
318  AssertDimension(lsize,
319  end - begin +
320  static_cast<PetscInt>(ghost_indices.n_elements()));
321  }
322 # endif
323 
324 
325  // in PETSc versions up to 3.5, VecCreateGhost zeroed out the locally
326  // owned vector elements but forgot about the ghost elements. we need to
327  // do this ourselves
328  //
329  // see https://code.google.com/p/dealii/issues/detail?id=233
330 # if DEAL_II_PETSC_VERSION_LT(3, 6, 0)
332  zero.reinit(communicator, this->size(), local_size);
333  *this = zero;
334 # endif
335  }
336 
337 
338 
339  bool
341  {
342  unsigned int has_nonzero = VectorBase::all_zero() ? 0 : 1;
343 # ifdef DEAL_II_WITH_MPI
344  // in parallel, check that the vector
345  // is zero on _all_ processors.
346  unsigned int num_nonzero = Utilities::MPI::sum(has_nonzero, communicator);
347  return num_nonzero == 0;
348 # else
349  return has_nonzero == 0;
350 # endif
351  }
352 
353 
354  void
355  Vector::print(std::ostream & out,
356  const unsigned int precision,
357  const bool scientific,
358  const bool across) const
359  {
360  AssertThrow(out, ExcIO());
361 
362  // get a representation of the vector and
363  // loop over all the elements
364  PetscScalar *val;
365  PetscInt nlocal, istart, iend;
366 
367  PetscErrorCode ierr = VecGetArray(vector, &val);
368  AssertThrow(ierr == 0, ExcPETScError(ierr));
369 
370  ierr = VecGetLocalSize(vector, &nlocal);
371  AssertThrow(ierr == 0, ExcPETScError(ierr));
372 
373  ierr = VecGetOwnershipRange(vector, &istart, &iend);
374  AssertThrow(ierr == 0, ExcPETScError(ierr));
375 
376  // save the state of out stream
377  std::ios::fmtflags old_flags = out.flags();
378  unsigned int old_precision = out.precision(precision);
379 
380  out.precision(precision);
381  if (scientific)
382  out.setf(std::ios::scientific, std::ios::floatfield);
383  else
384  out.setf(std::ios::fixed, std::ios::floatfield);
385 
386  // let each processor produce its output in turn. this requires
387  // synchronizing output between processors using a barrier --
388  // which is clearly slow, but nobody is going to print a whole
389  // matrix this way on a regular basis for production runs, so
390  // the slowness of the barrier doesn't matter
391  for (unsigned int i = 0;
393  i++)
394  {
395  const int mpi_ierr = MPI_Barrier(communicator);
396  AssertThrowMPI(mpi_ierr);
397 
399  {
400  if (across)
401  {
402  out << "[Proc" << i << " " << istart << "-" << iend - 1 << "]"
403  << ' ';
404  for (PetscInt i = 0; i < nlocal; ++i)
405  out << val[i] << ' ';
406  }
407  else
408  {
409  out << "[Proc " << i << " " << istart << "-" << iend - 1
410  << "]" << std::endl;
411  for (PetscInt i = 0; i < nlocal; ++i)
412  out << val[i] << std::endl;
413  }
414  out << std::endl;
415  }
416  }
417  // reset output format
418  out.flags(old_flags);
419  out.precision(old_precision);
420 
421  // restore the representation of the
422  // vector
423  ierr = VecRestoreArray(vector, &val);
424  AssertThrow(ierr == 0, ExcPETScError(ierr));
425 
426  AssertThrow(out, ExcIO());
427  }
428 
429  } // namespace MPI
430 
431 } // namespace PETScWrappers
432 
433 DEAL_II_NAMESPACE_CLOSE
434 
435 #endif // DEAL_II_WITH_PETSC
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
static ::ExceptionBase & ExcIO()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
virtual void clear() override
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
bool is_ascending_and_one_to_one(const MPI_Comm &communicator) const
Definition: index_set.cc:664
size_type size() const
Definition: index_set.h:1600
static ::ExceptionBase & ExcMessage(std::string arg1)
T sum(const T &t, const MPI_Comm &mpi_communicator)
void subtract_set(const IndexSet &other)
Definition: index_set.cc:267
VectorOperation::values last_action
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
types::global_dof_index size_type
Definition: petsc_vector.h:164
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:71
void fill_index_vector(std::vector< size_type > &indices) const
Definition: index_set.cc:505
Vector & operator=(const Vector &v)
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1695
bool has_ghost_elements() const
void reinit(const MPI_Comm &communicator, const size_type N, const size_type local_size, const bool omit_zeroing_entries=false)
IndexSet locally_owned_elements() const
virtual void create_vector(const size_type n, const size_type local_size)
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:82
static ::ExceptionBase & ExcNotImplemented()
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
size_type n_elements() const
Definition: index_set.h:1732