deal.II version GIT relicensing-3512-g0a98d4ed9f 2025-06-14 14:10: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));
153 }
154
155
156
158 {
159 const PetscErrorCode ierr = VecDestroy(&vector);
160 AssertNothrow(ierr == 0, ExcPETScError(ierr));
161 }
162
163
164
165 void
167 {
169 ExcMessage("Cannot assign a new Vec"));
170 PetscErrorCode ierr =
171 PetscObjectReference(reinterpret_cast<PetscObject>(v));
172 AssertThrow(ierr == 0, ExcPETScError(ierr));
173 ierr = VecDestroy(&vector);
174 AssertThrow(ierr == 0, ExcPETScError(ierr));
175 vector = v;
177 }
178
179
180
181 namespace
182 {
183 template <typename Iterator, typename OutType>
184 class ConvertingIterator
185 {
186 Iterator m_iterator;
187
188 public:
189 using difference_type =
190 typename std::iterator_traits<Iterator>::difference_type;
191 using value_type = OutType;
192 using pointer = OutType *;
193 using reference = OutType &;
194 using iterator_category = std::forward_iterator_tag;
195
196 ConvertingIterator(const Iterator &iterator)
197 : m_iterator(iterator)
198 {}
199
200 OutType
201 operator*() const
202 {
203 return static_cast<OutType>(std::real(*m_iterator));
204 }
205
206 ConvertingIterator &
207 operator++()
208 {
209 ++m_iterator;
210 return *this;
211 }
212
213 ConvertingIterator
214 operator++(int)
215 {
216 ConvertingIterator old = *this;
217 ++m_iterator;
218 return old;
219 }
220
221 bool
222 operator==(const ConvertingIterator &other) const
223 {
224 return this->m_iterator == other.m_iterator;
225 }
226
227 bool
228 operator!=(const ConvertingIterator &other) const
229 {
230 return this->m_iterator != other.m_iterator;
231 }
232 };
233 } // namespace
234
235
236
237 void
239 {
240 // Reset ghost data
241 ghosted = false;
243
244 // There's no API to infer ghost indices from a PETSc Vec which
245 // unfortunately doesn't allow integer entries. We use the
246 // "ConvertingIterator" class above to do an implicit conversion when
247 // sorting and adding ghost indices below.
248 PetscErrorCode ierr;
249 Vec ghosted_vec;
250 ierr = VecGhostGetLocalForm(vector, &ghosted_vec);
251 AssertThrow(ierr == 0, ExcPETScError(ierr));
252 if (ghosted_vec && ghosted_vec != vector)
253 {
254 Vec tvector;
255 PetscScalar *array;
256 PetscInt ghost_start_index, end_index, n_elements_stored_locally;
257
258 ierr = VecGhostRestoreLocalForm(vector, &ghosted_vec);
259 AssertThrow(ierr == 0, ExcPETScError(ierr));
260
261 ierr = VecGetOwnershipRange(vector, &ghost_start_index, &end_index);
262 AssertThrow(ierr == 0, ExcPETScError(ierr));
263 ierr = VecDuplicate(vector, &tvector);
264 AssertThrow(ierr == 0, ExcPETScError(ierr));
265 ierr = VecGetArray(tvector, &array);
266 AssertThrow(ierr == 0, ExcPETScError(ierr));
267
268 // Store the indices we care about in the vector, so that we can then
269 // exchange this information between processes. It is unfortunate that
270 // we have to store integers in floating point numbers. Let's at least
271 // make sure we do that in a way that ensures that when we get these
272 // numbers back as integers later on, we get the same thing.
273 for (PetscInt i = 0; i < end_index - ghost_start_index; i++)
274 {
275 Assert(static_cast<PetscInt>(std::real(static_cast<PetscScalar>(
276 ghost_start_index + i))) == (ghost_start_index + i),
278 array[i] = ghost_start_index + i;
279 }
280
281 ierr = VecRestoreArray(tvector, &array);
282 AssertThrow(ierr == 0, ExcPETScError(ierr));
283 ierr = VecGhostUpdateBegin(tvector, INSERT_VALUES, SCATTER_FORWARD);
284 AssertThrow(ierr == 0, ExcPETScError(ierr));
285 ierr = VecGhostUpdateEnd(tvector, INSERT_VALUES, SCATTER_FORWARD);
286 AssertThrow(ierr == 0, ExcPETScError(ierr));
287 ierr = VecGhostGetLocalForm(tvector, &ghosted_vec);
288 AssertThrow(ierr == 0, ExcPETScError(ierr));
289 ierr = VecGetLocalSize(ghosted_vec, &n_elements_stored_locally);
290 AssertThrow(ierr == 0, ExcPETScError(ierr));
291 ierr = VecGetArrayRead(ghosted_vec, (const PetscScalar **)&array);
292 AssertThrow(ierr == 0, ExcPETScError(ierr));
293
294 // Populate the 'ghosted' flag and the ghost_indices variable. The
295 // latter is an index set that is most efficiently filled by
296 // sorting the indices to add. At the same time, we don't want to
297 // sort the indices stored in a PETSc-owned array; so if the array
298 // is already sorted, pass that to the IndexSet variable, and if
299 // not then copy the indices, sort them, and then add those.
300 ghosted = true;
301 ghost_indices.set_size(this->size());
302
303 ConvertingIterator<PetscScalar *, types::global_dof_index> begin_ghosts(
304 &array[end_index - ghost_start_index]);
305 ConvertingIterator<PetscScalar *, types::global_dof_index> end_ghosts(
306 &array[n_elements_stored_locally]);
307 if (std::is_sorted(&array[end_index - ghost_start_index],
308 &array[n_elements_stored_locally],
309 [](PetscScalar left, PetscScalar right) {
310 return static_cast<PetscInt>(std::real(left)) <
311 static_cast<PetscInt>(std::real(right));
312 }))
313 {
314 ghost_indices.add_indices(begin_ghosts, end_ghosts);
315 }
316 else
317 {
318 std::vector<PetscInt> sorted_indices(begin_ghosts, end_ghosts);
319 std::sort(sorted_indices.begin(), sorted_indices.end());
320 ghost_indices.add_indices(sorted_indices.begin(),
321 sorted_indices.end());
322 }
324
325 ierr = VecRestoreArrayRead(ghosted_vec, (const PetscScalar **)&array);
326 AssertThrow(ierr == 0, ExcPETScError(ierr));
327 ierr = VecGhostRestoreLocalForm(tvector, &ghosted_vec);
328 AssertThrow(ierr == 0, ExcPETScError(ierr));
329 ierr = VecDestroy(&tvector);
330 AssertThrow(ierr == 0, ExcPETScError(ierr));
331 }
332 else
333 {
334 ierr = VecGhostRestoreLocalForm(vector, &ghosted_vec);
335 AssertThrow(ierr == 0, ExcPETScError(ierr));
336 }
338
339
340 void
342 {
343 const PetscErrorCode ierr = VecDestroy(&vector);
344 AssertThrow(ierr == 0, ExcPETScError(ierr));
345
346 ghosted = false;
349 }
350
351
352
353 VectorBase &
355 {
356 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
357
358 PetscErrorCode ierr = VecCopy(v, vector);
359 AssertThrow(ierr == 0, ExcPETScError(ierr));
360
361 return *this;
362 }
363
364
365
366 VectorBase &
367 VectorBase::operator=(const PetscScalar s)
368 {
370
371 // TODO[TH]: assert(is_compressed())
372
373 PetscErrorCode ierr = VecSet(vector, s);
374 AssertThrow(ierr == 0, ExcPETScError(ierr));
375
376 if (has_ghost_elements())
377 {
378 Vec ghost = nullptr;
379 ierr = VecGhostGetLocalForm(vector, &ghost);
380 AssertThrow(ierr == 0, ExcPETScError(ierr));
381
382 ierr = VecSet(ghost, s);
383 AssertThrow(ierr == 0, ExcPETScError(ierr));
384
385 ierr = VecGhostRestoreLocalForm(vector, &ghost);
386 AssertThrow(ierr == 0, ExcPETScError(ierr));
387 }
388
389 return *this;
390 }
391
392
393
394 bool
396 {
397 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
398
399 PetscBool flag;
400 const PetscErrorCode ierr = VecEqual(vector, v.vector, &flag);
401 AssertThrow(ierr == 0, ExcPETScError(ierr));
402
403 return (flag == PETSC_TRUE);
404 }
405
406
407
408 bool
410 {
411 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
412
413 PetscBool flag;
414 const PetscErrorCode ierr = VecEqual(vector, v.vector, &flag);
415 AssertThrow(ierr == 0, ExcPETScError(ierr));
416
417 return (flag == PETSC_FALSE);
418 }
419
420
421
424 {
425 PetscInt sz;
426 const PetscErrorCode ierr = VecGetSize(vector, &sz);
427 AssertThrow(ierr == 0, ExcPETScError(ierr));
428
429 return sz;
430 }
431
432
433
436 {
437 PetscInt sz;
438 const PetscErrorCode ierr = VecGetLocalSize(vector, &sz);
439 AssertThrow(ierr == 0, ExcPETScError(ierr));
440
441 return sz;
442 }
443
444
445
446 std::pair<VectorBase::size_type, VectorBase::size_type>
448 {
449 PetscInt begin, end;
450 const PetscErrorCode ierr =
451 VecGetOwnershipRange(static_cast<const Vec &>(vector), &begin, &end);
452 AssertThrow(ierr == 0, ExcPETScError(ierr));
453
454 return std::make_pair(begin, end);
455 }
456
457
458
459 void
460 VectorBase::set(const std::vector<size_type> &indices,
461 const std::vector<PetscScalar> &values)
462 {
463 Assert(indices.size() == values.size(),
464 ExcMessage("Function called with arguments of different sizes"));
465 do_set_add_operation(indices.size(), indices.data(), values.data(), false);
466 }
467
468
469
470 void
471 VectorBase::add(const std::vector<size_type> &indices,
472 const std::vector<PetscScalar> &values)
473 {
474 Assert(indices.size() == values.size(),
475 ExcMessage("Function called with arguments of different sizes"));
476 do_set_add_operation(indices.size(), indices.data(), values.data(), true);
477 }
478
479
480
481 void
482 VectorBase::add(const std::vector<size_type> &indices,
483 const ::Vector<PetscScalar> &values)
484 {
485 Assert(indices.size() == values.size(),
486 ExcMessage("Function called with arguments of different sizes"));
487 do_set_add_operation(indices.size(), indices.data(), values.begin(), true);
488 }
489
490
491
492 void
493 VectorBase::add(const size_type n_elements,
494 const size_type *indices,
495 const PetscScalar *values)
496 {
497 do_set_add_operation(n_elements, indices, values, true);
498 }
499
500
501
502 PetscScalar
504 {
505 Assert(size() == vec.size(), ExcDimensionMismatch(size(), vec.size()));
506
507 PetscScalar result;
508
509 // For complex vectors, VecDot() computes
510 // val = (x,y) = y^H x,
511 // where y^H denotes the conjugate transpose of y.
512 // Note that this corresponds to the usual "mathematicians'"
513 // complex inner product where the SECOND argument gets the
514 // complex conjugate, which is also how we document this
515 // operation.
516 const PetscErrorCode ierr = VecDot(vec.vector, vector, &result);
517 AssertThrow(ierr == 0, ExcPETScError(ierr));
518
519 return result;
520 }
521
522
523
524 PetscScalar
525 VectorBase::add_and_dot(const PetscScalar a,
526 const VectorBase &V,
527 const VectorBase &W)
528 {
529 this->add(a, V);
530 return *this * W;
531 }
532
533
534
535 void
537 {
538 Assert(has_ghost_elements() == false,
539 ExcMessage("Calling compress() is only useful if a vector "
540 "has been written into, but this is a vector with ghost "
541 "elements and consequently is read-only. It does "
542 "not make sense to call compress() for such "
543 "vectors."));
544
545 {
546 if constexpr (running_in_debug_mode())
547 {
548 // Check that all processors agree that last_action is the same (or
549 // none!)
550
551 int my_int_last_action = last_action;
552 int all_int_last_action;
553
554 const int ierr = MPI_Allreduce(&my_int_last_action,
555 &all_int_last_action,
556 1,
557 MPI_INT,
558 MPI_BOR,
560 AssertThrowMPI(ierr);
561
562 AssertThrow(all_int_last_action !=
565 "Error: not all processors agree on the last "
566 "VectorOperation before this compress() call."));
567 }
568 }
569
573 "Missing compress() or calling with wrong VectorOperation argument."));
574
575 // note that one may think that
576 // we only need to do something
577 // if in fact the state is
578 // anything but
579 // last_action::unknown. but
580 // that's not true: one
581 // frequently gets into
582 // situations where only one
583 // processor (or a subset of
584 // processors) actually writes
585 // something into a vector, but
586 // we still need to call
587 // VecAssemblyBegin/End on all
588 // processors.
589 PetscErrorCode ierr = VecAssemblyBegin(vector);
590 AssertThrow(ierr == 0, ExcPETScError(ierr));
591 ierr = VecAssemblyEnd(vector);
592 AssertThrow(ierr == 0, ExcPETScError(ierr));
593
594 // reset the last action field to
595 // indicate that we're back to a
596 // pristine state
598 }
599
600
601
604 {
605 const real_type d = l2_norm();
606 return d * d;
607 }
608
609
610
611 PetscScalar
613 {
614 // We can only use our more efficient
615 // routine in the serial case.
616 if (dynamic_cast<const PETScWrappers::MPI::Vector *>(this) != nullptr)
617 {
618 PetscScalar sum;
619 const PetscErrorCode ierr = VecSum(vector, &sum);
620 AssertThrow(ierr == 0, ExcPETScError(ierr));
621 return sum / static_cast<PetscReal>(size());
622 }
623
624 // get a representation of the vector and
625 // loop over all the elements
626 const PetscScalar *start_ptr;
627 PetscErrorCode ierr = VecGetArrayRead(vector, &start_ptr);
628 AssertThrow(ierr == 0, ExcPETScError(ierr));
629
630 PetscScalar mean = 0;
631 {
632 PetscScalar sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0;
633
634 // use modern processors better by
635 // allowing pipelined commands to be
636 // executed in parallel
637 const PetscScalar *ptr = start_ptr;
638 const PetscScalar *eptr = ptr + (locally_owned_size() / 4) * 4;
639 while (ptr != eptr)
640 {
641 sum0 += *ptr++;
642 sum1 += *ptr++;
643 sum2 += *ptr++;
644 sum3 += *ptr++;
645 }
646 // add up remaining elements
647 while (ptr != start_ptr + locally_owned_size())
648 sum0 += *ptr++;
649
650 mean =
651 Utilities::MPI::sum(sum0 + sum1 + sum2 + sum3, get_mpi_communicator()) /
652 static_cast<PetscReal>(size());
653 }
654
655 // restore the representation of the
656 // vector
657 ierr = VecRestoreArrayRead(vector, &start_ptr);
658 AssertThrow(ierr == 0, ExcPETScError(ierr));
659
660 return mean;
661 }
662
663
666 {
667 real_type d;
668
669 const PetscErrorCode ierr = VecNorm(vector, NORM_1, &d);
670 AssertThrow(ierr == 0, ExcPETScError(ierr));
671
672 return d;
673 }
674
675
676
679 {
680 real_type d;
681
682 const PetscErrorCode ierr = VecNorm(vector, NORM_2, &d);
683 AssertThrow(ierr == 0, ExcPETScError(ierr));
684
685 return d;
686 }
687
688
689
692 {
693 // get a representation of the vector and
694 // loop over all the elements
695 const PetscScalar *start_ptr;
696 PetscErrorCode ierr = VecGetArrayRead(vector, &start_ptr);
697 AssertThrow(ierr == 0, ExcPETScError(ierr));
698
699 real_type norm = 0;
700 {
701 real_type sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0;
702
703 // use modern processors better by
704 // allowing pipelined commands to be
705 // executed in parallel
706 const PetscScalar *ptr = start_ptr;
707 const PetscScalar *eptr = ptr + (locally_owned_size() / 4) * 4;
708 while (ptr != eptr)
709 {
714 }
715 // add up remaining elements
716 while (ptr != start_ptr + locally_owned_size())
718
719 norm = std::pow(Utilities::MPI::sum(sum0 + sum1 + sum2 + sum3,
721 1. / p);
722 }
723
724 // restore the representation of the
725 // vector
726 ierr = VecRestoreArrayRead(vector, &start_ptr);
727 AssertThrow(ierr == 0, ExcPETScError(ierr));
728
729 return norm;
730 }
731
732
733
736 {
737 real_type d;
738
739 const PetscErrorCode ierr = VecNorm(vector, NORM_INFINITY, &d);
740 AssertThrow(ierr == 0, ExcPETScError(ierr));
741
742 return d;
743 }
744
745
746
747 bool
749 {
750 // get a representation of the vector and
751 // loop over all the elements
752 const PetscScalar *start_ptr;
753 PetscErrorCode ierr = VecGetArrayRead(vector, &start_ptr);
754 AssertThrow(ierr == 0, ExcPETScError(ierr));
755
756 const PetscScalar *ptr = start_ptr,
757 *eptr = start_ptr + locally_owned_size();
758 bool flag = true;
759 while (ptr != eptr)
760 {
761 if (*ptr != value_type())
762 {
763 flag = false;
764 break;
765 }
766 ++ptr;
767 }
768
769 // restore the representation of the
770 // vector
771 ierr = VecRestoreArrayRead(vector, &start_ptr);
772 AssertThrow(ierr == 0, ExcPETScError(ierr));
773
774 return flag;
775 }
776
777
778 namespace internal
779 {
780 template <typename T>
781 bool
782 is_non_negative(const T &t)
783 {
784 return t >= 0;
785 }
786
787
788
789 template <typename T>
790 bool
791 is_non_negative(const std::complex<T> &)
792 {
793 Assert(false,
794 ExcMessage("You can't ask a complex value "
795 "whether it is non-negative."));
796 return true;
797 }
798 } // namespace internal
799
800
801
802 VectorBase &
803 VectorBase::operator*=(const PetscScalar a)
804 {
807
808 const PetscErrorCode ierr = VecScale(vector, a);
809 AssertThrow(ierr == 0, ExcPETScError(ierr));
810
811 return *this;
812 }
813
814
815
816 VectorBase &
817 VectorBase::operator/=(const PetscScalar a)
818 {
821
822 const PetscScalar factor = 1. / a;
823 AssertIsFinite(factor);
824
825 const PetscErrorCode ierr = VecScale(vector, factor);
826 AssertThrow(ierr == 0, ExcPETScError(ierr));
827
828 return *this;
829 }
830
831
832
833 VectorBase &
835 {
837 const PetscErrorCode ierr = VecAXPY(vector, 1, v);
838 AssertThrow(ierr == 0, ExcPETScError(ierr));
839
840 return *this;
841 }
842
843
844
845 VectorBase &
847 {
849 const PetscErrorCode ierr = VecAXPY(vector, -1, v);
850 AssertThrow(ierr == 0, ExcPETScError(ierr));
851
852 return *this;
853 }
854
855
856
857 void
858 VectorBase::add(const PetscScalar s)
859 {
862
863 const PetscErrorCode ierr = VecShift(vector, s);
864 AssertThrow(ierr == 0, ExcPETScError(ierr));
865 }
866
867
868
869 void
870 VectorBase::add(const PetscScalar a, const VectorBase &v)
871 {
874
875 const PetscErrorCode ierr = VecAXPY(vector, a, v);
876 AssertThrow(ierr == 0, ExcPETScError(ierr));
877 }
878
879
880
881 void
882 VectorBase::add(const PetscScalar a,
883 const VectorBase &v,
884 const PetscScalar b,
885 const VectorBase &w)
886 {
890
891 const PetscScalar weights[2] = {a, b};
892 Vec addends[2] = {v.vector, w.vector};
893
894 const PetscErrorCode ierr = VecMAXPY(vector, 2, weights, addends);
895 AssertThrow(ierr == 0, ExcPETScError(ierr));
896 }
897
898
899
900 void
901 VectorBase::sadd(const PetscScalar s, const VectorBase &v)
902 {
905
906 const PetscErrorCode ierr = VecAYPX(vector, s, v);
907 AssertThrow(ierr == 0, ExcPETScError(ierr));
908 }
909
910
911
912 void
913 VectorBase::sadd(const PetscScalar s,
914 const PetscScalar a,
915 const VectorBase &v)
916 {
920
921 // there is nothing like a AXPAY
922 // operation in PETSc, so do it in two
923 // steps
924 *this *= s;
925 add(a, v);
926 }
927
928
929
930 void
932 {
934 const PetscErrorCode ierr = VecPointwiseMult(vector, factors, vector);
935 AssertThrow(ierr == 0, ExcPETScError(ierr));
936 }
937
938
939
940 void
941 VectorBase::equ(const PetscScalar a, const VectorBase &v)
942 {
945
946 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
947
948 const PetscErrorCode ierr = VecAXPBY(vector, a, 0.0, v.vector);
949 AssertThrow(ierr == 0, ExcPETScError(ierr));
950 }
951
952
953
954 void
955 VectorBase::write_ascii(const PetscViewerFormat format)
956 {
957 // TODO[TH]:assert(is_compressed())
958 MPI_Comm comm = PetscObjectComm((PetscObject)vector);
959
960 // Set options
961 PetscErrorCode ierr =
962 PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(comm), format);
963 AssertThrow(ierr == 0, ExcPETScError(ierr));
964
965 // Write to screen
966 ierr = VecView(vector, PETSC_VIEWER_STDOUT_(comm));
967 AssertThrow(ierr == 0, ExcPETScError(ierr));
968 ierr = PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(comm));
969 AssertThrow(ierr == 0, ExcPETScError(ierr));
970 }
971
972
973
974 void
975 VectorBase::print(std::ostream &out,
976 const unsigned int precision,
977 const bool scientific,
978 const bool across) const
979 {
980 AssertThrow(out.fail() == false, ExcIO());
981
982 // get a representation of the vector and
983 // loop over all the elements
984 const PetscScalar *val;
985 PetscErrorCode ierr = VecGetArrayRead(vector, &val);
986
987 AssertThrow(ierr == 0, ExcPETScError(ierr));
988
989 // save the state of out stream
990 const std::ios::fmtflags old_flags = out.flags();
991 const unsigned int old_precision = out.precision(precision);
992
993 out.precision(precision);
994 if (scientific)
995 out.setf(std::ios::scientific, std::ios::floatfield);
996 else
997 out.setf(std::ios::fixed, std::ios::floatfield);
998
999 if (across)
1000 for (size_type i = 0; i < locally_owned_size(); ++i)
1001 out << val[i] << ' ';
1002 else
1003 for (size_type i = 0; i < locally_owned_size(); ++i)
1004 out << val[i] << std::endl;
1005 out << std::endl;
1006
1007 // reset output format
1008 out.flags(old_flags);
1009 out.precision(old_precision);
1010
1011 // restore the representation of the
1012 // vector
1013 ierr = VecRestoreArrayRead(vector, &val);
1014 AssertThrow(ierr == 0, ExcPETScError(ierr));
1015
1016 AssertThrow(out.fail() == false, ExcIO());
1017 }
1018
1019
1020
1021 void
1023 {
1024 std::swap(this->vector, v.vector);
1025 std::swap(this->ghosted, v.ghosted);
1026 std::swap(this->last_action, v.last_action);
1027 // missing swap for IndexSet
1028 IndexSet t(this->ghost_indices);
1029 this->ghost_indices = v.ghost_indices;
1030 v.ghost_indices = t;
1031 }
1032
1033
1034 VectorBase::operator const Vec &() const
1035 {
1036 return vector;
1037 }
1038
1039
1040 Vec &
1042 {
1043 return vector;
1044 }
1045
1046
1047 std::size_t
1049 {
1050 std::size_t mem = sizeof(Vec) + sizeof(last_action) +
1053
1054 // TH: I am relatively sure that PETSc is
1055 // storing the local data in a contiguous
1056 // block without indices:
1057 mem += locally_owned_size() * sizeof(PetscScalar);
1058 // assume that PETSc is storing one index
1059 // and one double per ghost element
1060 if (ghosted)
1061 mem += ghost_indices.n_elements() * (sizeof(PetscScalar) + sizeof(int));
1062
1063 // TODO[TH]: size of constant memory for PETSc?
1064 return mem;
1065 }
1066
1067
1068
1069 void
1071 const size_type *indices,
1072 const PetscScalar *values,
1073 const bool add_values)
1074 {
1078 internal::VectorReference::ExcWrongMode(action, last_action));
1080
1081 std::vector<PetscInt> petsc_indices(n_elements);
1082 for (size_type i = 0; i < n_elements; ++i)
1083 {
1084 const auto petsc_index = static_cast<PetscInt>(indices[i]);
1085 AssertIntegerConversion(petsc_index, indices[i]);
1086 petsc_indices[i] = petsc_index;
1087 }
1088
1089 const InsertMode mode = (add_values ? ADD_VALUES : INSERT_VALUES);
1090 const PetscErrorCode ierr = VecSetValues(
1091 vector, petsc_indices.size(), petsc_indices.data(), values, mode);
1092 AssertThrow(ierr == 0, ExcPETScError(ierr));
1093
1094 last_action = action;
1095 }
1096
1097} // namespace PETScWrappers
1098
1100
1101#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:35
constexpr bool running_in_debug_mode()
Definition config.h:73
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#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:928
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)