deal.II version GIT relicensing-2791-g2e2a880805 2025-03-11 02:00: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_vector_base.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
16
17#ifdef DEAL_II_WITH_PETSC
18
20
24
25# include <cmath>
26
28
29namespace PETScWrappers
30{
31 namespace internal
32 {
33# ifndef DOXYGEN
34 VectorReference::operator PetscScalar() const
35 {
36 AssertIndexRange(index, vector.size());
37
38 // The vector may have ghost entries. In that case, we first need to
39 // figure out which elements we own locally, then get a pointer to the
40 // elements that are stored here (both the ones we own as well as the
41 // ghost elements). In this array, the locally owned elements come first
42 // followed by the ghost elements whose position we can get from an
43 // index set.
44 if (vector.ghosted)
45 {
46 PetscInt begin, end;
47 PetscErrorCode ierr =
48 VecGetOwnershipRange(vector.vector, &begin, &end);
49 AssertThrow(ierr == 0, ExcPETScError(ierr));
50
51 Vec locally_stored_elements = nullptr;
52 ierr = VecGhostGetLocalForm(vector.vector, &locally_stored_elements);
53 AssertThrow(ierr == 0, ExcPETScError(ierr));
54
55 PetscInt lsize;
56 ierr = VecGetSize(locally_stored_elements, &lsize);
57 AssertThrow(ierr == 0, ExcPETScError(ierr));
58
59 const PetscScalar *ptr;
60 ierr = VecGetArrayRead(locally_stored_elements, &ptr);
61 AssertThrow(ierr == 0, ExcPETScError(ierr));
62
63 PetscScalar value;
64
65 if (index >= static_cast<size_type>(begin) &&
66 index < static_cast<size_type>(end))
67 {
68 // local entry
69 value = *(ptr + index - begin);
70 }
71 else
72 {
73 // ghost entry
74 Assert(vector.ghost_indices.is_element(index),
76 "You are trying to access an element of a vector "
77 "that is neither a locally owned element nor a "
78 "ghost element of the vector."));
79 const size_type ghostidx =
80 vector.ghost_indices.index_within_set(index);
81
82 AssertIndexRange(ghostidx + end - begin, lsize);
83 value = *(ptr + ghostidx + end - begin);
84 }
85
86 ierr = VecRestoreArrayRead(locally_stored_elements, &ptr);
87 AssertThrow(ierr == 0, ExcPETScError(ierr));
88
89 ierr =
90 VecGhostRestoreLocalForm(vector.vector, &locally_stored_elements);
91 AssertThrow(ierr == 0, ExcPETScError(ierr));
92
93 return value;
94 }
95
96
97 // first verify that the requested
98 // element is actually locally
99 // available
100 PetscInt begin, end;
101
102 PetscErrorCode ierr = VecGetOwnershipRange(vector.vector, &begin, &end);
103 AssertThrow(ierr == 0, ExcPETScError(ierr));
104
105 AssertThrow((index >= static_cast<size_type>(begin)) &&
106 (index < static_cast<size_type>(end)),
107 ExcAccessToNonlocalElement(index, begin, end - 1));
108
109 const PetscScalar *ptr;
110 PetscScalar value;
111 ierr = VecGetArrayRead(vector.vector, &ptr);
112 AssertThrow(ierr == 0, ExcPETScError(ierr));
113 value = *(ptr + index - begin);
114 ierr = VecRestoreArrayRead(vector.vector, &ptr);
115 AssertThrow(ierr == 0, ExcPETScError(ierr));
116
117 return value;
118 }
119# endif
120 } // namespace internal
121
123 : vector(nullptr)
124 , ghosted(false)
125 , last_action(VectorOperation::unknown)
126 {}
127
128
129
131 : ghosted(v.ghosted)
132 , ghost_indices(v.ghost_indices)
133 , last_action(VectorOperation::unknown)
134 {
135 PetscErrorCode ierr = VecDuplicate(v.vector, &vector);
136 AssertThrow(ierr == 0, ExcPETScError(ierr));
137
138 ierr = VecCopy(v.vector, vector);
139 AssertThrow(ierr == 0, ExcPETScError(ierr));
140 }
141
142
143
145 : vector(v)
146 , ghosted(false)
147 , last_action(VectorOperation::unknown)
148 {
149 const PetscErrorCode ierr =
150 PetscObjectReference(reinterpret_cast<PetscObject>(vector));
151 AssertNothrow(ierr == 0, ExcPETScError(ierr));
152 (void)ierr;
154 }
155
156
157
159 {
160 const PetscErrorCode ierr = VecDestroy(&vector);
161 AssertNothrow(ierr == 0, ExcPETScError(ierr));
162 (void)ierr;
163 }
164
165
166
167 void
169 {
171 ExcMessage("Cannot assign a new Vec"));
172 PetscErrorCode ierr =
173 PetscObjectReference(reinterpret_cast<PetscObject>(v));
174 AssertThrow(ierr == 0, ExcPETScError(ierr));
175 ierr = VecDestroy(&vector);
176 AssertThrow(ierr == 0, ExcPETScError(ierr));
177 vector = v;
179 }
180
181
182
183 namespace
184 {
185 template <typename Iterator, typename OutType>
186 class ConvertingIterator
187 {
188 Iterator m_iterator;
189
190 public:
191 using difference_type =
192 typename std::iterator_traits<Iterator>::difference_type;
193 using value_type = OutType;
194 using pointer = OutType *;
195 using reference = OutType &;
196 using iterator_category = std::forward_iterator_tag;
197
198 ConvertingIterator(const Iterator &iterator)
199 : m_iterator(iterator)
200 {}
201
202 OutType
203 operator*() const
204 {
205 return static_cast<OutType>(std::real(*m_iterator));
206 }
207
208 ConvertingIterator &
209 operator++()
210 {
211 ++m_iterator;
212 return *this;
213 }
214
215 ConvertingIterator
216 operator++(int)
217 {
218 ConvertingIterator old = *this;
219 ++m_iterator;
220 return old;
221 }
222
223 bool
224 operator==(const ConvertingIterator &other) const
225 {
226 return this->m_iterator == other.m_iterator;
227 }
228
229 bool
230 operator!=(const ConvertingIterator &other) const
231 {
232 return this->m_iterator != other.m_iterator;
233 }
234 };
235 } // namespace
236
237
238
239 void
241 {
242 // Reset ghost data
243 ghosted = false;
245
246 // There's no API to infer ghost indices from a PETSc Vec which
247 // unfortunately doesn't allow integer entries. We use the
248 // "ConvertingIterator" class above to do an implicit conversion when
249 // sorting and adding ghost indices below.
250 PetscErrorCode ierr;
251 Vec ghosted_vec;
252 ierr = VecGhostGetLocalForm(vector, &ghosted_vec);
253 AssertThrow(ierr == 0, ExcPETScError(ierr));
254 if (ghosted_vec && ghosted_vec != vector)
255 {
256 Vec tvector;
257 PetscScalar *array;
258 PetscInt ghost_start_index, end_index, n_elements_stored_locally;
259
260 ierr = VecGhostRestoreLocalForm(vector, &ghosted_vec);
261 AssertThrow(ierr == 0, ExcPETScError(ierr));
262
263 ierr = VecGetOwnershipRange(vector, &ghost_start_index, &end_index);
264 AssertThrow(ierr == 0, ExcPETScError(ierr));
265 ierr = VecDuplicate(vector, &tvector);
266 AssertThrow(ierr == 0, ExcPETScError(ierr));
267 ierr = VecGetArray(tvector, &array);
268 AssertThrow(ierr == 0, ExcPETScError(ierr));
270 // Store the indices we care about in the vector, so that we can then
271 // exchange this information between processes. It is unfortunate that
272 // we have to store integers in floating point numbers. Let's at least
273 // make sure we do that in a way that ensures that when we get these
274 // numbers back as integers later on, we get the same thing.
275 for (PetscInt i = 0; i < end_index - ghost_start_index; i++)
276 {
277 Assert(static_cast<PetscInt>(std::real(static_cast<PetscScalar>(
278 ghost_start_index + i))) == (ghost_start_index + i),
280 array[i] = ghost_start_index + i;
282
283 ierr = VecRestoreArray(tvector, &array);
284 AssertThrow(ierr == 0, ExcPETScError(ierr));
285 ierr = VecGhostUpdateBegin(tvector, INSERT_VALUES, SCATTER_FORWARD);
286 AssertThrow(ierr == 0, ExcPETScError(ierr));
287 ierr = VecGhostUpdateEnd(tvector, INSERT_VALUES, SCATTER_FORWARD);
288 AssertThrow(ierr == 0, ExcPETScError(ierr));
289 ierr = VecGhostGetLocalForm(tvector, &ghosted_vec);
290 AssertThrow(ierr == 0, ExcPETScError(ierr));
291 ierr = VecGetLocalSize(ghosted_vec, &n_elements_stored_locally);
292 AssertThrow(ierr == 0, ExcPETScError(ierr));
293 ierr = VecGetArrayRead(ghosted_vec, (const PetscScalar **)&array);
294 AssertThrow(ierr == 0, ExcPETScError(ierr));
295
296 // Populate the 'ghosted' flag and the ghost_indices variable. The
297 // latter is an index set that is most efficiently filled by
298 // sorting the indices to add. At the same time, we don't want to
299 // sort the indices stored in a PETSc-owned array; so if the array
300 // is already sorted, pass that to the IndexSet variable, and if
301 // not then copy the indices, sort them, and then add those.
302 ghosted = true;
303 ghost_indices.set_size(this->size());
304
305 ConvertingIterator<PetscScalar *, types::global_dof_index> begin_ghosts(
306 &array[end_index - ghost_start_index]);
307 ConvertingIterator<PetscScalar *, types::global_dof_index> end_ghosts(
308 &array[n_elements_stored_locally]);
309 if (std::is_sorted(&array[end_index - ghost_start_index],
310 &array[n_elements_stored_locally],
311 [](PetscScalar left, PetscScalar right) {
312 return static_cast<PetscInt>(std::real(left)) <
313 static_cast<PetscInt>(std::real(right));
314 }))
315 {
316 ghost_indices.add_indices(begin_ghosts, end_ghosts);
317 }
318 else
319 {
320 std::vector<PetscInt> sorted_indices(begin_ghosts, end_ghosts);
321 std::sort(sorted_indices.begin(), sorted_indices.end());
322 ghost_indices.add_indices(sorted_indices.begin(),
323 sorted_indices.end());
324 }
326
327 ierr = VecRestoreArrayRead(ghosted_vec, (const PetscScalar **)&array);
328 AssertThrow(ierr == 0, ExcPETScError(ierr));
329 ierr = VecGhostRestoreLocalForm(tvector, &ghosted_vec);
330 AssertThrow(ierr == 0, ExcPETScError(ierr));
331 ierr = VecDestroy(&tvector);
332 AssertThrow(ierr == 0, ExcPETScError(ierr));
333 }
334 else
335 {
336 ierr = VecGhostRestoreLocalForm(vector, &ghosted_vec);
337 AssertThrow(ierr == 0, ExcPETScError(ierr));
338 }
339 }
340
341
342 void
344 {
345 const PetscErrorCode ierr = VecDestroy(&vector);
346 AssertThrow(ierr == 0, ExcPETScError(ierr));
347
348 ghosted = false;
351 }
352
353
354
355 VectorBase &
357 {
358 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
359
360 PetscErrorCode ierr = VecCopy(v, vector);
361 AssertThrow(ierr == 0, ExcPETScError(ierr));
362
363 return *this;
364 }
365
366
367
368 VectorBase &
369 VectorBase::operator=(const PetscScalar s)
370 {
372
373 // TODO[TH]: assert(is_compressed())
374
375 PetscErrorCode ierr = VecSet(vector, s);
376 AssertThrow(ierr == 0, ExcPETScError(ierr));
377
378 if (has_ghost_elements())
379 {
380 Vec ghost = nullptr;
381 ierr = VecGhostGetLocalForm(vector, &ghost);
382 AssertThrow(ierr == 0, ExcPETScError(ierr));
383
384 ierr = VecSet(ghost, s);
385 AssertThrow(ierr == 0, ExcPETScError(ierr));
386
387 ierr = VecGhostRestoreLocalForm(vector, &ghost);
388 AssertThrow(ierr == 0, ExcPETScError(ierr));
389 }
390
391 return *this;
392 }
393
394
395
396 bool
398 {
399 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
400
401 PetscBool flag;
402 const PetscErrorCode ierr = VecEqual(vector, v.vector, &flag);
403 AssertThrow(ierr == 0, ExcPETScError(ierr));
404
405 return (flag == PETSC_TRUE);
406 }
407
408
409
410 bool
412 {
413 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
414
415 PetscBool flag;
416 const PetscErrorCode ierr = VecEqual(vector, v.vector, &flag);
417 AssertThrow(ierr == 0, ExcPETScError(ierr));
418
419 return (flag == PETSC_FALSE);
420 }
421
422
423
426 {
427 PetscInt sz;
428 const PetscErrorCode ierr = VecGetSize(vector, &sz);
429 AssertThrow(ierr == 0, ExcPETScError(ierr));
430
431 return sz;
432 }
433
434
435
438 {
439 PetscInt sz;
440 const PetscErrorCode ierr = VecGetLocalSize(vector, &sz);
441 AssertThrow(ierr == 0, ExcPETScError(ierr));
442
443 return sz;
444 }
445
446
447
448 std::pair<VectorBase::size_type, VectorBase::size_type>
450 {
451 PetscInt begin, end;
452 const PetscErrorCode ierr =
453 VecGetOwnershipRange(static_cast<const Vec &>(vector), &begin, &end);
454 AssertThrow(ierr == 0, ExcPETScError(ierr));
455
456 return std::make_pair(begin, end);
457 }
458
459
460
461 void
462 VectorBase::set(const std::vector<size_type> &indices,
463 const std::vector<PetscScalar> &values)
464 {
465 Assert(indices.size() == values.size(),
466 ExcMessage("Function called with arguments of different sizes"));
467 do_set_add_operation(indices.size(), indices.data(), values.data(), false);
468 }
469
470
471
472 void
473 VectorBase::add(const std::vector<size_type> &indices,
474 const std::vector<PetscScalar> &values)
475 {
476 Assert(indices.size() == values.size(),
477 ExcMessage("Function called with arguments of different sizes"));
478 do_set_add_operation(indices.size(), indices.data(), values.data(), true);
479 }
480
481
482
483 void
484 VectorBase::add(const std::vector<size_type> &indices,
485 const ::Vector<PetscScalar> &values)
486 {
487 Assert(indices.size() == values.size(),
488 ExcMessage("Function called with arguments of different sizes"));
489 do_set_add_operation(indices.size(), indices.data(), values.begin(), true);
490 }
491
492
493
494 void
495 VectorBase::add(const size_type n_elements,
496 const size_type *indices,
497 const PetscScalar *values)
498 {
499 do_set_add_operation(n_elements, indices, values, true);
500 }
501
502
503
504 PetscScalar
506 {
507 Assert(size() == vec.size(), ExcDimensionMismatch(size(), vec.size()));
508
509 PetscScalar result;
510
511 // For complex vectors, VecDot() computes
512 // val = (x,y) = y^H x,
513 // where y^H denotes the conjugate transpose of y.
514 // Note that this corresponds to the usual "mathematicians'"
515 // complex inner product where the SECOND argument gets the
516 // complex conjugate, which is also how we document this
517 // operation.
518 const PetscErrorCode ierr = VecDot(vec.vector, vector, &result);
519 AssertThrow(ierr == 0, ExcPETScError(ierr));
520
521 return result;
522 }
523
524
525
526 PetscScalar
527 VectorBase::add_and_dot(const PetscScalar a,
528 const VectorBase &V,
529 const VectorBase &W)
530 {
531 this->add(a, V);
532 return *this * W;
533 }
534
535
536
537 void
539 {
540 Assert(has_ghost_elements() == false,
541 ExcMessage("Calling compress() is only useful if a vector "
542 "has been written into, but this is a vector with ghost "
543 "elements and consequently is read-only. It does "
544 "not make sense to call compress() for such "
545 "vectors."));
546
547 {
548 if constexpr (running_in_debug_mode())
549 {
550 // Check that all processors agree that last_action is the same (or
551 // none!)
552
553 int my_int_last_action = last_action;
554 int all_int_last_action;
555
556 const int ierr = MPI_Allreduce(&my_int_last_action,
557 &all_int_last_action,
558 1,
559 MPI_INT,
560 MPI_BOR,
562 AssertThrowMPI(ierr);
563
564 AssertThrow(all_int_last_action !=
567 "Error: not all processors agree on the last "
568 "VectorOperation before this compress() call."));
569 }
570 }
571
575 "Missing compress() or calling with wrong VectorOperation argument."));
576
577 // note that one may think that
578 // we only need to do something
579 // if in fact the state is
580 // anything but
581 // last_action::unknown. but
582 // that's not true: one
583 // frequently gets into
584 // situations where only one
585 // processor (or a subset of
586 // processors) actually writes
587 // something into a vector, but
588 // we still need to call
589 // VecAssemblyBegin/End on all
590 // processors.
591 PetscErrorCode ierr = VecAssemblyBegin(vector);
592 AssertThrow(ierr == 0, ExcPETScError(ierr));
593 ierr = VecAssemblyEnd(vector);
594 AssertThrow(ierr == 0, ExcPETScError(ierr));
595
596 // reset the last action field to
597 // indicate that we're back to a
598 // pristine state
600 }
601
602
603
606 {
607 const real_type d = l2_norm();
608 return d * d;
609 }
610
611
612
613 PetscScalar
615 {
616 // We can only use our more efficient
617 // routine in the serial case.
618 if (dynamic_cast<const PETScWrappers::MPI::Vector *>(this) != nullptr)
619 {
620 PetscScalar sum;
621 const PetscErrorCode ierr = VecSum(vector, &sum);
622 AssertThrow(ierr == 0, ExcPETScError(ierr));
623 return sum / static_cast<PetscReal>(size());
624 }
625
626 // get a representation of the vector and
627 // loop over all the elements
628 const PetscScalar *start_ptr;
629 PetscErrorCode ierr = VecGetArrayRead(vector, &start_ptr);
630 AssertThrow(ierr == 0, ExcPETScError(ierr));
631
632 PetscScalar mean = 0;
633 {
634 PetscScalar sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0;
635
636 // use modern processors better by
637 // allowing pipelined commands to be
638 // executed in parallel
639 const PetscScalar *ptr = start_ptr;
640 const PetscScalar *eptr = ptr + (locally_owned_size() / 4) * 4;
641 while (ptr != eptr)
642 {
643 sum0 += *ptr++;
644 sum1 += *ptr++;
645 sum2 += *ptr++;
646 sum3 += *ptr++;
647 }
648 // add up remaining elements
649 while (ptr != start_ptr + locally_owned_size())
650 sum0 += *ptr++;
651
652 mean =
653 Utilities::MPI::sum(sum0 + sum1 + sum2 + sum3, get_mpi_communicator()) /
654 static_cast<PetscReal>(size());
655 }
656
657 // restore the representation of the
658 // vector
659 ierr = VecRestoreArrayRead(vector, &start_ptr);
660 AssertThrow(ierr == 0, ExcPETScError(ierr));
661
662 return mean;
663 }
664
665
668 {
669 real_type d;
670
671 const PetscErrorCode ierr = VecNorm(vector, NORM_1, &d);
672 AssertThrow(ierr == 0, ExcPETScError(ierr));
673
674 return d;
675 }
676
677
678
681 {
682 real_type d;
683
684 const PetscErrorCode ierr = VecNorm(vector, NORM_2, &d);
685 AssertThrow(ierr == 0, ExcPETScError(ierr));
686
687 return d;
688 }
689
690
691
694 {
695 // get a representation of the vector and
696 // loop over all the elements
697 const PetscScalar *start_ptr;
698 PetscErrorCode ierr = VecGetArrayRead(vector, &start_ptr);
699 AssertThrow(ierr == 0, ExcPETScError(ierr));
700
701 real_type norm = 0;
702 {
703 real_type sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0;
704
705 // use modern processors better by
706 // allowing pipelined commands to be
707 // executed in parallel
708 const PetscScalar *ptr = start_ptr;
709 const PetscScalar *eptr = ptr + (locally_owned_size() / 4) * 4;
710 while (ptr != eptr)
711 {
716 }
717 // add up remaining elements
718 while (ptr != start_ptr + locally_owned_size())
720
721 norm = std::pow(Utilities::MPI::sum(sum0 + sum1 + sum2 + sum3,
723 1. / p);
724 }
725
726 // restore the representation of the
727 // vector
728 ierr = VecRestoreArrayRead(vector, &start_ptr);
729 AssertThrow(ierr == 0, ExcPETScError(ierr));
730
731 return norm;
732 }
733
734
735
738 {
739 real_type d;
740
741 const PetscErrorCode ierr = VecNorm(vector, NORM_INFINITY, &d);
742 AssertThrow(ierr == 0, ExcPETScError(ierr));
743
744 return d;
745 }
746
747
748
749 bool
751 {
752 // get a representation of the vector and
753 // loop over all the elements
754 const PetscScalar *start_ptr;
755 PetscErrorCode ierr = VecGetArrayRead(vector, &start_ptr);
756 AssertThrow(ierr == 0, ExcPETScError(ierr));
757
758 const PetscScalar *ptr = start_ptr,
759 *eptr = start_ptr + locally_owned_size();
760 bool flag = true;
761 while (ptr != eptr)
762 {
763 if (*ptr != value_type())
764 {
765 flag = false;
766 break;
767 }
768 ++ptr;
769 }
770
771 // restore the representation of the
772 // vector
773 ierr = VecRestoreArrayRead(vector, &start_ptr);
774 AssertThrow(ierr == 0, ExcPETScError(ierr));
775
776 return flag;
777 }
778
779
780 namespace internal
781 {
782 template <typename T>
783 bool
784 is_non_negative(const T &t)
785 {
786 return t >= 0;
787 }
788
789
790
791 template <typename T>
792 bool
793 is_non_negative(const std::complex<T> &)
794 {
795 Assert(false,
796 ExcMessage("You can't ask a complex value "
797 "whether it is non-negative."));
798 return true;
799 }
800 } // namespace internal
801
802
803
804 VectorBase &
805 VectorBase::operator*=(const PetscScalar a)
806 {
809
810 const PetscErrorCode ierr = VecScale(vector, a);
811 AssertThrow(ierr == 0, ExcPETScError(ierr));
812
813 return *this;
814 }
815
816
817
818 VectorBase &
819 VectorBase::operator/=(const PetscScalar a)
820 {
823
824 const PetscScalar factor = 1. / a;
825 AssertIsFinite(factor);
826
827 const PetscErrorCode ierr = VecScale(vector, factor);
828 AssertThrow(ierr == 0, ExcPETScError(ierr));
829
830 return *this;
831 }
832
833
834
835 VectorBase &
837 {
839 const PetscErrorCode ierr = VecAXPY(vector, 1, v);
840 AssertThrow(ierr == 0, ExcPETScError(ierr));
841
842 return *this;
843 }
844
845
846
847 VectorBase &
849 {
851 const PetscErrorCode ierr = VecAXPY(vector, -1, v);
852 AssertThrow(ierr == 0, ExcPETScError(ierr));
853
854 return *this;
855 }
856
857
858
859 void
860 VectorBase::add(const PetscScalar s)
861 {
864
865 const PetscErrorCode ierr = VecShift(vector, s);
866 AssertThrow(ierr == 0, ExcPETScError(ierr));
867 }
868
869
870
871 void
872 VectorBase::add(const PetscScalar a, const VectorBase &v)
873 {
876
877 const PetscErrorCode ierr = VecAXPY(vector, a, v);
878 AssertThrow(ierr == 0, ExcPETScError(ierr));
879 }
880
881
882
883 void
884 VectorBase::add(const PetscScalar a,
885 const VectorBase &v,
886 const PetscScalar b,
887 const VectorBase &w)
888 {
892
893 const PetscScalar weights[2] = {a, b};
894 Vec addends[2] = {v.vector, w.vector};
895
896 const PetscErrorCode ierr = VecMAXPY(vector, 2, weights, addends);
897 AssertThrow(ierr == 0, ExcPETScError(ierr));
898 }
899
900
901
902 void
903 VectorBase::sadd(const PetscScalar s, const VectorBase &v)
904 {
907
908 const PetscErrorCode ierr = VecAYPX(vector, s, v);
909 AssertThrow(ierr == 0, ExcPETScError(ierr));
910 }
911
912
913
914 void
915 VectorBase::sadd(const PetscScalar s,
916 const PetscScalar a,
917 const VectorBase &v)
918 {
922
923 // there is nothing like a AXPAY
924 // operation in PETSc, so do it in two
925 // steps
926 *this *= s;
927 add(a, v);
928 }
929
930
931
932 void
934 {
936 const PetscErrorCode ierr = VecPointwiseMult(vector, factors, vector);
937 AssertThrow(ierr == 0, ExcPETScError(ierr));
938 }
939
940
941
942 void
943 VectorBase::equ(const PetscScalar a, const VectorBase &v)
944 {
947
948 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
949
950 const PetscErrorCode ierr = VecAXPBY(vector, a, 0.0, v.vector);
951 AssertThrow(ierr == 0, ExcPETScError(ierr));
952 }
953
954
955
956 void
957 VectorBase::write_ascii(const PetscViewerFormat format)
958 {
959 // TODO[TH]:assert(is_compressed())
960 MPI_Comm comm = PetscObjectComm((PetscObject)vector);
961
962 // Set options
963 PetscErrorCode ierr =
964 PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(comm), format);
965 AssertThrow(ierr == 0, ExcPETScError(ierr));
966
967 // Write to screen
968 ierr = VecView(vector, PETSC_VIEWER_STDOUT_(comm));
969 AssertThrow(ierr == 0, ExcPETScError(ierr));
970 ierr = PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(comm));
971 AssertThrow(ierr == 0, ExcPETScError(ierr));
972 }
973
974
975
976 void
977 VectorBase::print(std::ostream &out,
978 const unsigned int precision,
979 const bool scientific,
980 const bool across) const
981 {
982 AssertThrow(out.fail() == false, ExcIO());
983
984 // get a representation of the vector and
985 // loop over all the elements
986 const PetscScalar *val;
987 PetscErrorCode ierr = VecGetArrayRead(vector, &val);
988
989 AssertThrow(ierr == 0, ExcPETScError(ierr));
990
991 // save the state of out stream
992 const std::ios::fmtflags old_flags = out.flags();
993 const unsigned int old_precision = out.precision(precision);
994
995 out.precision(precision);
996 if (scientific)
997 out.setf(std::ios::scientific, std::ios::floatfield);
998 else
999 out.setf(std::ios::fixed, std::ios::floatfield);
1000
1001 if (across)
1002 for (size_type i = 0; i < locally_owned_size(); ++i)
1003 out << val[i] << ' ';
1004 else
1005 for (size_type i = 0; i < locally_owned_size(); ++i)
1006 out << val[i] << std::endl;
1007 out << std::endl;
1008
1009 // reset output format
1010 out.flags(old_flags);
1011 out.precision(old_precision);
1012
1013 // restore the representation of the
1014 // vector
1015 ierr = VecRestoreArrayRead(vector, &val);
1016 AssertThrow(ierr == 0, ExcPETScError(ierr));
1017
1018 AssertThrow(out.fail() == false, ExcIO());
1019 }
1020
1021
1022
1023 void
1025 {
1026 std::swap(this->vector, v.vector);
1027 std::swap(this->ghosted, v.ghosted);
1028 std::swap(this->last_action, v.last_action);
1029 // missing swap for IndexSet
1030 IndexSet t(this->ghost_indices);
1031 this->ghost_indices = v.ghost_indices;
1032 v.ghost_indices = t;
1033 }
1034
1035
1036 VectorBase::operator const Vec &() const
1037 {
1038 return vector;
1039 }
1040
1041
1042 Vec &
1044 {
1045 return vector;
1046 }
1047
1048
1049 std::size_t
1051 {
1052 std::size_t mem = sizeof(Vec) + sizeof(last_action) +
1055
1056 // TH: I am relatively sure that PETSc is
1057 // storing the local data in a contiguous
1058 // block without indices:
1059 mem += locally_owned_size() * sizeof(PetscScalar);
1060 // assume that PETSc is storing one index
1061 // and one double per ghost element
1062 if (ghosted)
1063 mem += ghost_indices.n_elements() * (sizeof(PetscScalar) + sizeof(int));
1064
1065 // TODO[TH]: size of constant memory for PETSc?
1066 return mem;
1067 }
1068
1069
1070
1071 void
1073 const size_type *indices,
1074 const PetscScalar *values,
1075 const bool add_values)
1076 {
1080 internal::VectorReference::ExcWrongMode(action, last_action));
1082
1083 std::vector<PetscInt> petsc_indices(n_elements);
1084 for (size_type i = 0; i < n_elements; ++i)
1085 {
1086 const auto petsc_index = static_cast<PetscInt>(indices[i]);
1087 AssertIntegerConversion(petsc_index, indices[i]);
1088 petsc_indices[i] = petsc_index;
1089 }
1090
1091 const InsertMode mode = (add_values ? ADD_VALUES : INSERT_VALUES);
1092 const PetscErrorCode ierr = VecSetValues(
1093 vector, petsc_indices.size(), petsc_indices.data(), values, mode);
1094 AssertThrow(ierr == 0, ExcPETScError(ierr));
1095
1096 last_action = action;
1097 }
1098
1099} // namespace PETScWrappers
1100
1102
1103#endif // DEAL_II_WITH_PETSC
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
size_type n_elements() const
Definition index_set.h:1945
void set_size(const size_type size)
Definition index_set.h:1775
void clear()
Definition index_set.h:1763
void compress() const
Definition index_set.h:1795
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1842
real_type lp_norm(const real_type p) const
VectorBase & operator+=(const VectorBase &V)
PetscScalar mean_value() const
VectorOperation::values last_action
PetscScalar operator*(const VectorBase &vec) const
bool operator==(const VectorBase &v) const
VectorBase & operator*=(const PetscScalar factor)
VectorBase & operator-=(const VectorBase &V)
bool operator!=(const VectorBase &v) const
std::size_t memory_consumption() const
VectorBase & operator/=(const PetscScalar factor)
MPI_Comm get_mpi_communicator() const
void scale(const VectorBase &scaling_factors)
std::pair< size_type, size_type > local_range() const
void sadd(const PetscScalar s, const VectorBase &V)
void compress(const VectorOperation::values operation)
void set(const std::vector< size_type > &indices, const std::vector< PetscScalar > &values)
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
PetscScalar add_and_dot(const PetscScalar a, const VectorBase &V, const VectorBase &W)
void swap(VectorBase &v) noexcept
void add(const std::vector< size_type > &indices, const std::vector< PetscScalar > &values)
void write_ascii(const PetscViewerFormat format=PETSC_VIEWER_DEFAULT)
bool has_ghost_elements() const
size_type locally_owned_size() const
void equ(const PetscScalar a, const VectorBase &V)
void do_set_add_operation(const size_type n_elements, const size_type *indices, const PetscScalar *values, const bool add_values)
VectorBase & operator=(const VectorBase &)
size_type size() const override
std::enable_if_t< std::is_floating_point_v< T > &&std::is_floating_point_v< U > &&!std::is_same_v< T, U >, typename ProductType< std::complex< T >, std::complex< U > >::type > operator*(const std::complex< T > &left, const std::complex< U > &right)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
constexpr bool running_in_debug_mode()
Definition config.h:78
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
#define AssertIntegerConversion(index1, index2)
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcGhostsPresent()
#define Assert(cond, exc)
#define AssertIsFinite(number)
#define AssertThrowMPI(error_code)
#define AssertNothrow(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
const MPI_Comm comm
Definition mpi.cc:922
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
bool is_non_negative(const T &t)
T sum(const T &t, const MPI_Comm mpi_communicator)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
Iterator m_iterator
SynchronousIterators< Iterators > & operator++(SynchronousIterators< Iterators > &a)