Reference documentation for deal.II version GIT 574c7c8486 2023-09-22 10:20:02+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\}}\)
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 
22 # include <deal.II/lac/exceptions.h>
25 
26 # include <cmath>
27 
29 
30 namespace 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),
76  ExcMessage(
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;
156  this->determine_ghost_indices();
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;
181  this->determine_ghost_indices();
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>(std::real(static_cast<PetscScalar>(
281  ghost_start_index + i))) == (ghost_start_index + i),
282  ExcInternalError());
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));
336  }
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  {
374  AssertIsFinite(s);
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 
451  std::pair<VectorBase::size_type, VectorBase::size_type>
453  {
454  PetscInt begin, end;
455  const PetscErrorCode ierr =
456  VecGetOwnershipRange(static_cast<const Vec &>(vector), &begin, &end);
457  AssertThrow(ierr == 0, ExcPETScError(ierr));
458 
459  return std::make_pair(begin, end);
460  }
461 
462 
463 
464  void
465  VectorBase::set(const std::vector<size_type> &indices,
466  const std::vector<PetscScalar> &values)
467  {
468  Assert(indices.size() == values.size(),
469  ExcMessage("Function called with arguments of different sizes"));
470  do_set_add_operation(indices.size(), indices.data(), values.data(), false);
471  }
472 
473 
474 
475  void
476  VectorBase::add(const std::vector<size_type> &indices,
477  const std::vector<PetscScalar> &values)
478  {
479  Assert(indices.size() == values.size(),
480  ExcMessage("Function called with arguments of different sizes"));
481  do_set_add_operation(indices.size(), indices.data(), values.data(), true);
482  }
483 
484 
485 
486  void
487  VectorBase::add(const std::vector<size_type> &indices,
488  const ::Vector<PetscScalar> &values)
489  {
490  Assert(indices.size() == values.size(),
491  ExcMessage("Function called with arguments of different sizes"));
492  do_set_add_operation(indices.size(), indices.data(), values.begin(), true);
493  }
494 
495 
496 
497  void
498  VectorBase::add(const size_type n_elements,
499  const size_type *indices,
500  const PetscScalar *values)
501  {
502  do_set_add_operation(n_elements, indices, values, true);
503  }
504 
505 
506 
507  PetscScalar
509  {
510  Assert(size() == vec.size(), ExcDimensionMismatch(size(), vec.size()));
511 
512  PetscScalar result;
513 
514  // For complex vectors, VecDot() computes
515  // val = (x,y) = y^H x,
516  // where y^H denotes the conjugate transpose of y.
517  // Note that this corresponds to the usual "mathematicians'"
518  // complex inner product where the SECOND argument gets the
519  // complex conjugate, which is also how we document this
520  // operation.
521  const PetscErrorCode ierr = VecDot(vec.vector, vector, &result);
522  AssertThrow(ierr == 0, ExcPETScError(ierr));
523 
524  return result;
525  }
526 
527 
528 
529  PetscScalar
530  VectorBase::add_and_dot(const PetscScalar a,
531  const VectorBase &V,
532  const VectorBase &W)
533  {
534  this->add(a, V);
535  return *this * W;
536  }
537 
538 
539 
540  void
542  {
543  {
544 # ifdef DEBUG
545  // Check that all processors agree that last_action is the same (or none!)
546 
547  int my_int_last_action = last_action;
548  int all_int_last_action;
549 
550  const int ierr = MPI_Allreduce(&my_int_last_action,
551  &all_int_last_action,
552  1,
553  MPI_INT,
554  MPI_BOR,
556  AssertThrowMPI(ierr);
557 
558  AssertThrow(all_int_last_action !=
560  ExcMessage("Error: not all processors agree on the last "
561  "VectorOperation before this compress() call."));
562 # endif
563  }
564 
565  AssertThrow(
567  ExcMessage(
568  "Missing compress() or calling with wrong VectorOperation argument."));
569 
570  // note that one may think that
571  // we only need to do something
572  // if in fact the state is
573  // anything but
574  // last_action::unknown. but
575  // that's not true: one
576  // frequently gets into
577  // situations where only one
578  // processor (or a subset of
579  // processors) actually writes
580  // something into a vector, but
581  // we still need to call
582  // VecAssemblyBegin/End on all
583  // processors.
584  PetscErrorCode ierr = VecAssemblyBegin(vector);
585  AssertThrow(ierr == 0, ExcPETScError(ierr));
586  ierr = VecAssemblyEnd(vector);
587  AssertThrow(ierr == 0, ExcPETScError(ierr));
588 
589  // reset the last action field to
590  // indicate that we're back to a
591  // pristine state
593  }
594 
595 
596 
599  {
600  const real_type d = l2_norm();
601  return d * d;
602  }
603 
604 
605 
606  PetscScalar
608  {
609  // We can only use our more efficient
610  // routine in the serial case.
611  if (dynamic_cast<const PETScWrappers::MPI::Vector *>(this) != nullptr)
612  {
613  PetscScalar sum;
614  const PetscErrorCode ierr = VecSum(vector, &sum);
615  AssertThrow(ierr == 0, ExcPETScError(ierr));
616  return sum / static_cast<PetscReal>(size());
617  }
618 
619  // get a representation of the vector and
620  // loop over all the elements
621  const PetscScalar *start_ptr;
622  PetscErrorCode ierr = VecGetArrayRead(vector, &start_ptr);
623  AssertThrow(ierr == 0, ExcPETScError(ierr));
624 
625  PetscScalar mean = 0;
626  {
627  PetscScalar sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0;
628 
629  // use modern processors better by
630  // allowing pipelined commands to be
631  // executed in parallel
632  const PetscScalar *ptr = start_ptr;
633  const PetscScalar *eptr = ptr + (size() / 4) * 4;
634  while (ptr != eptr)
635  {
636  sum0 += *ptr++;
637  sum1 += *ptr++;
638  sum2 += *ptr++;
639  sum3 += *ptr++;
640  }
641  // add up remaining elements
642  while (ptr != start_ptr + size())
643  sum0 += *ptr++;
644 
645  mean = (sum0 + sum1 + sum2 + sum3) / static_cast<PetscReal>(size());
646  }
647 
648  // restore the representation of the
649  // vector
650  ierr = VecRestoreArrayRead(vector, &start_ptr);
651  AssertThrow(ierr == 0, ExcPETScError(ierr));
652 
653  return mean;
654  }
655 
656 
659  {
660  real_type d;
661 
662  const PetscErrorCode ierr = VecNorm(vector, NORM_1, &d);
663  AssertThrow(ierr == 0, ExcPETScError(ierr));
664 
665  return d;
666  }
667 
668 
669 
672  {
673  real_type d;
674 
675  const PetscErrorCode ierr = VecNorm(vector, NORM_2, &d);
676  AssertThrow(ierr == 0, ExcPETScError(ierr));
677 
678  return d;
679  }
680 
681 
682 
685  {
686  // get a representation of the vector and
687  // loop over all the elements
688  const PetscScalar *start_ptr;
689  PetscErrorCode ierr = VecGetArrayRead(vector, &start_ptr);
690  AssertThrow(ierr == 0, ExcPETScError(ierr));
691 
692  real_type norm = 0;
693  {
694  real_type sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0;
695 
696  // use modern processors better by
697  // allowing pipelined commands to be
698  // executed in parallel
699  const PetscScalar *ptr = start_ptr;
700  const PetscScalar *eptr = ptr + (size() / 4) * 4;
701  while (ptr != eptr)
702  {
703  sum0 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
704  sum1 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
705  sum2 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
706  sum3 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
707  }
708  // add up remaining elements
709  while (ptr != start_ptr + size())
710  sum0 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
711 
712  norm = std::pow(sum0 + sum1 + sum2 + sum3, 1. / p);
713  }
714 
715  // restore the representation of the
716  // vector
717  ierr = VecRestoreArrayRead(vector, &start_ptr);
718  AssertThrow(ierr == 0, ExcPETScError(ierr));
719 
720  return norm;
721  }
722 
723 
724 
727  {
728  real_type d;
729 
730  const PetscErrorCode ierr = VecNorm(vector, NORM_INFINITY, &d);
731  AssertThrow(ierr == 0, ExcPETScError(ierr));
732 
733  return d;
734  }
735 
736 
737 
738  bool
740  {
741  // get a representation of the vector and
742  // loop over all the elements
743  const PetscScalar *start_ptr;
744  PetscErrorCode ierr = VecGetArrayRead(vector, &start_ptr);
745  AssertThrow(ierr == 0, ExcPETScError(ierr));
746 
747  const PetscScalar *ptr = start_ptr,
748  *eptr = start_ptr + locally_owned_size();
749  bool flag = true;
750  while (ptr != eptr)
751  {
752  if (*ptr != value_type())
753  {
754  flag = false;
755  break;
756  }
757  ++ptr;
758  }
759 
760  // restore the representation of the
761  // vector
762  ierr = VecRestoreArrayRead(vector, &start_ptr);
763  AssertThrow(ierr == 0, ExcPETScError(ierr));
764 
765  return flag;
766  }
767 
768 
769  namespace internal
770  {
771  template <typename T>
772  bool
773  is_non_negative(const T &t)
774  {
775  return t >= 0;
776  }
777 
778 
779 
780  template <typename T>
781  bool
782  is_non_negative(const std::complex<T> &)
783  {
784  Assert(false,
785  ExcMessage("You can't ask a complex value "
786  "whether it is non-negative.")) return true;
787  }
788  } // namespace internal
789 
790 
791 
792  VectorBase &
793  VectorBase::operator*=(const PetscScalar a)
794  {
796  AssertIsFinite(a);
797 
798  const PetscErrorCode ierr = VecScale(vector, a);
799  AssertThrow(ierr == 0, ExcPETScError(ierr));
800 
801  return *this;
802  }
803 
804 
805 
806  VectorBase &
807  VectorBase::operator/=(const PetscScalar a)
808  {
810  AssertIsFinite(a);
811 
812  const PetscScalar factor = 1. / a;
813  AssertIsFinite(factor);
814 
815  const PetscErrorCode ierr = VecScale(vector, factor);
816  AssertThrow(ierr == 0, ExcPETScError(ierr));
817 
818  return *this;
819  }
820 
821 
822 
823  VectorBase &
825  {
827  const PetscErrorCode ierr = VecAXPY(vector, 1, v);
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  void
848  VectorBase::add(const PetscScalar s)
849  {
851  AssertIsFinite(s);
852 
853  const PetscErrorCode ierr = VecShift(vector, s);
854  AssertThrow(ierr == 0, ExcPETScError(ierr));
855  }
856 
857 
858 
859  void
860  VectorBase::add(const PetscScalar a, const VectorBase &v)
861  {
863  AssertIsFinite(a);
864 
865  const PetscErrorCode ierr = VecAXPY(vector, a, v);
866  AssertThrow(ierr == 0, ExcPETScError(ierr));
867  }
868 
869 
870 
871  void
872  VectorBase::add(const PetscScalar a,
873  const VectorBase &v,
874  const PetscScalar b,
875  const VectorBase &w)
876  {
878  AssertIsFinite(a);
879  AssertIsFinite(b);
880 
881  const PetscScalar weights[2] = {a, b};
882  Vec addends[2] = {v.vector, w.vector};
883 
884  const PetscErrorCode ierr = VecMAXPY(vector, 2, weights, addends);
885  AssertThrow(ierr == 0, ExcPETScError(ierr));
886  }
887 
888 
889 
890  void
891  VectorBase::sadd(const PetscScalar s, const VectorBase &v)
892  {
894  AssertIsFinite(s);
895 
896  const PetscErrorCode ierr = VecAYPX(vector, s, v);
897  AssertThrow(ierr == 0, ExcPETScError(ierr));
898  }
899 
900 
901 
902  void
903  VectorBase::sadd(const PetscScalar s,
904  const PetscScalar a,
905  const VectorBase &v)
906  {
908  AssertIsFinite(s);
909  AssertIsFinite(a);
910 
911  // there is nothing like a AXPAY
912  // operation in PETSc, so do it in two
913  // steps
914  *this *= s;
915  add(a, v);
916  }
917 
918 
919 
920  void
922  {
924  const PetscErrorCode ierr = VecPointwiseMult(vector, factors, vector);
925  AssertThrow(ierr == 0, ExcPETScError(ierr));
926  }
927 
928 
929 
930  void
931  VectorBase::equ(const PetscScalar a, const VectorBase &v)
932  {
934  AssertIsFinite(a);
935 
936  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
937 
938  const PetscErrorCode ierr = VecAXPBY(vector, a, 0.0, v.vector);
939  AssertThrow(ierr == 0, ExcPETScError(ierr));
940  }
941 
942 
943 
944  void
945  VectorBase::write_ascii(const PetscViewerFormat format)
946  {
947  // TODO[TH]:assert(is_compressed())
948  MPI_Comm comm = PetscObjectComm((PetscObject)vector);
949 
950  // Set options
951  PetscErrorCode ierr =
952  PetscViewerSetFormat(PETSC_VIEWER_STDOUT_(comm), format);
953  AssertThrow(ierr == 0, ExcPETScError(ierr));
954 
955  // Write to screen
956  ierr = VecView(vector, PETSC_VIEWER_STDOUT_(comm));
957  AssertThrow(ierr == 0, ExcPETScError(ierr));
958  }
959 
960 
961 
962  void
963  VectorBase::print(std::ostream &out,
964  const unsigned int precision,
965  const bool scientific,
966  const bool across) const
967  {
968  AssertThrow(out.fail() == false, ExcIO());
969 
970  // get a representation of the vector and
971  // loop over all the elements
972  const PetscScalar *val;
973  PetscErrorCode ierr = VecGetArrayRead(vector, &val);
974 
975  AssertThrow(ierr == 0, ExcPETScError(ierr));
976 
977  // save the state of out stream
978  const std::ios::fmtflags old_flags = out.flags();
979  const unsigned int old_precision = out.precision(precision);
980 
981  out.precision(precision);
982  if (scientific)
983  out.setf(std::ios::scientific, std::ios::floatfield);
984  else
985  out.setf(std::ios::fixed, std::ios::floatfield);
986 
987  if (across)
988  for (size_type i = 0; i < locally_owned_size(); ++i)
989  out << val[i] << ' ';
990  else
991  for (size_type i = 0; i < locally_owned_size(); ++i)
992  out << val[i] << std::endl;
993  out << std::endl;
994 
995  // reset output format
996  out.flags(old_flags);
997  out.precision(old_precision);
998 
999  // restore the representation of the
1000  // vector
1001  ierr = VecRestoreArrayRead(vector, &val);
1002  AssertThrow(ierr == 0, ExcPETScError(ierr));
1003 
1004  AssertThrow(out.fail() == false, ExcIO());
1005  }
1006 
1007 
1008 
1009  void
1011  {
1012  std::swap(this->vector, v.vector);
1013  std::swap(this->ghosted, v.ghosted);
1014  std::swap(this->last_action, v.last_action);
1015  // missing swap for IndexSet
1016  IndexSet t(this->ghost_indices);
1017  this->ghost_indices = v.ghost_indices;
1018  v.ghost_indices = t;
1019  }
1020 
1021 
1022  VectorBase::operator const Vec &() const
1023  {
1024  return vector;
1025  }
1026 
1027 
1028  Vec &
1030  {
1031  return vector;
1032  }
1033 
1034 
1035  std::size_t
1037  {
1038  std::size_t mem = sizeof(Vec) + sizeof(last_action) +
1041 
1042  // TH: I am relatively sure that PETSc is
1043  // storing the local data in a contiguous
1044  // block without indices:
1045  mem += locally_owned_size() * sizeof(PetscScalar);
1046  // assume that PETSc is storing one index
1047  // and one double per ghost element
1048  if (ghosted)
1049  mem += ghost_indices.n_elements() * (sizeof(PetscScalar) + sizeof(int));
1050 
1051  // TODO[TH]: size of constant memory for PETSc?
1052  return mem;
1053  }
1054 
1055 
1056 
1057  void
1059  const size_type *indices,
1060  const PetscScalar *values,
1061  const bool add_values)
1062  {
1063  VectorOperation::values action =
1066  internal::VectorReference::ExcWrongMode(action, last_action));
1068  // VecSetValues complains if we
1069  // come with an empty
1070  // vector. however, it is not a
1071  // collective operation, so we
1072  // can skip the call if necessary
1073  // (unlike the above calls)
1074  if (n_elements != 0)
1075  {
1076  const PetscInt *petsc_indices =
1077  reinterpret_cast<const PetscInt *>(indices);
1078 
1079  const InsertMode mode = (add_values ? ADD_VALUES : INSERT_VALUES);
1080  const PetscErrorCode ierr =
1081  VecSetValues(vector, n_elements, petsc_indices, values, mode);
1082  AssertThrow(ierr == 0, ExcPETScError(ierr));
1083  }
1084 
1085  // set the mode here, independent of whether we have actually
1086  // written elements or whether the list was empty
1087  last_action = action;
1088  }
1089 
1090 } // namespace PETScWrappers
1091 
1093 
1094 #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:1854
void set_size(const size_type size)
Definition: index_set.h:1684
void clear()
Definition: index_set.h:1672
void compress() const
Definition: index_set.h:1704
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1751
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)
virtual ~VectorBase() override
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 &)
types::global_dof_index size_type
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:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1616
#define AssertIsFinite(number)
Definition: exceptions.h:1884
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1916
#define AssertNothrow(cond, exc)
Definition: exceptions.h:1658
#define AssertIndexRange(index, range)
Definition: exceptions.h:1857
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcGhostsPresent()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
static const char T
static const char V
types::global_dof_index size_type
Definition: cuda_kernels.h:45
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:472
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
void swap(MemorySpaceData< T, MemorySpace > &u, MemorySpaceData< T, MemorySpace > &v)
bool is_non_negative(const T &t)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
Iterator m_iterator
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
const MPI_Comm comm