Reference documentation for deal.II version GIT 8e09676776 2023-03-27 21:15:01+00:00
\(\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 - 2022 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  {
33  // virtual functions called in constructors and destructors never use the
34  // override in a derived class
35  // for clarity be explicit on which function is called
36  Vector::create_vector(MPI_COMM_SELF, 0, 0);
37  }
38 
39 
40 
41  Vector::Vector(const MPI_Comm &communicator,
42  const size_type n,
43  const size_type locally_owned_size)
44  {
45  Vector::create_vector(communicator, n, locally_owned_size);
46  }
47 
48 
49 
50  Vector::Vector(const IndexSet &local,
51  const IndexSet &ghost,
52  const MPI_Comm &communicator)
53  {
54  Assert(local.is_ascending_and_one_to_one(communicator),
56 
57  IndexSet ghost_set = ghost;
58  ghost_set.subtract_set(local);
59 
60  Vector::create_vector(communicator,
61  local.size(),
62  local.n_elements(),
63  ghost_set);
64  }
65 
66 
67 
69  : VectorBase()
70  {
71  if (v.has_ghost_elements())
73  v.size(),
75  v.ghost_indices);
76  else
78  v.size(),
79  v.locally_owned_size());
80 
81  this->operator=(v);
82  }
83 
84 
85 
86  Vector::Vector(const IndexSet &local, const MPI_Comm &communicator)
87  {
88  Assert(local.is_ascending_and_one_to_one(communicator),
90  Vector::create_vector(communicator, local.size(), local.n_elements());
91  }
92 
93 
94 
95  Vector &
97  {
98  // make sure left- and right-hand side of the assignment are
99  // compress()'ed:
101  internal::VectorReference::ExcWrongMode(VectorOperation::unknown,
102  v.last_action));
104  internal::VectorReference::ExcWrongMode(VectorOperation::unknown,
105  last_action));
106 
107  // if the vectors have different sizes,
108  // then first resize the present one
109  if (size() != v.size())
110  {
111  if (v.has_ghost_elements())
113  v.ghost_indices,
115  else
117  v.size(),
118  v.locally_owned_size(),
119  true);
120  }
121 
122  PetscErrorCode ierr = VecCopy(v.vector, vector);
123  AssertThrow(ierr == 0, ExcPETScError(ierr));
124 
125  if (has_ghost_elements())
126  {
127  ierr = VecGhostUpdateBegin(vector, INSERT_VALUES, SCATTER_FORWARD);
128  AssertThrow(ierr == 0, ExcPETScError(ierr));
129  ierr = VecGhostUpdateEnd(vector, INSERT_VALUES, SCATTER_FORWARD);
130  AssertThrow(ierr == 0, ExcPETScError(ierr));
131  }
132  return *this;
133  }
134 
135 
136 
137  void
139  {
141 
142  create_vector(MPI_COMM_SELF, 0, 0);
143  }
144 
145 
146 
147  void
148  Vector::reinit(const MPI_Comm &communicator,
149  const size_type n,
150  const size_type local_sz,
151  const bool omit_zeroing_entries)
152  {
153  // only do something if the sizes
154  // mismatch (may not be true for every proc)
155 
156  int k_global, k = ((size() != n) || (locally_owned_size() != local_sz));
157  {
158  const int ierr =
159  MPI_Allreduce(&k, &k_global, 1, MPI_INT, MPI_LOR, communicator);
160  AssertThrowMPI(ierr);
161  }
162 
163  if (k_global || has_ghost_elements())
164  {
165  // FIXME: I'd like to use this here,
166  // but somehow it leads to odd errors
167  // somewhere down the line in some of
168  // the tests:
169  // const PetscErrorCode ierr = VecSetSizes (vector, n, n);
170  // AssertThrow (ierr == 0, ExcPETScError(ierr));
171 
172  // so let's go the slow way:
173 
174  const PetscErrorCode ierr = VecDestroy(&vector);
175  AssertThrow(ierr == 0, ExcPETScError(ierr));
176 
177  create_vector(communicator, n, local_sz);
178  }
179 
180  // finally clear the new vector if so
181  // desired
182  if (omit_zeroing_entries == false)
183  *this = 0;
184  }
185 
186 
187 
188  void
189  Vector::reinit(const Vector &v, const bool omit_zeroing_entries)
190  {
191  if (v.has_ghost_elements())
192  {
194  v.ghost_indices,
196  if (!omit_zeroing_entries)
197  {
198  const PetscErrorCode ierr = VecSet(vector, 0.0);
199  AssertThrow(ierr == 0, ExcPETScError(ierr));
200  }
201  }
202  else
204  v.size(),
205  v.locally_owned_size(),
206  omit_zeroing_entries);
207  }
208 
209 
210 
211  void
212  Vector::reinit(const IndexSet &local,
213  const IndexSet &ghost,
214  const MPI_Comm &comm)
215  {
216  const PetscErrorCode ierr = VecDestroy(&vector);
217  AssertThrow(ierr == 0, ExcPETScError(ierr));
218 
220 
221  IndexSet ghost_set = ghost;
222  ghost_set.subtract_set(local);
223 
224  create_vector(comm, local.size(), local.n_elements(), ghost_set);
225  }
226 
227  void
228  Vector::reinit(const IndexSet &local, const MPI_Comm &comm)
229  {
230  const PetscErrorCode ierr = VecDestroy(&vector);
231  AssertThrow(ierr == 0, ExcPETScError(ierr));
232 
234  Assert(local.size() > 0, ExcMessage("can not create vector of size 0."));
235  create_vector(comm, local.size(), local.n_elements());
236  }
237 
238  void
240  const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner)
241  {
242  this->reinit(partitioner->locally_owned_range(),
243  partitioner->ghost_indices(),
244  partitioner->get_mpi_communicator());
245  }
246 
247 
248  void
249  Vector::create_vector(const MPI_Comm &communicator,
250  const size_type n,
251  const size_type locally_owned_size)
252  {
253  (void)n;
255  ghosted = false;
256 
257  const PetscErrorCode ierr = VecCreateMPI(communicator,
259  PETSC_DETERMINE,
260  &vector);
261  AssertThrow(ierr == 0, ExcPETScError(ierr));
262 
263  Assert(size() == n, ExcDimensionMismatch(size(), n));
264  }
265 
266 
267 
268  void
269  Vector::create_vector(const MPI_Comm &communicator,
270  const size_type n,
271  const size_type locally_owned_size,
272  const IndexSet &ghostnodes)
273  {
274  (void)n;
276  ghosted = true;
277  ghost_indices = ghostnodes;
278 
279  const std::vector<size_type> ghostindices = ghostnodes.get_index_vector();
280 
281  const PetscInt *ptr =
282  (ghostindices.size() > 0 ?
283  reinterpret_cast<const PetscInt *>(ghostindices.data()) :
284  nullptr);
285 
286  PetscErrorCode ierr = VecCreateGhost(communicator,
288  PETSC_DETERMINE,
289  ghostindices.size(),
290  ptr,
291  &vector);
292  AssertThrow(ierr == 0, ExcPETScError(ierr));
293 
294  Assert(size() == n, ExcDimensionMismatch(size(), n));
295 
296 # ifdef DEBUG
297  {
298  // test ghost allocation in debug mode
299  PetscInt begin, end;
300 
301  ierr = VecGetOwnershipRange(vector, &begin, &end);
302  AssertThrow(ierr == 0, ExcPETScError(ierr));
303 
305  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 
326 
327  bool
329  {
330  unsigned int has_nonzero = VectorBase::all_zero() ? 0 : 1;
331 # ifdef DEAL_II_WITH_MPI
332  // in parallel, check that the vector
333  // is zero on _all_ processors.
334  unsigned int num_nonzero =
335  Utilities::MPI::sum(has_nonzero, this->get_mpi_communicator());
336  return num_nonzero == 0;
337 # else
338  return has_nonzero == 0;
339 # endif
340  }
341 
342 
343  void
344  Vector::print(std::ostream & out,
345  const unsigned int precision,
346  const bool scientific,
347  const bool across) const
348  {
349  AssertThrow(out.fail() == false, ExcIO());
350 
351  // get a representation of the vector and
352  // loop over all the elements
353  const PetscScalar *val;
354  PetscInt nlocal, istart, iend;
355 
356  PetscErrorCode ierr = VecGetArrayRead(vector, &val);
357  AssertThrow(ierr == 0, ExcPETScError(ierr));
358 
359  ierr = VecGetLocalSize(vector, &nlocal);
360  AssertThrow(ierr == 0, ExcPETScError(ierr));
361 
362  ierr = VecGetOwnershipRange(vector, &istart, &iend);
363  AssertThrow(ierr == 0, ExcPETScError(ierr));
364 
365  // save the state of out stream
366  std::ios::fmtflags old_flags = out.flags();
367  unsigned int old_precision = out.precision(precision);
368 
369  out.precision(precision);
370  if (scientific)
371  out.setf(std::ios::scientific, std::ios::floatfield);
372  else
373  out.setf(std::ios::fixed, std::ios::floatfield);
374 
375  // let each processor produce its output in turn. this requires
376  // synchronizing output between processors using a barrier --
377  // which is clearly slow, but nobody is going to print a whole
378  // matrix this way on a regular basis for production runs, so
379  // the slowness of the barrier doesn't matter
380  MPI_Comm communicator = this->get_mpi_communicator();
381  for (unsigned int i = 0;
382  i < Utilities::MPI::n_mpi_processes(communicator);
383  i++)
384  {
385  const int mpi_ierr = MPI_Barrier(communicator);
386  AssertThrowMPI(mpi_ierr);
387 
388  if (i == Utilities::MPI::this_mpi_process(communicator))
389  {
390  if (across)
391  {
392  out << "[Proc" << i << " " << istart << "-" << iend - 1 << "]"
393  << ' ';
394  for (PetscInt i = 0; i < nlocal; ++i)
395  out << val[i] << ' ';
396  }
397  else
398  {
399  out << "[Proc " << i << " " << istart << "-" << iend - 1
400  << "]" << std::endl;
401  for (PetscInt i = 0; i < nlocal; ++i)
402  out << val[i] << std::endl;
403  }
404  out << std::endl;
405  }
406  }
407  // reset output format
408  out.flags(old_flags);
409  out.precision(old_precision);
410 
411  // restore the representation of the
412  // vector
413  ierr = VecRestoreArrayRead(vector, &val);
414  AssertThrow(ierr == 0, ExcPETScError(ierr));
415 
416  AssertThrow(out.fail() == false, ExcIO());
417  }
418 
419  } // namespace MPI
420 
421 } // namespace PETScWrappers
422 
424 
425 #endif // DEAL_II_WITH_PETSC
size_type size() const
Definition: index_set.h:1648
size_type n_elements() const
Definition: index_set.h:1796
void subtract_set(const IndexSet &other)
Definition: index_set.cc:268
bool is_ascending_and_one_to_one(const MPI_Comm &communicator) const
Definition: index_set.cc:878
std::vector< size_type > get_index_vector() const
Definition: index_set.cc:690
types::global_dof_index size_type
Definition: petsc_vector.h:164
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
virtual void clear() override
Vector & operator=(const Vector &v)
virtual void create_vector(const MPI_Comm &comm, const size_type n, const size_type locally_owned_size)
VectorOperation::values last_action
IndexSet locally_owned_elements() const
const MPI_Comm & get_mpi_communicator() const
bool has_ghost_elements() const
size_type locally_owned_size() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1586
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1759
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1886
#define AssertIndexRange(index, range)
Definition: exceptions.h:1827
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1675
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:161
T sum(const T &t, const MPI_Comm &mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:150
const MPI_Comm & comm