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