Reference documentation for deal.II version GIT relicensing-214-g6e74dec06b 2024-03-27 18:10: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\}}\)
Loading...
Searching...
No Matches
petsc_parallel_vector.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2004 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#include <deal.II/base/mpi.h>
16
18
19#ifdef DEAL_II_WITH_PETSC
20
21# include <algorithm>
22# include <cmath>
23
25
26namespace PETScWrappers
27{
28 namespace MPI
29 {
30 // Check that the class we declare here satisfies the
31 // vector-space-vector concept. If we catch it here,
32 // any mistake in the vector class declaration would
33 // show up in uses of this class later on as well.
34# ifdef DEAL_II_HAVE_CXX20
36# endif
37
39 {
40 // virtual functions called in constructors and destructors never use the
41 // override in a derived class
42 // for clarity be explicit on which function is called
43 Vector::create_vector(MPI_COMM_SELF, 0, 0);
44 }
45
46
47
48 Vector::Vector(const MPI_Comm communicator,
49 const size_type n,
50 const size_type locally_owned_size)
51 {
53 }
54
55
56
58 const IndexSet &ghost,
59 const MPI_Comm communicator)
60 {
61 Assert(local.is_ascending_and_one_to_one(communicator),
63
64 IndexSet ghost_set = ghost;
65 ghost_set.subtract_set(local);
66
67 Vector::create_vector(communicator,
68 local.size(),
69 local.n_elements(),
70 ghost_set);
71 }
72
73
74
76 : VectorBase()
77 {
78 if (v.has_ghost_elements())
80 v.size(),
83 else
85 v.size(),
87
88 this->operator=(v);
89 }
90
91
92
93 Vector::Vector(const IndexSet &local, const MPI_Comm communicator)
94 {
95 Assert(local.is_ascending_and_one_to_one(communicator),
97 Vector::create_vector(communicator, local.size(), local.n_elements());
98 }
99
100
101
102 Vector &
104 {
105 // make sure left- and right-hand side of the assignment are
106 // compress()'ed:
108 internal::VectorReference::ExcWrongMode(VectorOperation::unknown,
109 v.last_action));
111 internal::VectorReference::ExcWrongMode(VectorOperation::unknown,
112 last_action));
113
114 // if the vectors have different sizes,
115 // then first resize the present one
116 if (size() != v.size())
117 {
118 if (v.has_ghost_elements())
122 else
124 v.size(),
126 true);
127 }
128
129 PetscErrorCode ierr = VecCopy(v.vector, vector);
130 AssertThrow(ierr == 0, ExcPETScError(ierr));
131
132 if (has_ghost_elements())
133 {
134 ierr = VecGhostUpdateBegin(vector, INSERT_VALUES, SCATTER_FORWARD);
135 AssertThrow(ierr == 0, ExcPETScError(ierr));
136 ierr = VecGhostUpdateEnd(vector, INSERT_VALUES, SCATTER_FORWARD);
137 AssertThrow(ierr == 0, ExcPETScError(ierr));
138 }
139 return *this;
140 }
141
142
143
144 void
146 {
148
149 create_vector(MPI_COMM_SELF, 0, 0);
150 }
151
152
153
154 void
155 Vector::reinit(const MPI_Comm communicator,
156 const size_type n,
157 const size_type local_sz,
158 const bool omit_zeroing_entries)
159 {
160 // only do something if the sizes
161 // mismatch (may not be true for every proc)
162
163 int k_global, k = ((size() != n) || (locally_owned_size() != local_sz));
164 {
165 const int ierr =
166 MPI_Allreduce(&k, &k_global, 1, MPI_INT, MPI_LOR, communicator);
167 AssertThrowMPI(ierr);
168 }
169
170 if (k_global || has_ghost_elements())
171 {
172 // FIXME: I'd like to use this here,
173 // but somehow it leads to odd errors
174 // somewhere down the line in some of
175 // the tests:
176 // const PetscErrorCode ierr = VecSetSizes (vector, n, n);
177 // AssertThrow (ierr == 0, ExcPETScError(ierr));
178
179 // so let's go the slow way:
180
181 const PetscErrorCode ierr = VecDestroy(&vector);
182 AssertThrow(ierr == 0, ExcPETScError(ierr));
183
184 create_vector(communicator, n, local_sz);
185 }
186
187 // finally clear the new vector if so
188 // desired
189 if (omit_zeroing_entries == false)
190 *this = 0;
191 }
192
193
194
195 void
196 Vector::reinit(const Vector &v, const bool omit_zeroing_entries)
197 {
198 if (v.has_ghost_elements())
199 {
203 if (!omit_zeroing_entries)
204 {
205 const PetscErrorCode ierr = VecSet(vector, 0.0);
206 AssertThrow(ierr == 0, ExcPETScError(ierr));
207 }
208 }
209 else
211 v.size(),
213 omit_zeroing_entries);
214 }
215
216
217
218 void
220 const IndexSet &ghost,
221 const MPI_Comm comm)
222 {
223 const PetscErrorCode ierr = VecDestroy(&vector);
224 AssertThrow(ierr == 0, ExcPETScError(ierr));
225
227
228 IndexSet ghost_set = ghost;
229 ghost_set.subtract_set(local);
230
231 create_vector(comm, local.size(), local.n_elements(), ghost_set);
232 }
233
234 void
236 {
237 const PetscErrorCode ierr = VecDestroy(&vector);
238 AssertThrow(ierr == 0, ExcPETScError(ierr));
239
241 Assert(local.size() > 0, ExcMessage("can not create vector of size 0."));
242 create_vector(comm, local.size(), local.n_elements());
243 }
244
245 void
247 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
248 const bool make_ghosted)
249 {
250 if (make_ghosted)
251 {
252 Assert(partitioner->ghost_indices_initialized(),
253 ExcMessage("You asked to create a ghosted vector, but the "
254 "partitioner does not provide ghost indices."));
255
256 this->reinit(partitioner->locally_owned_range(),
257 partitioner->ghost_indices(),
258 partitioner->get_mpi_communicator());
259 }
260 else
261 {
262 this->reinit(partitioner->locally_owned_range(),
263 partitioner->get_mpi_communicator());
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 {
273 (void)n;
275 ghosted = false;
276
277 const PetscErrorCode ierr = VecCreateMPI(communicator,
279 PETSC_DETERMINE,
280 &vector);
281 AssertThrow(ierr == 0, ExcPETScError(ierr));
282
283 Assert(size() == n, ExcDimensionMismatch(size(), n));
284 }
285
286
287
288 void
289 Vector::create_vector(const MPI_Comm communicator,
290 const size_type n,
291 const size_type locally_owned_size,
292 const IndexSet &ghostnodes)
293 {
294 (void)n;
296 ghosted = true;
297 ghost_indices = ghostnodes;
298
299 const std::vector<size_type> ghostindices = ghostnodes.get_index_vector();
300
301 const PetscInt *ptr =
302 (ghostindices.size() > 0 ?
303 reinterpret_cast<const PetscInt *>(ghostindices.data()) :
304 nullptr);
305
306 PetscErrorCode ierr = VecCreateGhost(communicator,
308 PETSC_DETERMINE,
309 ghostindices.size(),
310 ptr,
311 &vector);
312 AssertThrow(ierr == 0, ExcPETScError(ierr));
313
314 Assert(size() == n, ExcDimensionMismatch(size(), n));
315
316# ifdef DEBUG
317 {
318 // test ghost allocation in debug mode
319 PetscInt begin, end;
320
321 ierr = VecGetOwnershipRange(vector, &begin, &end);
322 AssertThrow(ierr == 0, ExcPETScError(ierr));
323
325 static_cast<size_type>(end - begin));
326
327 Vec l;
328 ierr = VecGhostGetLocalForm(vector, &l);
329 AssertThrow(ierr == 0, ExcPETScError(ierr));
330
331 PetscInt lsize;
332 ierr = VecGetSize(l, &lsize);
333 AssertThrow(ierr == 0, ExcPETScError(ierr));
334
335 ierr = VecGhostRestoreLocalForm(vector, &l);
336 AssertThrow(ierr == 0, ExcPETScError(ierr));
337
338 AssertDimension(lsize,
339 end - begin +
340 static_cast<PetscInt>(ghost_indices.n_elements()));
341 }
342# endif
343 }
344
345
346
347 bool
349 {
350 unsigned int has_nonzero = VectorBase::all_zero() ? 0 : 1;
351 // in parallel, check that the vector
352 // is zero on _all_ processors.
353 unsigned int num_nonzero =
354 Utilities::MPI::sum(has_nonzero, this->get_mpi_communicator());
355 return num_nonzero == 0;
356 }
357
358
359 void
360 Vector::print(std::ostream &out,
361 const unsigned int precision,
362 const bool scientific,
363 const bool across) const
364 {
365 AssertThrow(out.fail() == false, ExcIO());
366
367 // get a representation of the vector and
368 // loop over all the elements
369 const PetscScalar *val;
370 PetscInt nlocal, istart, iend;
371
372 PetscErrorCode ierr = VecGetArrayRead(vector, &val);
373 AssertThrow(ierr == 0, ExcPETScError(ierr));
374
375 ierr = VecGetLocalSize(vector, &nlocal);
376 AssertThrow(ierr == 0, ExcPETScError(ierr));
377
378 ierr = VecGetOwnershipRange(vector, &istart, &iend);
379 AssertThrow(ierr == 0, ExcPETScError(ierr));
380
381 // save the state of out stream
382 std::ios::fmtflags old_flags = out.flags();
383 unsigned int old_precision = out.precision(precision);
384
385 out.precision(precision);
386 if (scientific)
387 out.setf(std::ios::scientific, std::ios::floatfield);
388 else
389 out.setf(std::ios::fixed, std::ios::floatfield);
390
391 // let each processor produce its output in turn. this requires
392 // synchronizing output between processors using a barrier --
393 // which is clearly slow, but nobody is going to print a whole
394 // matrix this way on a regular basis for production runs, so
395 // the slowness of the barrier doesn't matter
396 MPI_Comm communicator = this->get_mpi_communicator();
397 for (unsigned int i = 0;
398 i < Utilities::MPI::n_mpi_processes(communicator);
399 i++)
400 {
401 const int mpi_ierr = MPI_Barrier(communicator);
402 AssertThrowMPI(mpi_ierr);
403
404 if (i == Utilities::MPI::this_mpi_process(communicator))
405 {
406 if (across)
407 {
408 out << "[Proc" << i << " " << istart << "-" << iend - 1 << "]"
409 << ' ';
410 for (PetscInt i = 0; i < nlocal; ++i)
411 out << val[i] << ' ';
412 }
413 else
414 {
415 out << "[Proc " << i << " " << istart << "-" << iend - 1
416 << "]" << std::endl;
417 for (PetscInt i = 0; i < nlocal; ++i)
418 out << val[i] << std::endl;
419 }
420 out << std::endl;
421 }
422 }
423 // reset output format
424 out.flags(old_flags);
425 out.precision(old_precision);
426
427 // restore the representation of the
428 // vector
429 ierr = VecRestoreArrayRead(vector, &val);
430 AssertThrow(ierr == 0, ExcPETScError(ierr));
431
432 AssertThrow(out.fail() == false, ExcIO());
433 }
434
435 } // namespace MPI
436
437} // namespace PETScWrappers
438
440
441#endif // DEAL_II_WITH_PETSC
bool is_ascending_and_one_to_one(const MPI_Comm communicator) const
size_type size() const
Definition index_set.h:1765
size_type n_elements() const
Definition index_set.h:1923
void subtract_set(const IndexSet &other)
Definition index_set.cc:470
std::vector< size_type > get_index_vector() const
Definition index_set.cc:923
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
Vector & operator=(const Vector &v)
virtual void create_vector(const MPI_Comm comm, const size_type n, const size_type locally_owned_size)
void reinit(const MPI_Comm communicator, const size_type N, const size_type locally_owned_size, const bool omit_zeroing_entries=false)
VectorOperation::values last_action
IndexSet locally_owned_elements() const
MPI_Comm get_mpi_communicator() const
bool has_ghost_elements() const
size_type locally_owned_size() const
size_type size() const override
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:128
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:143
*braid_SplitCommworld & comm