deal.II version GIT relicensing-2017-g7ae82fad8b 2024-10-22 19:20:00+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 // If the size of the index set can be converted to a PetscInt then every
297 // index can also be converted
298 AssertThrowIntegerConversion(static_cast<PetscInt>(n), n);
299 ghosted = true;
300 ghost_indices = ghostnodes;
301
302 std::size_t i = 0;
303 std::vector<PetscInt> petsc_ghost_indices(ghostnodes.n_elements());
304 for (const auto &index : ghostnodes)
305 {
306 petsc_ghost_indices[i] = static_cast<PetscInt>(index);
307 ++i;
308 }
309
310 PetscErrorCode ierr = VecCreateGhost(communicator,
312 PETSC_DETERMINE,
313 petsc_ghost_indices.size(),
314 petsc_ghost_indices.data(),
315 &vector);
316 AssertThrow(ierr == 0, ExcPETScError(ierr));
317
318 Assert(size() == n, ExcDimensionMismatch(size(), n));
319
320# ifdef DEBUG
321 {
322 // test ghost allocation in debug mode
323 PetscInt begin, end;
324
325 ierr = VecGetOwnershipRange(vector, &begin, &end);
326 AssertThrow(ierr == 0, ExcPETScError(ierr));
327
329 static_cast<size_type>(end - begin));
330
331 Vec l;
332 ierr = VecGhostGetLocalForm(vector, &l);
333 AssertThrow(ierr == 0, ExcPETScError(ierr));
334
335 PetscInt lsize;
336 ierr = VecGetSize(l, &lsize);
337 AssertThrow(ierr == 0, ExcPETScError(ierr));
338
339 ierr = VecGhostRestoreLocalForm(vector, &l);
340 AssertThrow(ierr == 0, ExcPETScError(ierr));
341
342 AssertDimension(lsize,
343 end - begin +
344 static_cast<PetscInt>(ghost_indices.n_elements()));
345 }
346# endif
347 }
348
349
350
351 bool
353 {
354 unsigned int has_nonzero = VectorBase::all_zero() ? 0 : 1;
355 // in parallel, check that the vector
356 // is zero on _all_ processors.
357 unsigned int num_nonzero =
358 Utilities::MPI::sum(has_nonzero, this->get_mpi_communicator());
359 return num_nonzero == 0;
360 }
361
362
363 void
364 Vector::print(std::ostream &out,
365 const unsigned int precision,
366 const bool scientific,
367 const bool across) const
368 {
369 AssertThrow(out.fail() == false, ExcIO());
370
371 // get a representation of the vector and
372 // loop over all the elements
373 const PetscScalar *val;
374 PetscInt nlocal, istart, iend;
375
376 PetscErrorCode ierr = VecGetArrayRead(vector, &val);
377 AssertThrow(ierr == 0, ExcPETScError(ierr));
378
379 ierr = VecGetLocalSize(vector, &nlocal);
380 AssertThrow(ierr == 0, ExcPETScError(ierr));
381
382 ierr = VecGetOwnershipRange(vector, &istart, &iend);
383 AssertThrow(ierr == 0, ExcPETScError(ierr));
384
385 // save the state of out stream
386 std::ios::fmtflags old_flags = out.flags();
387 unsigned int old_precision = out.precision(precision);
388
389 out.precision(precision);
390 if (scientific)
391 out.setf(std::ios::scientific, std::ios::floatfield);
392 else
393 out.setf(std::ios::fixed, std::ios::floatfield);
394
395 // let each processor produce its output in turn. this requires
396 // synchronizing output between processors using a barrier --
397 // which is clearly slow, but nobody is going to print a whole
398 // matrix this way on a regular basis for production runs, so
399 // the slowness of the barrier doesn't matter
400 MPI_Comm communicator = this->get_mpi_communicator();
401 for (unsigned int i = 0;
402 i < Utilities::MPI::n_mpi_processes(communicator);
403 i++)
404 {
405 const int mpi_ierr = MPI_Barrier(communicator);
406 AssertThrowMPI(mpi_ierr);
407
408 if (i == Utilities::MPI::this_mpi_process(communicator))
409 {
410 if (across)
411 {
412 out << "[Proc" << i << " " << istart << "-" << iend - 1 << "]"
413 << ' ';
414 for (PetscInt i = 0; i < nlocal; ++i)
415 out << val[i] << ' ';
416 }
417 else
418 {
419 out << "[Proc " << i << " " << istart << "-" << iend - 1
420 << "]" << std::endl;
421 for (PetscInt i = 0; i < nlocal; ++i)
422 out << val[i] << std::endl;
423 }
424 out << std::endl;
425 }
426 }
427 // reset output format
428 out.flags(old_flags);
429 out.precision(old_precision);
430
431 // restore the representation of the
432 // vector
433 ierr = VecRestoreArrayRead(vector, &val);
434 AssertThrow(ierr == 0, ExcPETScError(ierr));
435
436 AssertThrow(out.fail() == false, ExcIO());
437 }
438
439 } // namespace MPI
440
441} // namespace PETScWrappers
442
444
445#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:1776
size_type n_elements() const
Definition index_set.h:1934
void subtract_set(const IndexSet &other)
Definition index_set.cc:473
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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
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)
#define AssertThrowIntegerConversion(index1, index2)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:92
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:107
*braid_SplitCommworld & comm