Reference documentation for deal.II version 9.3.0
\(\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 
258  void
260  {
261  (void)n;
262  AssertIndexRange(locally_owned_size, n + 1);
263  ghosted = false;
264 
265  const PetscErrorCode ierr = VecCreateMPI(communicator,
266  locally_owned_size,
267  PETSC_DETERMINE,
268  &vector);
269  AssertThrow(ierr == 0, ExcPETScError(ierr));
270 
271  Assert(size() == n, ExcDimensionMismatch(size(), n));
272  }
273 
274 
275 
276  void
279  const IndexSet &ghostnodes)
280  {
281  (void)n;
282  AssertIndexRange(locally_owned_size, n + 1);
283  ghosted = true;
284  ghost_indices = ghostnodes;
285 
286  std::vector<size_type> ghostindices;
287  ghostnodes.fill_index_vector(ghostindices);
288 
289  const PetscInt *ptr =
290  (ghostindices.size() > 0 ?
291  reinterpret_cast<const PetscInt *>(ghostindices.data()) :
292  nullptr);
293 
294  PetscErrorCode ierr = VecCreateGhost(communicator,
295  locally_owned_size,
296  PETSC_DETERMINE,
297  ghostindices.size(),
298  ptr,
299  &vector);
300  AssertThrow(ierr == 0, ExcPETScError(ierr));
301 
302  Assert(size() == n, ExcDimensionMismatch(size(), n));
303 
304 # if DEBUG
305  {
306  // test ghost allocation in debug mode
307  PetscInt begin, end;
308 
309  ierr = VecGetOwnershipRange(vector, &begin, &end);
310  AssertThrow(ierr == 0, ExcPETScError(ierr));
311 
312  AssertDimension(locally_owned_size,
313  static_cast<size_type>(end - begin));
314 
315  Vec l;
316  ierr = VecGhostGetLocalForm(vector, &l);
317  AssertThrow(ierr == 0, ExcPETScError(ierr));
318 
319  PetscInt lsize;
320  ierr = VecGetSize(l, &lsize);
321  AssertThrow(ierr == 0, ExcPETScError(ierr));
322 
323  ierr = VecGhostRestoreLocalForm(vector, &l);
324  AssertThrow(ierr == 0, ExcPETScError(ierr));
325 
326  AssertDimension(lsize,
327  end - begin +
328  static_cast<PetscInt>(ghost_indices.n_elements()));
329  }
330 # endif
331 
332 
333  // in PETSc versions up to 3.5, VecCreateGhost zeroed out the locally
334  // owned vector elements but forgot about the ghost elements. we need to
335  // do this ourselves
336  //
337  // see https://code.google.com/p/dealii/issues/detail?id=233
338 # if DEAL_II_PETSC_VERSION_LT(3, 6, 0)
340  zero.reinit(communicator, this->size(), locally_owned_size);
341  *this = zero;
342 # endif
343  }
344 
345 
346 
347  bool
349  {
350  unsigned int has_nonzero = VectorBase::all_zero() ? 0 : 1;
351 # ifdef DEAL_II_WITH_MPI
352  // in parallel, check that the vector
353  // is zero on _all_ processors.
354  unsigned int num_nonzero = Utilities::MPI::sum(has_nonzero, communicator);
355  return num_nonzero == 0;
356 # else
357  return has_nonzero == 0;
358 # endif
359  }
360 
361 
362  void
363  Vector::print(std::ostream & out,
364  const unsigned int precision,
365  const bool scientific,
366  const bool across) const
367  {
368  AssertThrow(out, ExcIO());
369 
370  // get a representation of the vector and
371  // loop over all the elements
372  PetscScalar *val;
373  PetscInt nlocal, istart, iend;
374 
375  PetscErrorCode ierr = VecGetArray(vector, &val);
376  AssertThrow(ierr == 0, ExcPETScError(ierr));
377 
378  ierr = VecGetLocalSize(vector, &nlocal);
379  AssertThrow(ierr == 0, ExcPETScError(ierr));
380 
381  ierr = VecGetOwnershipRange(vector, &istart, &iend);
382  AssertThrow(ierr == 0, ExcPETScError(ierr));
383 
384  // save the state of out stream
385  std::ios::fmtflags old_flags = out.flags();
386  unsigned int old_precision = out.precision(precision);
387 
388  out.precision(precision);
389  if (scientific)
390  out.setf(std::ios::scientific, std::ios::floatfield);
391  else
392  out.setf(std::ios::fixed, std::ios::floatfield);
393 
394  // let each processor produce its output in turn. this requires
395  // synchronizing output between processors using a barrier --
396  // which is clearly slow, but nobody is going to print a whole
397  // matrix this way on a regular basis for production runs, so
398  // the slowness of the barrier doesn't matter
399  for (unsigned int i = 0;
401  i++)
402  {
403  const int mpi_ierr = MPI_Barrier(communicator);
404  AssertThrowMPI(mpi_ierr);
405 
407  {
408  if (across)
409  {
410  out << "[Proc" << i << " " << istart << "-" << iend - 1 << "]"
411  << ' ';
412  for (PetscInt i = 0; i < nlocal; ++i)
413  out << val[i] << ' ';
414  }
415  else
416  {
417  out << "[Proc " << i << " " << istart << "-" << iend - 1
418  << "]" << std::endl;
419  for (PetscInt i = 0; i < nlocal; ++i)
420  out << val[i] << std::endl;
421  }
422  out << std::endl;
423  }
424  }
425  // reset output format
426  out.flags(old_flags);
427  out.precision(old_precision);
428 
429  // restore the representation of the
430  // vector
431  ierr = VecRestoreArray(vector, &val);
432  AssertThrow(ierr == 0, ExcPETScError(ierr));
433 
434  AssertThrow(out, ExcIO());
435  }
436 
437  } // namespace MPI
438 
439 } // namespace PETScWrappers
440 
442 
443 #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:1622
static ::ExceptionBase & ExcIO()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
size_type locally_owned_size() const
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
virtual void clear() override
bool is_ascending_and_one_to_one(const MPI_Comm &communicator) const
Definition: index_set.cc:666
size_type size() const
Definition: index_set.h:1634
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:258
VectorOperation::values last_action
#define Assert(cond, exc)
Definition: exceptions.h:1465
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
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:507
Vector & operator=(const Vector &v)
*braid_SplitCommworld & comm
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1746
bool has_ghost_elements() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
VectorType::value_type * begin(VectorType &V)
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:1832
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)