Reference documentation for deal.II version Git 3f1f337db3 2021-10-23 13:19:02 -0600
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
petsc_parallel_vector.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2021 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 
19 
20 #ifdef DEAL_II_WITH_PETSC
21 
22 # include <algorithm>
23 # include <cmath>
24 
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 
43  const size_type n,
45  : communicator(communicator)
46  {
47  Vector::create_vector(n, locally_owned_size);
48  }
49 
50 
51 
53  const VectorBase &v,
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 locally_owned_size is ultimately unused.
67  (void)locally_owned_size;
68  }
69 
70 
71 
72  Vector::Vector(const IndexSet &local,
73  const IndexSet &ghost,
74  const MPI_Comm &communicator)
75  : communicator(communicator)
76  {
77  Assert(local.is_ascending_and_one_to_one(communicator),
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()
91  {
92  if (v.has_ghost_elements())
95  v.ghost_indices);
96  else
98 
99  this->operator=(v);
100  }
101 
102 
103 
105  : communicator(communicator)
106  {
107  Assert(local.is_ascending_and_one_to_one(communicator),
109  Vector::create_vector(local.size(), local.n_elements());
110  }
111 
112 
113 
114  Vector &
116  {
117  // make sure left- and right-hand side of the assignment are
118  // compress()'ed:
120  internal::VectorReference::ExcWrongMode(VectorOperation::unknown,
121  v.last_action));
123  internal::VectorReference::ExcWrongMode(VectorOperation::unknown,
124  last_action));
125 
126  // if the vectors have different sizes,
127  // then first resize the present one
128  if (size() != v.size())
129  {
130  if (v.has_ghost_elements())
132  else
133  reinit(v.communicator, v.size(), v.locally_owned_size(), true);
134  }
135 
136  PetscErrorCode ierr = VecCopy(v.vector, vector);
137  AssertThrow(ierr == 0, ExcPETScError(ierr));
138 
139  if (has_ghost_elements())
140  {
141  ierr = VecGhostUpdateBegin(vector, INSERT_VALUES, SCATTER_FORWARD);
142  AssertThrow(ierr == 0, ExcPETScError(ierr));
143  ierr = VecGhostUpdateEnd(vector, INSERT_VALUES, SCATTER_FORWARD);
144  AssertThrow(ierr == 0, ExcPETScError(ierr));
145  }
146  return *this;
147  }
148 
149 
150 
151  void
153  {
154  obtained_ownership = true;
156 
157  create_vector(0, 0);
158  }
159 
160 
161 
162  void
164  const size_type n,
165  const size_type local_sz,
166  const bool omit_zeroing_entries)
167  {
168  communicator = comm;
169 
170  // only do something if the sizes
171  // mismatch (may not be true for every proc)
172 
173  int k_global, k = ((size() != n) || (locally_owned_size() != local_sz));
174  {
175  const int ierr =
176  MPI_Allreduce(&k, &k_global, 1, MPI_INT, MPI_LOR, communicator);
177  AssertThrowMPI(ierr);
178  }
179 
180  if (k_global || has_ghost_elements())
181  {
182  // FIXME: I'd like to use this here,
183  // but somehow it leads to odd errors
184  // somewhere down the line in some of
185  // the tests:
186  // const PetscErrorCode ierr = VecSetSizes (vector, n, n);
187  // AssertThrow (ierr == 0, ExcPETScError(ierr));
188 
189  // so let's go the slow way:
190 
191  const PetscErrorCode ierr = VecDestroy(&vector);
192  AssertThrow(ierr == 0, ExcPETScError(ierr));
193 
194  create_vector(n, local_sz);
195  }
196 
197  // finally clear the new vector if so
198  // desired
199  if (omit_zeroing_entries == false)
200  *this = 0;
201  }
202 
203 
204 
205  void
206  Vector::reinit(const Vector &v, const bool omit_zeroing_entries)
207  {
208  if (v.has_ghost_elements())
209  {
211  if (!omit_zeroing_entries)
212  {
213  const PetscErrorCode ierr = VecSet(vector, 0.0);
214  AssertThrow(ierr == 0, ExcPETScError(ierr));
215  }
216  }
217  else
219  v.size(),
220  v.locally_owned_size(),
221  omit_zeroing_entries);
222  }
223 
224 
225 
226  void
227  Vector::reinit(const IndexSet &local,
228  const IndexSet &ghost,
229  const MPI_Comm &comm)
230  {
231  const PetscErrorCode ierr = VecDestroy(&vector);
232  AssertThrow(ierr == 0, ExcPETScError(ierr));
233 
234  communicator = comm;
235 
237 
238  IndexSet ghost_set = ghost;
239  ghost_set.subtract_set(local);
240 
241  create_vector(local.size(), local.n_elements(), ghost_set);
242  }
243 
244  void
245  Vector::reinit(const IndexSet &local, const MPI_Comm &comm)
246  {
247  const PetscErrorCode ierr = VecDestroy(&vector);
248  AssertThrow(ierr == 0, ExcPETScError(ierr));
249 
250  communicator = comm;
251 
253  Assert(local.size() > 0, ExcMessage("can not create vector of size 0."));
254  create_vector(local.size(), local.n_elements());
255  }
256 
257  void
259  const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner)
260  {
261  this->reinit(partitioner->locally_owned_range(),
262  partitioner->ghost_indices(),
263  partitioner->get_mpi_communicator());
264  }
265 
266 
267  void
269  {
270  (void)n;
271  AssertIndexRange(locally_owned_size, n + 1);
272  ghosted = false;
273 
274  const PetscErrorCode ierr = VecCreateMPI(communicator,
275  locally_owned_size,
276  PETSC_DETERMINE,
277  &vector);
278  AssertThrow(ierr == 0, ExcPETScError(ierr));
279 
280  Assert(size() == n, ExcDimensionMismatch(size(), n));
281  }
282 
283 
284 
285  void
288  const IndexSet &ghostnodes)
289  {
290  (void)n;
291  AssertIndexRange(locally_owned_size, n + 1);
292  ghosted = true;
293  ghost_indices = ghostnodes;
294 
295  std::vector<size_type> ghostindices;
296  ghostnodes.fill_index_vector(ghostindices);
297 
298  const PetscInt *ptr =
299  (ghostindices.size() > 0 ?
300  reinterpret_cast<const PetscInt *>(ghostindices.data()) :
301  nullptr);
302 
303  PetscErrorCode ierr = VecCreateGhost(communicator,
304  locally_owned_size,
305  PETSC_DETERMINE,
306  ghostindices.size(),
307  ptr,
308  &vector);
309  AssertThrow(ierr == 0, ExcPETScError(ierr));
310 
311  Assert(size() == n, ExcDimensionMismatch(size(), n));
312 
313 # if DEBUG
314  {
315  // test ghost allocation in debug mode
316  PetscInt begin, end;
317 
318  ierr = VecGetOwnershipRange(vector, &begin, &end);
319  AssertThrow(ierr == 0, ExcPETScError(ierr));
320 
321  AssertDimension(locally_owned_size,
322  static_cast<size_type>(end - begin));
323 
324  Vec l;
325  ierr = VecGhostGetLocalForm(vector, &l);
326  AssertThrow(ierr == 0, ExcPETScError(ierr));
327 
328  PetscInt lsize;
329  ierr = VecGetSize(l, &lsize);
330  AssertThrow(ierr == 0, ExcPETScError(ierr));
331 
332  ierr = VecGhostRestoreLocalForm(vector, &l);
333  AssertThrow(ierr == 0, ExcPETScError(ierr));
334 
335  AssertDimension(lsize,
336  end - begin +
337  static_cast<PetscInt>(ghost_indices.n_elements()));
338  }
339 # endif
340 
341 
342  // in PETSc versions up to 3.5, VecCreateGhost zeroed out the locally
343  // owned vector elements but forgot about the ghost elements. we need to
344  // do this ourselves
345  //
346  // see https://code.google.com/p/dealii/issues/detail?id=233
347 # if DEAL_II_PETSC_VERSION_LT(3, 6, 0)
349  zero.reinit(communicator, this->size(), locally_owned_size);
350  *this = zero;
351 # endif
352  }
353 
354 
355 
356  bool
358  {
359  unsigned int has_nonzero = VectorBase::all_zero() ? 0 : 1;
360 # ifdef DEAL_II_WITH_MPI
361  // in parallel, check that the vector
362  // is zero on _all_ processors.
363  unsigned int num_nonzero = Utilities::MPI::sum(has_nonzero, communicator);
364  return num_nonzero == 0;
365 # else
366  return has_nonzero == 0;
367 # endif
368  }
369 
370 
371  void
372  Vector::print(std::ostream & out,
373  const unsigned int precision,
374  const bool scientific,
375  const bool across) const
376  {
377  AssertThrow(out, ExcIO());
378 
379  // get a representation of the vector and
380  // loop over all the elements
381  PetscScalar *val;
382  PetscInt nlocal, istart, iend;
383 
384  PetscErrorCode ierr = VecGetArray(vector, &val);
385  AssertThrow(ierr == 0, ExcPETScError(ierr));
386 
387  ierr = VecGetLocalSize(vector, &nlocal);
388  AssertThrow(ierr == 0, ExcPETScError(ierr));
389 
390  ierr = VecGetOwnershipRange(vector, &istart, &iend);
391  AssertThrow(ierr == 0, ExcPETScError(ierr));
392 
393  // save the state of out stream
394  std::ios::fmtflags old_flags = out.flags();
395  unsigned int old_precision = out.precision(precision);
396 
397  out.precision(precision);
398  if (scientific)
399  out.setf(std::ios::scientific, std::ios::floatfield);
400  else
401  out.setf(std::ios::fixed, std::ios::floatfield);
402 
403  // let each processor produce its output in turn. this requires
404  // synchronizing output between processors using a barrier --
405  // which is clearly slow, but nobody is going to print a whole
406  // matrix this way on a regular basis for production runs, so
407  // the slowness of the barrier doesn't matter
408  for (unsigned int i = 0;
410  i++)
411  {
412  const int mpi_ierr = MPI_Barrier(communicator);
413  AssertThrowMPI(mpi_ierr);
414 
416  {
417  if (across)
418  {
419  out << "[Proc" << i << " " << istart << "-" << iend - 1 << "]"
420  << ' ';
421  for (PetscInt i = 0; i < nlocal; ++i)
422  out << val[i] << ' ';
423  }
424  else
425  {
426  out << "[Proc " << i << " " << istart << "-" << iend - 1
427  << "]" << std::endl;
428  for (PetscInt i = 0; i < nlocal; ++i)
429  out << val[i] << std::endl;
430  }
431  out << std::endl;
432  }
433  }
434  // reset output format
435  out.flags(old_flags);
436  out.precision(old_precision);
437 
438  // restore the representation of the
439  // vector
440  ierr = VecRestoreArray(vector, &val);
441  AssertThrow(ierr == 0, ExcPETScError(ierr));
442 
443  AssertThrow(out, ExcIO());
444  }
445 
446  } // namespace MPI
447 
448 } // namespace PETScWrappers
449 
451 
452 #endif // DEAL_II_WITH_PETSC
virtual void create_vector(const size_type n, const size_type locally_owned_size)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1655
static ::ExceptionBase & ExcIO()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1720
size_type locally_owned_size() const
#define AssertThrow(cond, exc)
Definition: exceptions.h:1571
virtual void clear() override
bool is_ascending_and_one_to_one(const MPI_Comm &communicator) const
Definition: index_set.cc:674
size_type size() const
Definition: index_set.h:1639
static ::ExceptionBase & ExcMessage(std::string arg1)
void reinit(const MPI_Comm &communicator, const size_type N, const size_type locally_owned_size, const bool omit_zeroing_entries=false)
T sum(const T &t, const MPI_Comm &mpi_communicator)
void subtract_set(const IndexSet &other)
Definition: index_set.cc:266
VectorOperation::values last_action
#define Assert(cond, exc)
Definition: exceptions.h:1461
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
VectorType::value_type * end(VectorType &V)
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:117
void fill_index_vector(std::vector< size_type > &indices) const
Definition: index_set.cc:515
Vector & operator=(const Vector &v)
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1778
bool has_ghost_elements() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
VectorType::value_type * begin(VectorType &V)
const MPI_Comm & comm
IndexSet locally_owned_elements() const
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:128
static ::ExceptionBase & ExcNotImplemented()
static const types::blas_int zero
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:1837
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)