Reference documentation for deal.II version 9.5.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// Copyright (C) 2004 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
17
18#ifdef DEAL_II_WITH_PETSC
19
21
25
26# include <cmath>
27
29
30namespace PETScWrappers
31{
32 namespace internal
33 {
34# ifndef DOXYGEN
35 VectorReference::operator PetscScalar() const
36 {
37 AssertIndexRange(index, vector.size());
38
39 // The vector may have ghost entries. In that case, we first need to
40 // figure out which elements we own locally, then get a pointer to the
41 // elements that are stored here (both the ones we own as well as the
42 // ghost elements). In this array, the locally owned elements come first
43 // followed by the ghost elements whose position we can get from an
44 // index set.
45 if (vector.ghosted)
46 {
47 PetscInt begin, end;
48 PetscErrorCode ierr =
49 VecGetOwnershipRange(vector.vector, &begin, &end);
50 AssertThrow(ierr == 0, ExcPETScError(ierr));
51
52 Vec locally_stored_elements = nullptr;
53 ierr = VecGhostGetLocalForm(vector.vector, &locally_stored_elements);
54 AssertThrow(ierr == 0, ExcPETScError(ierr));
55
56 PetscInt lsize;
57 ierr = VecGetSize(locally_stored_elements, &lsize);
58 AssertThrow(ierr == 0, ExcPETScError(ierr));
59
60 const PetscScalar *ptr;
61 ierr = VecGetArrayRead(locally_stored_elements, &ptr);
62 AssertThrow(ierr == 0, ExcPETScError(ierr));
63
64 PetscScalar value;
65
66 if (index >= static_cast<size_type>(begin) &&
67 index < static_cast<size_type>(end))
68 {
69 // local entry
70 value = *(ptr + index - begin);
71 }
72 else
73 {
74 // ghost entry
75 Assert(vector.ghost_indices.is_element(index),
77 "You are trying to access an element of a vector "
78 "that is neither a locally owned element nor a "
79 "ghost element of the vector."));
80 const size_type ghostidx =
81 vector.ghost_indices.index_within_set(index);
82
83 AssertIndexRange(ghostidx + end - begin, lsize);
84 value = *(ptr + ghostidx + end - begin);
85 }
86
87 ierr = VecRestoreArrayRead(locally_stored_elements, &ptr);
88 AssertThrow(ierr == 0, ExcPETScError(ierr));
89
90 ierr =
91 VecGhostRestoreLocalForm(vector.vector, &locally_stored_elements);
92 AssertThrow(ierr == 0, ExcPETScError(ierr));
93
94 return value;
95 }
96
97
98 // first verify that the requested
99 // element is actually locally
100 // available
101 PetscInt begin, end;
102
103 PetscErrorCode ierr = VecGetOwnershipRange(vector.vector, &begin, &end);
104 AssertThrow(ierr == 0, ExcPETScError(ierr));
105
106 AssertThrow((index >= static_cast<size_type>(begin)) &&
107 (index < static_cast<size_type>(end)),
108 ExcAccessToNonlocalElement(index, begin, end - 1));
109
110 const PetscScalar *ptr;
111 PetscScalar value;
112 ierr = VecGetArrayRead(vector.vector, &ptr);
113 AssertThrow(ierr == 0, ExcPETScError(ierr));
114 value = *(ptr + index - begin);
115 ierr = VecRestoreArrayRead(vector.vector, &ptr);
116 AssertThrow(ierr == 0, ExcPETScError(ierr));
117
118 return value;
119 }
120# endif
121 } // namespace internal
122
124 : vector(nullptr)
125 , ghosted(false)
126 , last_action(VectorOperation::unknown)
127 {}
128
129
130
132 : Subscriptor()
133 , ghosted(v.ghosted)
134 , ghost_indices(v.ghost_indices)
135 , last_action(VectorOperation::unknown)
136 {
137 PetscErrorCode ierr = VecDuplicate(v.vector, &vector);
138 AssertThrow(ierr == 0, ExcPETScError(ierr));
139
140 ierr = VecCopy(v.vector, vector);
141 AssertThrow(ierr == 0, ExcPETScError(ierr));
142 }
143
144
145
147 : Subscriptor()
148 , vector(v)
149 , ghosted(false)
150 , last_action(VectorOperation::unknown)
151 {
152 const PetscErrorCode ierr =
153 PetscObjectReference(reinterpret_cast<PetscObject>(vector));
154 AssertNothrow(ierr == 0, ExcPETScError(ierr));
155 (void)ierr;
157 }
158
159
160
162 {
163 const PetscErrorCode ierr = VecDestroy(&vector);
164 AssertNothrow(ierr == 0, ExcPETScError(ierr));
165 (void)ierr;
166 }
167
168
169
170 void
172 {
174 ExcMessage("Cannot assign a new Vec"));
175 PetscErrorCode ierr =
176 PetscObjectReference(reinterpret_cast<PetscObject>(v));
177 AssertThrow(ierr == 0, ExcPETScError(ierr));
178 ierr = VecDestroy(&vector);
179 AssertThrow(ierr == 0, ExcPETScError(ierr));
180 vector = v;
182 }
183
184
185
186 namespace
187 {
188 template <typename Iterator, typename OutType>
189 class ConvertingIterator
190 {
191 Iterator m_iterator;
192
193 public:
194 using difference_type =
195 typename std::iterator_traits<Iterator>::difference_type;
196 using value_type = OutType;
197 using pointer = OutType *;
198 using reference = OutType &;
199 using iterator_category = std::forward_iterator_tag;
200
201 ConvertingIterator(const Iterator &iterator)
202 : m_iterator(iterator)
203 {}
204
205 OutType
206 operator*() const
207 {
208 return static_cast<OutType>(std::real(*m_iterator));
209 }
210
211 ConvertingIterator &
212 operator++()
213 {
214 ++m_iterator;
215 return *this;
216 }
217
218 ConvertingIterator
219 operator++(int)
220 {
221 ConvertingIterator old = *this;
222 ++m_iterator;
223 return old;
224 }
225
226 bool
227 operator==(const ConvertingIterator &other) const
228 {
229 return this->m_iterator == other.m_iterator;
230 }
231
232 bool
233 operator!=(const ConvertingIterator &other) const
234 {
235 return this->m_iterator != other.m_iterator;
236 }
237 };
238 } // namespace
239
240
241
242 void
244 {
245 // Reset ghost data
246 ghosted = false;
248
249 // There's no API to infer ghost indices from a PETSc Vec which
250 // unfortunately doesn't allow integer entries. We use the
251 // "ConvertingIterator" class above to do an implicit conversion when
252 // sorting and adding ghost indices below.
253 PetscErrorCode ierr;
254 Vec ghosted_vec;
255 ierr = VecGhostGetLocalForm(vector, &ghosted_vec);
256 AssertThrow(ierr == 0, ExcPETScError(ierr));
257 if (ghosted_vec && ghosted_vec != vector)
258 {
259 Vec tvector;
260 PetscScalar *array;
261 PetscInt ghost_start_index, end_index, n_elements_stored_locally;
262
263 ierr = VecGhostRestoreLocalForm(vector, &ghosted_vec);
264 AssertThrow(ierr == 0, ExcPETScError(ierr));
265
266 ierr = VecGetOwnershipRange(vector, &ghost_start_index, &end_index);
267 AssertThrow(ierr == 0, ExcPETScError(ierr));
268 ierr = VecDuplicate(vector, &tvector);
269 AssertThrow(ierr == 0, ExcPETScError(ierr));
270 ierr = VecGetArray(tvector, &array);
271 AssertThrow(ierr == 0, ExcPETScError(ierr));
272
273 // Store the indices we care about in the vector, so that we can then
274 // exchange this information between processes. It is unfortunate that
275 // we have to store integers in floating point numbers. Let's at least
276 // make sure we do that in a way that ensures that when we get these
277 // numbers back as integers later on, we get the same thing.
278 for (PetscInt i = 0; i < end_index - ghost_start_index; i++)
279 {
280 Assert(static_cast<PetscInt>(static_cast<PetscScalar>(
281 ghost_start_index + i)) == (ghost_start_index + i),
283 array[i] = ghost_start_index + i;
284 }
285
286 ierr = VecRestoreArray(tvector, &array);
287 AssertThrow(ierr == 0, ExcPETScError(ierr));
288 ierr = VecGhostUpdateBegin(tvector, INSERT_VALUES, SCATTER_FORWARD);
289 AssertThrow(ierr == 0, ExcPETScError(ierr));
290 ierr = VecGhostUpdateEnd(tvector, INSERT_VALUES, SCATTER_FORWARD);
291 AssertThrow(ierr == 0, ExcPETScError(ierr));
292 ierr = VecGhostGetLocalForm(tvector, &ghosted_vec);
293 AssertThrow(ierr == 0, ExcPETScError(ierr));
294 ierr = VecGetLocalSize(ghosted_vec, &n_elements_stored_locally);
295 AssertThrow(ierr == 0, ExcPETScError(ierr));
296 ierr = VecGetArrayRead(ghosted_vec, (const PetscScalar **)&array);
297 AssertThrow(ierr == 0, ExcPETScError(ierr));
298
299 // Populate the 'ghosted' flag and the ghost_indices variable. The
300 // latter is an index set that is most efficiently filled by
301 // sorting the indices to add. At the same time, we don't want to
302 // sort the indices stored in a PETSc-owned array; so if the array
303 // is already sorted, pass that to the IndexSet variable, and if
304 // not then copy the indices, sort them, and then add those.
305 ghosted = true;
306 ghost_indices.set_size(this->size());
307
308 ConvertingIterator<PetscScalar *, types::global_dof_index> begin_ghosts(
309 &array[end_index - ghost_start_index]);
310 ConvertingIterator<PetscScalar *, types::global_dof_index> end_ghosts(
311 &array[n_elements_stored_locally]);
312 if (std::is_sorted(&array[end_index - ghost_start_index],
313 &array[n_elements_stored_locally],
314 [](PetscScalar left, PetscScalar right) {
315 return static_cast<PetscInt>(std::real(left)) <
316 static_cast<PetscInt>(std::real(right));
317 }))
318 {
319 ghost_indices.add_indices(begin_ghosts, end_ghosts);
320 }
321 else
322 {
323 std::vector<PetscInt> sorted_indices(begin_ghosts, end_ghosts);
324 std::sort(sorted_indices.begin(), sorted_indices.end());
325 ghost_indices.add_indices(sorted_indices.begin(),
326 sorted_indices.end());
327 }
329
330 ierr = VecRestoreArrayRead(ghosted_vec, (const PetscScalar **)&array);
331 AssertThrow(ierr == 0, ExcPETScError(ierr));
332 ierr = VecGhostRestoreLocalForm(tvector, &ghosted_vec);
333 AssertThrow(ierr == 0, ExcPETScError(ierr));
334 ierr = VecDestroy(&tvector);
335 AssertThrow(ierr == 0, ExcPETScError(ierr));
337 else
338 {
339 ierr = VecGhostRestoreLocalForm(vector, &ghosted_vec);
340 AssertThrow(ierr == 0, ExcPETScError(ierr));
341 }
342 }
343
344
345 void
347 {
348 const PetscErrorCode ierr = VecDestroy(&vector);
349 AssertThrow(ierr == 0, ExcPETScError(ierr));
350
351 ghosted = false;
354 }
355
356
357
358 VectorBase &
360 {
361 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
362
363 PetscErrorCode ierr = VecCopy(v, vector);
364 AssertThrow(ierr == 0, ExcPETScError(ierr));
365
366 return *this;
367 }
368
369
370
371 VectorBase &
372 VectorBase::operator=(const PetscScalar s)
373 {
375
376 // TODO[TH]: assert(is_compressed())
377
378 PetscErrorCode ierr = VecSet(vector, s);
379 AssertThrow(ierr == 0, ExcPETScError(ierr));
380
381 if (has_ghost_elements())
382 {
383 Vec ghost = nullptr;
384 ierr = VecGhostGetLocalForm(vector, &ghost);
385 AssertThrow(ierr == 0, ExcPETScError(ierr));
386
387 ierr = VecSet(ghost, s);
388 AssertThrow(ierr == 0, ExcPETScError(ierr));
389
390 ierr = VecGhostRestoreLocalForm(vector, &ghost);
391 AssertThrow(ierr == 0, ExcPETScError(ierr));
392 }
393
394 return *this;
395 }
396
397
398
399 bool
401 {
402 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
403
404 PetscBool flag;
405 const PetscErrorCode ierr = VecEqual(vector, v.vector, &flag);
406 AssertThrow(ierr == 0, ExcPETScError(ierr));
407
408 return (flag == PETSC_TRUE);
409 }
410
411
412
413 bool
415 {
416 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
417
418 PetscBool flag;
419 const PetscErrorCode ierr = VecEqual(vector, v.vector, &flag);
420 AssertThrow(ierr == 0, ExcPETScError(ierr));
421
422 return (flag == PETSC_FALSE);
423 }
424
425
426
429 {
430 PetscInt sz;
431 const PetscErrorCode ierr = VecGetSize(vector, &sz);
432 AssertThrow(ierr == 0, ExcPETScError(ierr));
433
434 return sz;
435 }
436
437
438
441 {
442 PetscInt sz;
443 const PetscErrorCode ierr = VecGetLocalSize(vector, &sz);
444 AssertThrow(ierr == 0, ExcPETScError(ierr));
445
446 return sz;
447 }
448
449
450
453 {
454 PetscInt sz;
455 const PetscErrorCode ierr = VecGetLocalSize(vector, &sz);
456 AssertThrow(ierr == 0, ExcPETScError(ierr));
457
458 return sz;
459 }
460
461
462
463 std::pair<VectorBase::size_type, VectorBase::size_type>
465 {
466 PetscInt begin, end;
467 const PetscErrorCode ierr =
468 VecGetOwnershipRange(static_cast<const Vec &>(vector), &begin, &end);
469 AssertThrow(ierr == 0, ExcPETScError(ierr));
470
471 return std::make_pair(begin, end);
472 }
473
474
475
476 void
477 VectorBase::set(const std::vector<size_type> & indices,
478 const std::vector<PetscScalar> &values)
479 {
480 Assert(indices.size() == values.size(),
481 ExcMessage("Function called with arguments of different sizes"));
482 do_set_add_operation(indices.size(), indices.data(), values.data(), false);
483 }
484
485
486
487 void
488 VectorBase::add(const std::vector<size_type> & indices,
489 const std::vector<PetscScalar> &values)
490 {
491 Assert(indices.size() == values.size(),
492 ExcMessage("Function called with arguments of different sizes"));
493 do_set_add_operation(indices.size(), indices.data(), values.data(), true);
494 }
495
496
497
498 void
499 VectorBase::add(const std::vector<size_type> & indices,
500 const ::Vector<PetscScalar> &values)
501 {
502 Assert(indices.size() == values.size(),
503 ExcMessage("Function called with arguments of different sizes"));
504 do_set_add_operation(indices.size(), indices.data(), values.begin(), true);
505 }
506
507
508
509 void
510 VectorBase::add(const size_type n_elements,
511 const size_type * indices,
512 const PetscScalar *values)
513 {
514 do_set_add_operation(n_elements, indices, values, true);
515 }
516
517
518
519 PetscScalar
521 {
522 Assert(size() == vec.size(), ExcDimensionMismatch(size(), vec.size()));
523
524 PetscScalar result;
525
526 // For complex vectors, VecDot() computes
527 // val = (x,y) = y^H x,
528 // where y^H denotes the conjugate transpose of y.
529 // Note that this corresponds to the usual "mathematicians'"
530 // complex inner product where the SECOND argument gets the
531 // complex conjugate, which is also how we document this
532 // operation.
533 const PetscErrorCode ierr = VecDot(vec.vector, vector, &result);
534 AssertThrow(ierr == 0, ExcPETScError(ierr));
535
536 return result;
537 }
538
539
540
541 PetscScalar
542 VectorBase::add_and_dot(const PetscScalar a,
543 const VectorBase &V,
544 const VectorBase &W)
545 {
546 this->add(a, V);
547 return *this * W;
548 }
549
550
551
552 void
554 {
555 {
556# ifdef DEBUG
557# ifdef DEAL_II_WITH_MPI
558 // Check that all processors agree that last_action is the same (or none!)
559
560 int my_int_last_action = last_action;
561 int all_int_last_action;
562
563 const int ierr = MPI_Allreduce(&my_int_last_action,
564 &all_int_last_action,
565 1,
566 MPI_INT,
567 MPI_BOR,
569 AssertThrowMPI(ierr);
570
571 AssertThrow(all_int_last_action !=
573 ExcMessage("Error: not all processors agree on the last "
574 "VectorOperation before this compress() call."));
575# endif
576# endif
577 }
578
582 "Missing compress() or calling with wrong VectorOperation argument."));
583
584 // note that one may think that
585 // we only need to do something
586 // if in fact the state is
587 // anything but
588 // last_action::unknown. but
589 // that's not true: one
590 // frequently gets into
591 // situations where only one
592 // processor (or a subset of
593 // processors) actually writes
594 // something into a vector, but
595 // we still need to call
596 // VecAssemblyBegin/End on all
597 // processors.
598 PetscErrorCode ierr = VecAssemblyBegin(vector);
599 AssertThrow(ierr == 0, ExcPETScError(ierr));
600 ierr = VecAssemblyEnd(vector);
601 AssertThrow(ierr == 0, ExcPETScError(ierr));
602
603 // reset the last action field to
604 // indicate that we're back to a
605 // pristine state
607 }
608
609
610
613 {
614 const real_type d = l2_norm();
615 return d * d;
616 }
617
618
619
620 PetscScalar
622 {
623 // We can only use our more efficient
624 // routine in the serial case.
625 if (dynamic_cast<const PETScWrappers::MPI::Vector *>(this) != nullptr)
626 {
627 PetscScalar sum;
628 const PetscErrorCode ierr = VecSum(vector, &sum);
629 AssertThrow(ierr == 0, ExcPETScError(ierr));
630 return sum / static_cast<PetscReal>(size());
631 }
632
633 // get a representation of the vector and
634 // loop over all the elements
635 const PetscScalar *start_ptr;
636 PetscErrorCode ierr = VecGetArrayRead(vector, &start_ptr);
637 AssertThrow(ierr == 0, ExcPETScError(ierr));
638
639 PetscScalar mean = 0;
640 {
641 PetscScalar sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0;
642
643 // use modern processors better by
644 // allowing pipelined commands to be
645 // executed in parallel
646 const PetscScalar *ptr = start_ptr;
647 const PetscScalar *eptr = ptr + (size() / 4) * 4;
648 while (ptr != eptr)
649 {
650 sum0 += *ptr++;
651 sum1 += *ptr++;
652 sum2 += *ptr++;
653 sum3 += *ptr++;
654 }
655 // add up remaining elements
656 while (ptr != start_ptr + size())
657 sum0 += *ptr++;
658
659 mean = (sum0 + sum1 + sum2 + sum3) / static_cast<PetscReal>(size());
660 }
661
662 // restore the representation of the
663 // vector
664 ierr = VecRestoreArrayRead(vector, &start_ptr);
665 AssertThrow(ierr == 0, ExcPETScError(ierr));
666
667 return mean;
668 }
669
670
673 {
674 real_type d;
675
676 const PetscErrorCode ierr = VecNorm(vector, NORM_1, &d);
677 AssertThrow(ierr == 0, ExcPETScError(ierr));
678
679 return d;
680 }
681
682
683
686 {
687 real_type d;
688
689 const PetscErrorCode ierr = VecNorm(vector, NORM_2, &d);
690 AssertThrow(ierr == 0, ExcPETScError(ierr));
691
692 return d;
693 }
694
695
696
699 {
700 // get a representation of the vector and
701 // loop over all the elements
702 const PetscScalar *start_ptr;
703 PetscErrorCode ierr = VecGetArrayRead(vector, &start_ptr);
704 AssertThrow(ierr == 0, ExcPETScError(ierr));
705
706 real_type norm = 0;
707 {
708 real_type sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0;
709
710 // use modern processors better by
711 // allowing pipelined commands to be
712 // executed in parallel
713 const PetscScalar *ptr = start_ptr;
714 const PetscScalar *eptr = ptr + (size() / 4) * 4;
715 while (ptr != eptr)
716 {
721 }
722 // add up remaining elements
723 while (ptr != start_ptr + size())
725
726 norm = std::pow(sum0 + sum1 + sum2 + sum3, 1. / p);
727 }
728
729 // restore the representation of the
730 // vector
731 ierr = VecRestoreArrayRead(vector, &start_ptr);
732 AssertThrow(ierr == 0, ExcPETScError(ierr));
733
734 return norm;
735 }
736
737
738
741 {
742 real_type d;
743
744 const PetscErrorCode ierr = VecNorm(vector, NORM_INFINITY, &d);
745 AssertThrow(ierr == 0, ExcPETScError(ierr));
746
747 return d;
748 }
749
750
751
754 {
755 PetscInt p;
756 real_type d;
757
758 const PetscErrorCode ierr = VecMin(vector, &p, &d);
759 AssertThrow(ierr == 0, ExcPETScError(ierr));
760
761 return d;
762 }
763
764
767 {
768 PetscInt p;
769 real_type d;
770
771 const PetscErrorCode ierr = VecMax(vector, &p, &d);
772 AssertThrow(ierr == 0, ExcPETScError(ierr));
773
774 return d;
775 }
776
777
778
779 bool
781 {
782 // get a representation of the vector and
783 // loop over all the elements
784 const PetscScalar *start_ptr;
785 PetscErrorCode ierr = VecGetArrayRead(vector, &start_ptr);
786 AssertThrow(ierr == 0, ExcPETScError(ierr));
787
788 const PetscScalar *ptr = start_ptr,
789 *eptr = start_ptr + locally_owned_size();
790 bool flag = true;
791 while (ptr != eptr)
792 {
793 if (*ptr != value_type())
794 {
795 flag = false;
796 break;
797 }
798 ++ptr;
799 }
800
801 // restore the representation of the
802 // vector
803 ierr = VecRestoreArrayRead(vector, &start_ptr);
804 AssertThrow(ierr == 0, ExcPETScError(ierr));
805
806 return flag;
807 }
808
809
810 namespace internal
811 {
812 template <typename T>
813 bool
814 is_non_negative(const T &t)
815 {
816 return t >= 0;
817 }
818
819
820
821 template <typename T>
822 bool
823 is_non_negative(const std::complex<T> &)
824 {
825 Assert(false,
826 ExcMessage("You can't ask a complex value "
827 "whether it is non-negative.")) return true;
828 }
829 } // namespace internal
830
831
832
833 bool
835 {
836 // get a representation of the vector and
837 // loop over all the elements
838 const PetscScalar *start_ptr;
839 PetscErrorCode ierr = VecGetArrayRead(vector, &start_ptr);
840 AssertThrow(ierr == 0, ExcPETScError(ierr));
841
842 const PetscScalar *ptr = start_ptr,
843 *eptr = start_ptr + locally_owned_size();
844 bool flag = true;
845 while (ptr != eptr)
846 {
847 if (!internal::is_non_negative(*ptr))
848 {
849 flag = false;
850 break;
851 }
852 ++ptr;
853 }
854
855 // restore the representation of the
856 // vector
857 ierr = VecRestoreArrayRead(vector, &start_ptr);
858 AssertThrow(ierr == 0, ExcPETScError(ierr));
859
860 return flag;
861 }
862
863
864
865 VectorBase &
866 VectorBase::operator*=(const PetscScalar a)
867 {
870
871 const PetscErrorCode ierr = VecScale(vector, a);
872 AssertThrow(ierr == 0, ExcPETScError(ierr));
873
874 return *this;
875 }
876
877
878
879 VectorBase &
880 VectorBase::operator/=(const PetscScalar a)
881 {
884
885 const PetscScalar factor = 1. / a;
886 AssertIsFinite(factor);
887
888 const PetscErrorCode ierr = VecScale(vector, factor);
889 AssertThrow(ierr == 0, ExcPETScError(ierr));
890
891 return *this;
892 }
893
894
895
896 VectorBase &
898 {
900 const PetscErrorCode ierr = VecAXPY(vector, 1, v);
901 AssertThrow(ierr == 0, ExcPETScError(ierr));
902
903 return *this;
904 }
905
906
907
908 VectorBase &
910 {
912 const PetscErrorCode ierr = VecAXPY(vector, -1, v);
913 AssertThrow(ierr == 0, ExcPETScError(ierr));
914
915 return *this;
916 }
917
918
919
920 void
921 VectorBase::add(const PetscScalar s)
922 {
925
926 const PetscErrorCode ierr = VecShift(vector, s);
927 AssertThrow(ierr == 0, ExcPETScError(ierr));
928 }
929
930
931
932 void
933 VectorBase::add(const PetscScalar a, const VectorBase &v)
934 {
937
938 const PetscErrorCode ierr = VecAXPY(vector, a, v);
939 AssertThrow(ierr == 0, ExcPETScError(ierr));
940 }
941
942
943
944 void
945 VectorBase::add(const PetscScalar a,
946 const VectorBase &v,
947 const PetscScalar b,
948 const VectorBase &w)
949 {
953
954 const PetscScalar weights[2] = {a, b};
955 Vec addends[2] = {v.vector, w.vector};
956
957 const PetscErrorCode ierr = VecMAXPY(vector, 2, weights, addends);
958 AssertThrow(ierr == 0, ExcPETScError(ierr));
959 }
960
961
962
963 void
964 VectorBase::sadd(const PetscScalar s, const VectorBase &v)
965 {
968
969 const PetscErrorCode ierr = VecAYPX(vector, s, v);
970 AssertThrow(ierr == 0, ExcPETScError(ierr));
971 }
972
973
974
975 void
976 VectorBase::sadd(const PetscScalar s,
977 const PetscScalar a,
978 const VectorBase &v)
979 {
983
984 // there is nothing like a AXPAY
985 // operation in PETSc, so do it in two
986 // steps
987 *this *= s;
988 add(a, v);
989 }
990
991
992
993 void
995 {
997 const PetscErrorCode ierr = VecPointwiseMult(vector, factors, vector);
998 AssertThrow(ierr == 0, ExcPETScError(ierr));
999 }
1000
1001
1002
1003 void
1004 VectorBase::equ(const PetscScalar a, const VectorBase &v)
1005 {
1007 AssertIsFinite(a);
1008
1009 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
1010
1011 const PetscErrorCode ierr = VecAXPBY(vector, a, 0.0, v.vector);
1012 AssertThrow(ierr == 0, ExcPETScError(ierr));
1013 }
1014
1015
1016
1017 void
1018 VectorBase::write_ascii(const PetscViewerFormat format)
1019 {
1020 // TODO[TH]:assert(is_compressed())
1021 MPI_Comm comm = PetscObjectComm((PetscObject)vector);
1022
1023 // Set options
1024 PetscErrorCode ierr =
1025 PetscViewerSetFormat(PETSC_VIEWER_STDOUT_(comm), format);
1026 AssertThrow(ierr == 0, ExcPETScError(ierr));
1027
1028 // Write to screen
1029 ierr = VecView(vector, PETSC_VIEWER_STDOUT_(comm));
1030 AssertThrow(ierr == 0, ExcPETScError(ierr));
1031 }
1032
1033
1034
1035 void
1036 VectorBase::print(std::ostream & out,
1037 const unsigned int precision,
1038 const bool scientific,
1039 const bool across) const
1040 {
1041 AssertThrow(out.fail() == false, ExcIO());
1042
1043 // get a representation of the vector and
1044 // loop over all the elements
1045 const PetscScalar *val;
1046 PetscErrorCode ierr = VecGetArrayRead(vector, &val);
1047
1048 AssertThrow(ierr == 0, ExcPETScError(ierr));
1049
1050 // save the state of out stream
1051 const std::ios::fmtflags old_flags = out.flags();
1052 const unsigned int old_precision = out.precision(precision);
1053
1054 out.precision(precision);
1055 if (scientific)
1056 out.setf(std::ios::scientific, std::ios::floatfield);
1057 else
1058 out.setf(std::ios::fixed, std::ios::floatfield);
1059
1060 if (across)
1061 for (size_type i = 0; i < locally_owned_size(); ++i)
1062 out << val[i] << ' ';
1063 else
1064 for (size_type i = 0; i < locally_owned_size(); ++i)
1065 out << val[i] << std::endl;
1066 out << std::endl;
1067
1068 // reset output format
1069 out.flags(old_flags);
1070 out.precision(old_precision);
1071
1072 // restore the representation of the
1073 // vector
1074 ierr = VecRestoreArrayRead(vector, &val);
1075 AssertThrow(ierr == 0, ExcPETScError(ierr));
1076
1077 AssertThrow(out.fail() == false, ExcIO());
1078 }
1079
1080
1081
1082 void
1084 {
1085 std::swap(this->vector, v.vector);
1086 std::swap(this->ghosted, v.ghosted);
1087 std::swap(this->last_action, v.last_action);
1088 // missing swap for IndexSet
1089 IndexSet t(this->ghost_indices);
1090 this->ghost_indices = v.ghost_indices;
1091 v.ghost_indices = t;
1092 }
1093
1094
1095 VectorBase::operator const Vec &() const
1096 {
1097 return vector;
1098 }
1099
1100
1101 Vec &
1103 {
1104 return vector;
1105 }
1106
1107
1108 std::size_t
1110 {
1111 std::size_t mem = sizeof(Vec) + sizeof(last_action) +
1114
1115 // TH: I am relatively sure that PETSc is
1116 // storing the local data in a contiguous
1117 // block without indices:
1118 mem += locally_owned_size() * sizeof(PetscScalar);
1119 // assume that PETSc is storing one index
1120 // and one double per ghost element
1121 if (ghosted)
1122 mem += ghost_indices.n_elements() * (sizeof(PetscScalar) + sizeof(int));
1123
1124 // TODO[TH]: size of constant memory for PETSc?
1125 return mem;
1126 }
1127
1128
1129
1130 void
1132 const size_type * indices,
1133 const PetscScalar *values,
1134 const bool add_values)
1135 {
1139 internal::VectorReference::ExcWrongMode(action, last_action));
1141 // VecSetValues complains if we
1142 // come with an empty
1143 // vector. however, it is not a
1144 // collective operation, so we
1145 // can skip the call if necessary
1146 // (unlike the above calls)
1147 if (n_elements != 0)
1148 {
1149 const PetscInt *petsc_indices =
1150 reinterpret_cast<const PetscInt *>(indices);
1151
1152 const InsertMode mode = (add_values ? ADD_VALUES : INSERT_VALUES);
1153 const PetscErrorCode ierr =
1154 VecSetValues(vector, n_elements, petsc_indices, values, mode);
1155 AssertThrow(ierr == 0, ExcPETScError(ierr));
1156 }
1157
1158 // set the mode here, independent of whether we have actually
1159 // written elements or whether the list was empty
1160 last_action = action;
1161 }
1162
1163} // namespace PETScWrappers
1164
1166
1167#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:1816
void set_size(const size_type size)
Definition index_set.h:1649
void clear()
Definition index_set.h:1637
void compress() const
Definition index_set.h:1669
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1727
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 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 &)
std::enable_if_t< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, 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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
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< T >::value, 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)
Iterator m_iterator
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
const MPI_Comm comm