Reference documentation for deal.II version Git 24041777c5 2021-12-02 20:43:25 -0700
\(\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 - 2021 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 
22 
23 # include <deal.II/lac/exceptions.h>
26 
27 # include <cmath>
28 
30 
31 namespace PETScWrappers
32 {
33  namespace internal
34  {
35 # ifndef DOXYGEN
36  VectorReference::operator PetscScalar() const
37  {
38  AssertIndexRange(index, vector.size());
39 
40  // The vector may have ghost entries. In that case, we first need to
41  // figure out which elements we own locally, then get a pointer to the
42  // elements that are stored here (both the ones we own as well as the
43  // ghost elements). In this array, the locally owned elements come first
44  // followed by the ghost elements whose position we can get from an
45  // index set.
46  if (vector.ghosted)
47  {
48  PetscInt begin, end;
49  PetscErrorCode ierr =
50  VecGetOwnershipRange(vector.vector, &begin, &end);
51  AssertThrow(ierr == 0, ExcPETScError(ierr));
52 
53  Vec locally_stored_elements = PETSC_NULL;
54  ierr = VecGhostGetLocalForm(vector.vector, &locally_stored_elements);
55  AssertThrow(ierr == 0, ExcPETScError(ierr));
56 
57  PetscInt lsize;
58  ierr = VecGetSize(locally_stored_elements, &lsize);
59  AssertThrow(ierr == 0, ExcPETScError(ierr));
60 
61  PetscScalar *ptr;
62  ierr = VecGetArray(locally_stored_elements, &ptr);
63  AssertThrow(ierr == 0, ExcPETScError(ierr));
64 
65  PetscScalar value;
66 
67  if (index >= static_cast<size_type>(begin) &&
68  index < static_cast<size_type>(end))
69  {
70  // local entry
71  value = *(ptr + index - begin);
72  }
73  else
74  {
75  // ghost entry
76  Assert(vector.ghost_indices.is_element(index),
77  ExcMessage(
78  "You are trying to access an element of a vector "
79  "that is neither a locally owned element nor a "
80  "ghost element of the vector."));
81  const size_type ghostidx =
82  vector.ghost_indices.index_within_set(index);
83 
84  AssertIndexRange(ghostidx + end - begin, lsize);
85  value = *(ptr + ghostidx + end - begin);
86  }
87 
88  ierr = VecRestoreArray(locally_stored_elements, &ptr);
89  AssertThrow(ierr == 0, ExcPETScError(ierr));
90 
91  ierr =
92  VecGhostRestoreLocalForm(vector.vector, &locally_stored_elements);
93  AssertThrow(ierr == 0, ExcPETScError(ierr));
94 
95  return value;
96  }
97 
98 
99  // first verify that the requested
100  // element is actually locally
101  // available
102  PetscInt begin, end;
103 
104  PetscErrorCode ierr = VecGetOwnershipRange(vector.vector, &begin, &end);
105  AssertThrow(ierr == 0, ExcPETScError(ierr));
106 
107  AssertThrow((index >= static_cast<size_type>(begin)) &&
108  (index < static_cast<size_type>(end)),
109  ExcAccessToNonlocalElement(index, begin, end - 1));
110 
111  PetscInt idx = index;
112  PetscScalar value;
113  ierr = VecGetValues(vector.vector, 1, &idx, &value);
114  AssertThrow(ierr == 0, ExcPETScError(ierr));
115 
116  return value;
117  }
118 # endif
119  } // namespace internal
120 
122  : vector(nullptr)
123  , ghosted(false)
124  , last_action(::VectorOperation::unknown)
125  , obtained_ownership(true)
126  {
128  ExcMessage("PETSc does not support multi-threaded access, set "
129  "the thread limit to 1 in MPI_InitFinalize()."));
130  }
131 
132 
133 
135  : Subscriptor()
136  , ghosted(v.ghosted)
138  , last_action(::VectorOperation::unknown)
139  , obtained_ownership(true)
140  {
142  ExcMessage("PETSc does not support multi-threaded access, set "
143  "the thread limit to 1 in MPI_InitFinalize()."));
144 
145  PetscErrorCode ierr = VecDuplicate(v.vector, &vector);
146  AssertThrow(ierr == 0, ExcPETScError(ierr));
147 
148  ierr = VecCopy(v.vector, vector);
149  AssertThrow(ierr == 0, ExcPETScError(ierr));
150  }
151 
152 
153 
155  : Subscriptor()
156  , vector(v)
157  , ghosted(false)
158  , last_action(::VectorOperation::unknown)
159  , obtained_ownership(false)
160  {
162  ExcMessage("PETSc does not support multi-threaded access, set "
163  "the thread limit to 1 in MPI_InitFinalize()."));
164  }
165 
166 
167 
169  {
170  if (obtained_ownership)
171  {
172  const PetscErrorCode ierr = VecDestroy(&vector);
173  AssertNothrow(ierr == 0, ExcPETScError(ierr));
174  (void)ierr;
175  }
176  }
177 
178 
179 
180  void
182  {
183  if (obtained_ownership)
184  {
185  const PetscErrorCode ierr = VecDestroy(&vector);
186  AssertThrow(ierr == 0, ExcPETScError(ierr));
187  }
188 
189  ghosted = false;
192  obtained_ownership = true;
193  }
194 
195 
196 
197  VectorBase &
198  VectorBase::operator=(const PetscScalar s)
199  {
200  AssertIsFinite(s);
201 
202  // TODO[TH]: assert(is_compressed())
203 
204  PetscErrorCode ierr = VecSet(vector, s);
205  AssertThrow(ierr == 0, ExcPETScError(ierr));
206 
207  if (has_ghost_elements())
208  {
209  Vec ghost = PETSC_NULL;
210  ierr = VecGhostGetLocalForm(vector, &ghost);
211  AssertThrow(ierr == 0, ExcPETScError(ierr));
212 
213  ierr = VecSet(ghost, s);
214  AssertThrow(ierr == 0, ExcPETScError(ierr));
215 
216  ierr = VecGhostRestoreLocalForm(vector, &ghost);
217  AssertThrow(ierr == 0, ExcPETScError(ierr));
218  }
219 
220  return *this;
221  }
222 
223 
224 
225  bool
227  {
228  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
229 
230  PetscBool flag;
231  const PetscErrorCode ierr = VecEqual(vector, v.vector, &flag);
232  AssertThrow(ierr == 0, ExcPETScError(ierr));
233 
234  return (flag == PETSC_TRUE);
235  }
236 
237 
238 
239  bool
241  {
242  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
243 
244  PetscBool flag;
245  const PetscErrorCode ierr = VecEqual(vector, v.vector, &flag);
246  AssertThrow(ierr == 0, ExcPETScError(ierr));
247 
248  return (flag == PETSC_FALSE);
249  }
250 
251 
252 
255  {
256  PetscInt sz;
257  const PetscErrorCode ierr = VecGetSize(vector, &sz);
258  AssertThrow(ierr == 0, ExcPETScError(ierr));
259 
260  return sz;
261  }
262 
263 
264 
267  {
268  PetscInt sz;
269  const PetscErrorCode ierr = VecGetLocalSize(vector, &sz);
270  AssertThrow(ierr == 0, ExcPETScError(ierr));
271 
272  return sz;
273  }
274 
275 
276 
279  {
280  PetscInt sz;
281  const PetscErrorCode ierr = VecGetLocalSize(vector, &sz);
282  AssertThrow(ierr == 0, ExcPETScError(ierr));
283 
284  return sz;
285  }
286 
287 
288 
289  std::pair<VectorBase::size_type, VectorBase::size_type>
291  {
292  PetscInt begin, end;
293  const PetscErrorCode ierr =
294  VecGetOwnershipRange(static_cast<const Vec &>(vector), &begin, &end);
295  AssertThrow(ierr == 0, ExcPETScError(ierr));
296 
297  return std::make_pair(begin, end);
298  }
299 
300 
301 
302  void
303  VectorBase::set(const std::vector<size_type> & indices,
304  const std::vector<PetscScalar> &values)
305  {
306  Assert(indices.size() == values.size(),
307  ExcMessage("Function called with arguments of different sizes"));
308  do_set_add_operation(indices.size(), indices.data(), values.data(), false);
309  }
310 
311 
312 
313  void
314  VectorBase::add(const std::vector<size_type> & indices,
315  const std::vector<PetscScalar> &values)
316  {
317  Assert(indices.size() == values.size(),
318  ExcMessage("Function called with arguments of different sizes"));
319  do_set_add_operation(indices.size(), indices.data(), values.data(), true);
320  }
321 
322 
323 
324  void
325  VectorBase::add(const std::vector<size_type> & indices,
326  const ::Vector<PetscScalar> &values)
327  {
328  Assert(indices.size() == values.size(),
329  ExcMessage("Function called with arguments of different sizes"));
330  do_set_add_operation(indices.size(), indices.data(), values.begin(), true);
331  }
332 
333 
334 
335  void
336  VectorBase::add(const size_type n_elements,
337  const size_type * indices,
338  const PetscScalar *values)
339  {
340  do_set_add_operation(n_elements, indices, values, true);
341  }
342 
343 
344 
345  PetscScalar
347  {
348  Assert(size() == vec.size(), ExcDimensionMismatch(size(), vec.size()));
349 
350  PetscScalar result;
351 
352  // For complex vectors, VecDot() computes
353  // val = (x,y) = y^H x,
354  // where y^H denotes the conjugate transpose of y.
355  // Note that this corresponds to the usual "mathematicians'"
356  // complex inner product where the SECOND argument gets the
357  // complex conjugate, which is also how we document this
358  // operation.
359  const PetscErrorCode ierr = VecDot(vec.vector, vector, &result);
360  AssertThrow(ierr == 0, ExcPETScError(ierr));
361 
362  return result;
363  }
364 
365 
366 
367  PetscScalar
368  VectorBase::add_and_dot(const PetscScalar a,
369  const VectorBase &V,
370  const VectorBase &W)
371  {
372  this->add(a, V);
373  return *this * W;
374  }
375 
376 
377 
378  void
380  {
381  {
382 # ifdef DEBUG
383 # ifdef DEAL_II_WITH_MPI
384  // Check that all processors agree that last_action is the same (or none!)
385 
386  int my_int_last_action = last_action;
387  int all_int_last_action;
388 
389  const int ierr = MPI_Allreduce(&my_int_last_action,
390  &all_int_last_action,
391  1,
392  MPI_INT,
393  MPI_BOR,
395  AssertThrowMPI(ierr);
396 
397  AssertThrow(all_int_last_action != (::VectorOperation::add |
399  ExcMessage("Error: not all processors agree on the last "
400  "VectorOperation before this compress() call."));
401 # endif
402 # endif
403  }
404 
405  AssertThrow(
407  last_action == operation,
408  ExcMessage(
409  "Missing compress() or calling with wrong VectorOperation argument."));
410 
411  // note that one may think that
412  // we only need to do something
413  // if in fact the state is
414  // anything but
415  // last_action::unknown. but
416  // that's not true: one
417  // frequently gets into
418  // situations where only one
419  // processor (or a subset of
420  // processors) actually writes
421  // something into a vector, but
422  // we still need to call
423  // VecAssemblyBegin/End on all
424  // processors.
425  PetscErrorCode ierr = VecAssemblyBegin(vector);
426  AssertThrow(ierr == 0, ExcPETScError(ierr));
427  ierr = VecAssemblyEnd(vector);
428  AssertThrow(ierr == 0, ExcPETScError(ierr));
429 
430  // reset the last action field to
431  // indicate that we're back to a
432  // pristine state
434  }
435 
436 
437 
440  {
441  const real_type d = l2_norm();
442  return d * d;
443  }
444 
445 
446 
447  PetscScalar
449  {
450  // We can only use our more efficient
451  // routine in the serial case.
452  if (dynamic_cast<const PETScWrappers::MPI::Vector *>(this) != nullptr)
453  {
454  PetscScalar sum;
455  const PetscErrorCode ierr = VecSum(vector, &sum);
456  AssertThrow(ierr == 0, ExcPETScError(ierr));
457  return sum / static_cast<PetscReal>(size());
458  }
459 
460  // get a representation of the vector and
461  // loop over all the elements
462  PetscScalar * start_ptr;
463  PetscErrorCode ierr = VecGetArray(vector, &start_ptr);
464  AssertThrow(ierr == 0, ExcPETScError(ierr));
465 
466  PetscScalar mean = 0;
467  {
468  PetscScalar sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0;
469 
470  // use modern processors better by
471  // allowing pipelined commands to be
472  // executed in parallel
473  const PetscScalar *ptr = start_ptr;
474  const PetscScalar *eptr = ptr + (size() / 4) * 4;
475  while (ptr != eptr)
476  {
477  sum0 += *ptr++;
478  sum1 += *ptr++;
479  sum2 += *ptr++;
480  sum3 += *ptr++;
481  }
482  // add up remaining elements
483  while (ptr != start_ptr + size())
484  sum0 += *ptr++;
485 
486  mean = (sum0 + sum1 + sum2 + sum3) / static_cast<PetscReal>(size());
487  }
488 
489  // restore the representation of the
490  // vector
491  ierr = VecRestoreArray(vector, &start_ptr);
492  AssertThrow(ierr == 0, ExcPETScError(ierr));
493 
494  return mean;
495  }
496 
497 
500  {
501  real_type d;
502 
503  const PetscErrorCode ierr = VecNorm(vector, NORM_1, &d);
504  AssertThrow(ierr == 0, ExcPETScError(ierr));
505 
506  return d;
507  }
508 
509 
510 
513  {
514  real_type d;
515 
516  const PetscErrorCode ierr = VecNorm(vector, NORM_2, &d);
517  AssertThrow(ierr == 0, ExcPETScError(ierr));
518 
519  return d;
520  }
521 
522 
523 
526  {
527  // get a representation of the vector and
528  // loop over all the elements
529  PetscScalar * start_ptr;
530  PetscErrorCode ierr = VecGetArray(vector, &start_ptr);
531  AssertThrow(ierr == 0, ExcPETScError(ierr));
532 
533  real_type norm = 0;
534  {
535  real_type sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0;
536 
537  // use modern processors better by
538  // allowing pipelined commands to be
539  // executed in parallel
540  const PetscScalar *ptr = start_ptr;
541  const PetscScalar *eptr = ptr + (size() / 4) * 4;
542  while (ptr != eptr)
543  {
544  sum0 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
545  sum1 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
546  sum2 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
547  sum3 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
548  }
549  // add up remaining elements
550  while (ptr != start_ptr + size())
551  sum0 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
552 
553  norm = std::pow(sum0 + sum1 + sum2 + sum3, 1. / p);
554  }
555 
556  // restore the representation of the
557  // vector
558  ierr = VecRestoreArray(vector, &start_ptr);
559  AssertThrow(ierr == 0, ExcPETScError(ierr));
560 
561  return norm;
562  }
563 
564 
565 
568  {
569  real_type d;
570 
571  const PetscErrorCode ierr = VecNorm(vector, NORM_INFINITY, &d);
572  AssertThrow(ierr == 0, ExcPETScError(ierr));
573 
574  return d;
575  }
576 
577 
578 
581  {
582  PetscInt p;
583  real_type d;
584 
585  const PetscErrorCode ierr = VecMin(vector, &p, &d);
586  AssertThrow(ierr == 0, ExcPETScError(ierr));
587 
588  return d;
589  }
590 
591 
594  {
595  PetscInt p;
596  real_type d;
597 
598  const PetscErrorCode ierr = VecMax(vector, &p, &d);
599  AssertThrow(ierr == 0, ExcPETScError(ierr));
600 
601  return d;
602  }
603 
604 
605 
606  bool
608  {
609  // get a representation of the vector and
610  // loop over all the elements
611  PetscScalar * start_ptr;
612  PetscErrorCode ierr = VecGetArray(vector, &start_ptr);
613  AssertThrow(ierr == 0, ExcPETScError(ierr));
614 
615  const PetscScalar *ptr = start_ptr,
616  *eptr = start_ptr + locally_owned_size();
617  bool flag = true;
618  while (ptr != eptr)
619  {
620  if (*ptr != value_type())
621  {
622  flag = false;
623  break;
624  }
625  ++ptr;
626  }
627 
628  // restore the representation of the
629  // vector
630  ierr = VecRestoreArray(vector, &start_ptr);
631  AssertThrow(ierr == 0, ExcPETScError(ierr));
632 
633  return flag;
634  }
635 
636 
637  namespace internal
638  {
639  template <typename T>
640  bool
641  is_non_negative(const T &t)
642  {
643  return t >= 0;
644  }
645 
646 
647 
648  template <typename T>
649  bool
650  is_non_negative(const std::complex<T> &)
651  {
652  Assert(false,
653  ExcMessage("You can't ask a complex value "
654  "whether it is non-negative.")) return true;
655  }
656  } // namespace internal
657 
658 
659 
660  bool
662  {
663  // get a representation of the vector and
664  // loop over all the elements
665  PetscScalar * start_ptr;
666  PetscErrorCode ierr = VecGetArray(vector, &start_ptr);
667  AssertThrow(ierr == 0, ExcPETScError(ierr));
668 
669  const PetscScalar *ptr = start_ptr,
670  *eptr = start_ptr + locally_owned_size();
671  bool flag = true;
672  while (ptr != eptr)
673  {
674  if (!internal::is_non_negative(*ptr))
675  {
676  flag = false;
677  break;
678  }
679  ++ptr;
680  }
681 
682  // restore the representation of the
683  // vector
684  ierr = VecRestoreArray(vector, &start_ptr);
685  AssertThrow(ierr == 0, ExcPETScError(ierr));
686 
687  return flag;
688  }
689 
690 
691 
692  VectorBase &
693  VectorBase::operator*=(const PetscScalar a)
694  {
696  AssertIsFinite(a);
697 
698  const PetscErrorCode ierr = VecScale(vector, a);
699  AssertThrow(ierr == 0, ExcPETScError(ierr));
700 
701  return *this;
702  }
703 
704 
705 
706  VectorBase &
707  VectorBase::operator/=(const PetscScalar a)
708  {
710  AssertIsFinite(a);
711 
712  const PetscScalar factor = 1. / a;
713  AssertIsFinite(factor);
714 
715  const PetscErrorCode ierr = VecScale(vector, factor);
716  AssertThrow(ierr == 0, ExcPETScError(ierr));
717 
718  return *this;
719  }
720 
721 
722 
723  VectorBase &
725  {
727  const PetscErrorCode ierr = VecAXPY(vector, 1, v);
728  AssertThrow(ierr == 0, ExcPETScError(ierr));
729 
730  return *this;
731  }
732 
733 
734 
735  VectorBase &
737  {
739  const PetscErrorCode ierr = VecAXPY(vector, -1, v);
740  AssertThrow(ierr == 0, ExcPETScError(ierr));
741 
742  return *this;
743  }
744 
745 
746 
747  void
748  VectorBase::add(const PetscScalar s)
749  {
751  AssertIsFinite(s);
752 
753  const PetscErrorCode ierr = VecShift(vector, s);
754  AssertThrow(ierr == 0, ExcPETScError(ierr));
755  }
756 
757 
758 
759  void
760  VectorBase::add(const PetscScalar a, const VectorBase &v)
761  {
763  AssertIsFinite(a);
764 
765  const PetscErrorCode ierr = VecAXPY(vector, a, v);
766  AssertThrow(ierr == 0, ExcPETScError(ierr));
767  }
768 
769 
770 
771  void
772  VectorBase::add(const PetscScalar a,
773  const VectorBase &v,
774  const PetscScalar b,
775  const VectorBase &w)
776  {
778  AssertIsFinite(a);
779  AssertIsFinite(b);
780 
781  const PetscScalar weights[2] = {a, b};
782  Vec addends[2] = {v.vector, w.vector};
783 
784  const PetscErrorCode ierr = VecMAXPY(vector, 2, weights, addends);
785  AssertThrow(ierr == 0, ExcPETScError(ierr));
786  }
787 
788 
789 
790  void
791  VectorBase::sadd(const PetscScalar s, const VectorBase &v)
792  {
794  AssertIsFinite(s);
795 
796  const PetscErrorCode ierr = VecAYPX(vector, s, v);
797  AssertThrow(ierr == 0, ExcPETScError(ierr));
798  }
799 
800 
801 
802  void
803  VectorBase::sadd(const PetscScalar s,
804  const PetscScalar a,
805  const VectorBase &v)
806  {
808  AssertIsFinite(s);
809  AssertIsFinite(a);
810 
811  // there is nothing like a AXPAY
812  // operation in Petsc, so do it in two
813  // steps
814  *this *= s;
815  add(a, v);
816  }
817 
818 
819 
820  void
822  {
824  const PetscErrorCode ierr = VecPointwiseMult(vector, factors, vector);
825  AssertThrow(ierr == 0, ExcPETScError(ierr));
826  }
827 
828 
829 
830  void
831  VectorBase::equ(const PetscScalar a, const VectorBase &v)
832  {
834  AssertIsFinite(a);
835 
836  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
837 
838  // there is no simple operation for this
839  // in PETSc. there are multiple ways to
840  // emulate it, we choose this one:
841  const PetscErrorCode ierr = VecCopy(v.vector, vector);
842  AssertThrow(ierr == 0, ExcPETScError(ierr));
843 
844  *this *= a;
845  }
846 
847 
848 
849  void
850  VectorBase::write_ascii(const PetscViewerFormat format)
851  {
852  // TODO[TH]:assert(is_compressed())
853 
854  // Set options
855  PetscErrorCode ierr =
856  PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD, format);
857  AssertThrow(ierr == 0, ExcPETScError(ierr));
858 
859  // Write to screen
860  ierr = VecView(vector, PETSC_VIEWER_STDOUT_WORLD);
861  AssertThrow(ierr == 0, ExcPETScError(ierr));
862  }
863 
864 
865 
866  void
867  VectorBase::print(std::ostream & out,
868  const unsigned int precision,
869  const bool scientific,
870  const bool across) const
871  {
872  AssertThrow(out, ExcIO());
873 
874  // get a representation of the vector and
875  // loop over all the elements
876  PetscScalar * val;
877  PetscErrorCode ierr = VecGetArray(vector, &val);
878 
879  AssertThrow(ierr == 0, ExcPETScError(ierr));
880 
881  // save the state of out stream
882  const std::ios::fmtflags old_flags = out.flags();
883  const unsigned int old_precision = out.precision(precision);
884 
885  out.precision(precision);
886  if (scientific)
887  out.setf(std::ios::scientific, std::ios::floatfield);
888  else
889  out.setf(std::ios::fixed, std::ios::floatfield);
890 
891  if (across)
892  for (size_type i = 0; i < locally_owned_size(); ++i)
893  out << val[i] << ' ';
894  else
895  for (size_type i = 0; i < locally_owned_size(); ++i)
896  out << val[i] << std::endl;
897  out << std::endl;
898 
899  // reset output format
900  out.flags(old_flags);
901  out.precision(old_precision);
902 
903  // restore the representation of the
904  // vector
905  ierr = VecRestoreArray(vector, &val);
906  AssertThrow(ierr == 0, ExcPETScError(ierr));
907 
908  AssertThrow(out, ExcIO());
909  }
910 
911 
912 
913  void
915  {
916  const PetscErrorCode ierr = VecSwap(vector, v.vector);
917  AssertThrow(ierr == 0, ExcPETScError(ierr));
918  }
919 
920 
921 
922  VectorBase::operator const Vec &() const
923  {
924  return vector;
925  }
926 
927 
928  std::size_t
930  {
931  std::size_t mem = sizeof(Vec) + sizeof(last_action) +
934 
935  // TH: I am relatively sure that PETSc is
936  // storing the local data in a contiguous
937  // block without indices:
938  mem += locally_owned_size() * sizeof(PetscScalar);
939  // assume that PETSc is storing one index
940  // and one double per ghost element
941  if (ghosted)
942  mem += ghost_indices.n_elements() * (sizeof(PetscScalar) + sizeof(int));
943 
944  // TODO[TH]: size of constant memory for PETSc?
945  return mem;
946  }
947 
948 
949 
950  void
952  const size_type * indices,
953  const PetscScalar *values,
954  const bool add_values)
955  {
957  (add_values ? ::VectorOperation::add :
959  Assert((last_action == action) ||
961  internal::VectorReference::ExcWrongMode(action, last_action));
963  // VecSetValues complains if we
964  // come with an empty
965  // vector. however, it is not a
966  // collective operation, so we
967  // can skip the call if necessary
968  // (unlike the above calls)
969  if (n_elements != 0)
970  {
971  const PetscInt *petsc_indices =
972  reinterpret_cast<const PetscInt *>(indices);
973 
974  const InsertMode mode = (add_values ? ADD_VALUES : INSERT_VALUES);
975  const PetscErrorCode ierr =
976  VecSetValues(vector, n_elements, petsc_indices, values, mode);
977  AssertThrow(ierr == 0, ExcPETScError(ierr));
978  }
979 
980  // set the mode here, independent of whether we have actually
981  // written elements or whether the list was empty
982  last_action = action;
983  }
984 
985 } // namespace PETScWrappers
986 
988 
989 #endif // DEAL_II_WITH_PETSC
#define AssertNothrow(cond, exc)
Definition: exceptions.h:1524
VectorBase & operator-=(const VectorBase &V)
types::global_dof_index size_type
Definition: cuda_kernels.h:45
static ::ExceptionBase & ExcIO()
void sadd(const PetscScalar s, const VectorBase &V)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1720
static const char V
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:472
size_type locally_owned_size() const
VectorBase & operator/=(const PetscScalar factor)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1571
void scale(const VectorBase &scaling_factors)
VectorBase & operator=(const VectorBase &)=delete
std::pair< size_type, size_type > local_range() const
real_type lp_norm(const real_type p) const
void clear()
Definition: index_set.h:1615
bool is_non_negative(const T &t)
static ::ExceptionBase & ExcMessage(std::string arg1)
void compress(const VectorOperation::values operation)
static const char T
T sum(const T &t, const MPI_Comm &mpi_communicator)
void do_set_add_operation(const size_type n_elements, const size_type *indices, const PetscScalar *values, const bool add_values)
VectorOperation::values last_action
#define Assert(cond, exc)
Definition: exceptions.h:1461
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static bool is_running_single_threaded()
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:404
PetscScalar add_and_dot(const PetscScalar a, const VectorBase &V, const VectorBase &W)
VectorType::value_type * end(VectorType &V)
std::size_t memory_consumption() const
void equ(const PetscScalar a, const VectorBase &V)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorBase & operator*=(const PetscScalar factor)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
static ::ExceptionBase & ExcGhostsPresent()
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorBase & operator+=(const VectorBase &V)
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1778
PetscScalar operator*(const VectorBase &vec) const
bool has_ghost_elements() const
PetscScalar mean_value() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:403
VectorType::value_type * begin(VectorType &V)
bool operator==(const VectorBase &v) const
void add(const std::vector< size_type > &indices, const std::vector< PetscScalar > &values)
void write_ascii(const PetscViewerFormat format=PETSC_VIEWER_DEFAULT)
virtual const MPI_Comm & get_mpi_communicator() const
virtual ~VectorBase() override
size_type n_elements() const
Definition: index_set.h:1837
#define AssertIsFinite(number)
Definition: exceptions.h:1746
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
bool operator!=(const VectorBase &v) const
void set(const std::vector< size_type > &indices, const std::vector< PetscScalar > &values)