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