Reference documentation for deal.II version GIT b8135fa6eb 2023-03-25 15:55: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 - 2022 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  /* TODO GHOSTED */
153  const PetscErrorCode ierr =
154  PetscObjectReference(reinterpret_cast<PetscObject>(vector));
155  AssertNothrow(ierr == 0, ExcPETScError(ierr));
156  (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  {
173  /* TODO GHOSTED */
175  ExcMessage("Cannot assign a new Vec"));
176  PetscErrorCode ierr =
177  PetscObjectReference(reinterpret_cast<PetscObject>(v));
178  AssertThrow(ierr == 0, ExcPETScError(ierr));
179  ierr = VecDestroy(&vector);
180  AssertThrow(ierr == 0, ExcPETScError(ierr));
181  vector = v;
182  }
183 
184  void
186  {
187  const PetscErrorCode ierr = VecDestroy(&vector);
188  AssertThrow(ierr == 0, ExcPETScError(ierr));
189 
190  ghosted = false;
193  }
194 
195 
196 
197  VectorBase &
199  {
200  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
201 
202  PetscErrorCode ierr = VecCopy(v, vector);
203  AssertThrow(ierr == 0, ExcPETScError(ierr));
204 
205  return *this;
206  }
207 
208 
209 
210  VectorBase &
211  VectorBase::operator=(const PetscScalar s)
212  {
213  AssertIsFinite(s);
214 
215  // TODO[TH]: assert(is_compressed())
216 
217  PetscErrorCode ierr = VecSet(vector, s);
218  AssertThrow(ierr == 0, ExcPETScError(ierr));
219 
220  if (has_ghost_elements())
221  {
222  Vec ghost = nullptr;
223  ierr = VecGhostGetLocalForm(vector, &ghost);
224  AssertThrow(ierr == 0, ExcPETScError(ierr));
225 
226  ierr = VecSet(ghost, s);
227  AssertThrow(ierr == 0, ExcPETScError(ierr));
228 
229  ierr = VecGhostRestoreLocalForm(vector, &ghost);
230  AssertThrow(ierr == 0, ExcPETScError(ierr));
231  }
232 
233  return *this;
234  }
235 
236 
237 
238  bool
240  {
241  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
242 
243  PetscBool flag;
244  const PetscErrorCode ierr = VecEqual(vector, v.vector, &flag);
245  AssertThrow(ierr == 0, ExcPETScError(ierr));
246 
247  return (flag == PETSC_TRUE);
248  }
249 
250 
251 
252  bool
254  {
255  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
256 
257  PetscBool flag;
258  const PetscErrorCode ierr = VecEqual(vector, v.vector, &flag);
259  AssertThrow(ierr == 0, ExcPETScError(ierr));
260 
261  return (flag == PETSC_FALSE);
262  }
263 
264 
265 
268  {
269  PetscInt sz;
270  const PetscErrorCode ierr = VecGetSize(vector, &sz);
271  AssertThrow(ierr == 0, ExcPETScError(ierr));
272 
273  return sz;
274  }
275 
276 
277 
280  {
281  PetscInt sz;
282  const PetscErrorCode ierr = VecGetLocalSize(vector, &sz);
283  AssertThrow(ierr == 0, ExcPETScError(ierr));
284 
285  return sz;
286  }
287 
288 
289 
292  {
293  PetscInt sz;
294  const PetscErrorCode ierr = VecGetLocalSize(vector, &sz);
295  AssertThrow(ierr == 0, ExcPETScError(ierr));
296 
297  return sz;
298  }
299 
300 
301 
302  std::pair<VectorBase::size_type, VectorBase::size_type>
304  {
305  PetscInt begin, end;
306  const PetscErrorCode ierr =
307  VecGetOwnershipRange(static_cast<const Vec &>(vector), &begin, &end);
308  AssertThrow(ierr == 0, ExcPETScError(ierr));
309 
310  return std::make_pair(begin, end);
311  }
312 
313 
314 
315  void
316  VectorBase::set(const std::vector<size_type> & indices,
317  const std::vector<PetscScalar> &values)
318  {
319  Assert(indices.size() == values.size(),
320  ExcMessage("Function called with arguments of different sizes"));
321  do_set_add_operation(indices.size(), indices.data(), values.data(), false);
322  }
323 
324 
325 
326  void
327  VectorBase::add(const std::vector<size_type> & indices,
328  const std::vector<PetscScalar> &values)
329  {
330  Assert(indices.size() == values.size(),
331  ExcMessage("Function called with arguments of different sizes"));
332  do_set_add_operation(indices.size(), indices.data(), values.data(), true);
333  }
334 
335 
336 
337  void
338  VectorBase::add(const std::vector<size_type> & indices,
339  const ::Vector<PetscScalar> &values)
340  {
341  Assert(indices.size() == values.size(),
342  ExcMessage("Function called with arguments of different sizes"));
343  do_set_add_operation(indices.size(), indices.data(), values.begin(), true);
344  }
345 
346 
347 
348  void
349  VectorBase::add(const size_type n_elements,
350  const size_type * indices,
351  const PetscScalar *values)
352  {
353  do_set_add_operation(n_elements, indices, values, true);
354  }
355 
356 
357 
358  PetscScalar
360  {
361  Assert(size() == vec.size(), ExcDimensionMismatch(size(), vec.size()));
362 
363  PetscScalar result;
364 
365  // For complex vectors, VecDot() computes
366  // val = (x,y) = y^H x,
367  // where y^H denotes the conjugate transpose of y.
368  // Note that this corresponds to the usual "mathematicians'"
369  // complex inner product where the SECOND argument gets the
370  // complex conjugate, which is also how we document this
371  // operation.
372  const PetscErrorCode ierr = VecDot(vec.vector, vector, &result);
373  AssertThrow(ierr == 0, ExcPETScError(ierr));
374 
375  return result;
376  }
377 
378 
379 
380  PetscScalar
381  VectorBase::add_and_dot(const PetscScalar a,
382  const VectorBase &V,
383  const VectorBase &W)
384  {
385  this->add(a, V);
386  return *this * W;
387  }
388 
389 
390 
391  void
393  {
394  {
395 # ifdef DEBUG
396 # ifdef DEAL_II_WITH_MPI
397  // Check that all processors agree that last_action is the same (or none!)
398 
399  int my_int_last_action = last_action;
400  int all_int_last_action;
401 
402  const int ierr = MPI_Allreduce(&my_int_last_action,
403  &all_int_last_action,
404  1,
405  MPI_INT,
406  MPI_BOR,
408  AssertThrowMPI(ierr);
409 
410  AssertThrow(all_int_last_action !=
412  ExcMessage("Error: not all processors agree on the last "
413  "VectorOperation before this compress() call."));
414 # endif
415 # endif
416  }
417 
418  AssertThrow(
420  ExcMessage(
421  "Missing compress() or calling with wrong VectorOperation argument."));
422 
423  // note that one may think that
424  // we only need to do something
425  // if in fact the state is
426  // anything but
427  // last_action::unknown. but
428  // that's not true: one
429  // frequently gets into
430  // situations where only one
431  // processor (or a subset of
432  // processors) actually writes
433  // something into a vector, but
434  // we still need to call
435  // VecAssemblyBegin/End on all
436  // processors.
437  PetscErrorCode ierr = VecAssemblyBegin(vector);
438  AssertThrow(ierr == 0, ExcPETScError(ierr));
439  ierr = VecAssemblyEnd(vector);
440  AssertThrow(ierr == 0, ExcPETScError(ierr));
441 
442  // reset the last action field to
443  // indicate that we're back to a
444  // pristine state
446  }
447 
448 
449 
452  {
453  const real_type d = l2_norm();
454  return d * d;
455  }
456 
457 
458 
459  PetscScalar
461  {
462  // We can only use our more efficient
463  // routine in the serial case.
464  if (dynamic_cast<const PETScWrappers::MPI::Vector *>(this) != nullptr)
465  {
466  PetscScalar sum;
467  const PetscErrorCode ierr = VecSum(vector, &sum);
468  AssertThrow(ierr == 0, ExcPETScError(ierr));
469  return sum / static_cast<PetscReal>(size());
470  }
471 
472  // get a representation of the vector and
473  // loop over all the elements
474  const PetscScalar *start_ptr;
475  PetscErrorCode ierr = VecGetArrayRead(vector, &start_ptr);
476  AssertThrow(ierr == 0, ExcPETScError(ierr));
477 
478  PetscScalar mean = 0;
479  {
480  PetscScalar sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0;
481 
482  // use modern processors better by
483  // allowing pipelined commands to be
484  // executed in parallel
485  const PetscScalar *ptr = start_ptr;
486  const PetscScalar *eptr = ptr + (size() / 4) * 4;
487  while (ptr != eptr)
488  {
489  sum0 += *ptr++;
490  sum1 += *ptr++;
491  sum2 += *ptr++;
492  sum3 += *ptr++;
493  }
494  // add up remaining elements
495  while (ptr != start_ptr + size())
496  sum0 += *ptr++;
497 
498  mean = (sum0 + sum1 + sum2 + sum3) / static_cast<PetscReal>(size());
499  }
500 
501  // restore the representation of the
502  // vector
503  ierr = VecRestoreArrayRead(vector, &start_ptr);
504  AssertThrow(ierr == 0, ExcPETScError(ierr));
505 
506  return mean;
507  }
508 
509 
512  {
513  real_type d;
514 
515  const PetscErrorCode ierr = VecNorm(vector, NORM_1, &d);
516  AssertThrow(ierr == 0, ExcPETScError(ierr));
517 
518  return d;
519  }
520 
521 
522 
525  {
526  real_type d;
527 
528  const PetscErrorCode ierr = VecNorm(vector, NORM_2, &d);
529  AssertThrow(ierr == 0, ExcPETScError(ierr));
530 
531  return d;
532  }
533 
534 
535 
538  {
539  // get a representation of the vector and
540  // loop over all the elements
541  const PetscScalar *start_ptr;
542  PetscErrorCode ierr = VecGetArrayRead(vector, &start_ptr);
543  AssertThrow(ierr == 0, ExcPETScError(ierr));
544 
545  real_type norm = 0;
546  {
547  real_type sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0;
548 
549  // use modern processors better by
550  // allowing pipelined commands to be
551  // executed in parallel
552  const PetscScalar *ptr = start_ptr;
553  const PetscScalar *eptr = ptr + (size() / 4) * 4;
554  while (ptr != eptr)
555  {
556  sum0 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
557  sum1 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
558  sum2 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
559  sum3 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
560  }
561  // add up remaining elements
562  while (ptr != start_ptr + size())
563  sum0 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
564 
565  norm = std::pow(sum0 + sum1 + sum2 + sum3, 1. / p);
566  }
567 
568  // restore the representation of the
569  // vector
570  ierr = VecRestoreArrayRead(vector, &start_ptr);
571  AssertThrow(ierr == 0, ExcPETScError(ierr));
572 
573  return norm;
574  }
575 
576 
577 
580  {
581  real_type d;
582 
583  const PetscErrorCode ierr = VecNorm(vector, NORM_INFINITY, &d);
584  AssertThrow(ierr == 0, ExcPETScError(ierr));
585 
586  return d;
587  }
588 
589 
590 
593  {
594  PetscInt p;
595  real_type d;
596 
597  const PetscErrorCode ierr = VecMin(vector, &p, &d);
598  AssertThrow(ierr == 0, ExcPETScError(ierr));
599 
600  return d;
601  }
602 
603 
606  {
607  PetscInt p;
608  real_type d;
609 
610  const PetscErrorCode ierr = VecMax(vector, &p, &d);
611  AssertThrow(ierr == 0, ExcPETScError(ierr));
612 
613  return d;
614  }
615 
616 
617 
618  bool
620  {
621  // get a representation of the vector and
622  // loop over all the elements
623  const PetscScalar *start_ptr;
624  PetscErrorCode ierr = VecGetArrayRead(vector, &start_ptr);
625  AssertThrow(ierr == 0, ExcPETScError(ierr));
626 
627  const PetscScalar *ptr = start_ptr,
628  *eptr = start_ptr + locally_owned_size();
629  bool flag = true;
630  while (ptr != eptr)
631  {
632  if (*ptr != value_type())
633  {
634  flag = false;
635  break;
636  }
637  ++ptr;
638  }
639 
640  // restore the representation of the
641  // vector
642  ierr = VecRestoreArrayRead(vector, &start_ptr);
643  AssertThrow(ierr == 0, ExcPETScError(ierr));
644 
645  return flag;
646  }
647 
648 
649  namespace internal
650  {
651  template <typename T>
652  bool
653  is_non_negative(const T &t)
654  {
655  return t >= 0;
656  }
657 
658 
659 
660  template <typename T>
661  bool
662  is_non_negative(const std::complex<T> &)
663  {
664  Assert(false,
665  ExcMessage("You can't ask a complex value "
666  "whether it is non-negative.")) return true;
667  }
668  } // namespace internal
669 
670 
671 
672  bool
674  {
675  // get a representation of the vector and
676  // loop over all the elements
677  const PetscScalar *start_ptr;
678  PetscErrorCode ierr = VecGetArrayRead(vector, &start_ptr);
679  AssertThrow(ierr == 0, ExcPETScError(ierr));
680 
681  const PetscScalar *ptr = start_ptr,
682  *eptr = start_ptr + locally_owned_size();
683  bool flag = true;
684  while (ptr != eptr)
685  {
686  if (!internal::is_non_negative(*ptr))
687  {
688  flag = false;
689  break;
690  }
691  ++ptr;
692  }
693 
694  // restore the representation of the
695  // vector
696  ierr = VecRestoreArrayRead(vector, &start_ptr);
697  AssertThrow(ierr == 0, ExcPETScError(ierr));
698 
699  return flag;
700  }
701 
702 
703 
704  VectorBase &
705  VectorBase::operator*=(const PetscScalar a)
706  {
708  AssertIsFinite(a);
709 
710  const PetscErrorCode ierr = VecScale(vector, a);
711  AssertThrow(ierr == 0, ExcPETScError(ierr));
712 
713  return *this;
714  }
715 
716 
717 
718  VectorBase &
719  VectorBase::operator/=(const PetscScalar a)
720  {
722  AssertIsFinite(a);
723 
724  const PetscScalar factor = 1. / a;
725  AssertIsFinite(factor);
726 
727  const PetscErrorCode ierr = VecScale(vector, factor);
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  VectorBase &
749  {
751  const PetscErrorCode ierr = VecAXPY(vector, -1, v);
752  AssertThrow(ierr == 0, ExcPETScError(ierr));
753 
754  return *this;
755  }
756 
757 
758 
759  void
760  VectorBase::add(const PetscScalar s)
761  {
763  AssertIsFinite(s);
764 
765  const PetscErrorCode ierr = VecShift(vector, s);
766  AssertThrow(ierr == 0, ExcPETScError(ierr));
767  }
768 
769 
770 
771  void
772  VectorBase::add(const PetscScalar a, const VectorBase &v)
773  {
775  AssertIsFinite(a);
776 
777  const PetscErrorCode ierr = VecAXPY(vector, a, v);
778  AssertThrow(ierr == 0, ExcPETScError(ierr));
779  }
780 
781 
782 
783  void
784  VectorBase::add(const PetscScalar a,
785  const VectorBase &v,
786  const PetscScalar b,
787  const VectorBase &w)
788  {
790  AssertIsFinite(a);
791  AssertIsFinite(b);
792 
793  const PetscScalar weights[2] = {a, b};
794  Vec addends[2] = {v.vector, w.vector};
795 
796  const PetscErrorCode ierr = VecMAXPY(vector, 2, weights, addends);
797  AssertThrow(ierr == 0, ExcPETScError(ierr));
798  }
799 
800 
801 
802  void
803  VectorBase::sadd(const PetscScalar s, const VectorBase &v)
804  {
806  AssertIsFinite(s);
807 
808  const PetscErrorCode ierr = VecAYPX(vector, s, v);
809  AssertThrow(ierr == 0, ExcPETScError(ierr));
810  }
811 
812 
813 
814  void
815  VectorBase::sadd(const PetscScalar s,
816  const PetscScalar a,
817  const VectorBase &v)
818  {
820  AssertIsFinite(s);
821  AssertIsFinite(a);
822 
823  // there is nothing like a AXPAY
824  // operation in Petsc, so do it in two
825  // steps
826  *this *= s;
827  add(a, v);
828  }
829 
830 
831 
832  void
834  {
836  const PetscErrorCode ierr = VecPointwiseMult(vector, factors, vector);
837  AssertThrow(ierr == 0, ExcPETScError(ierr));
838  }
839 
840 
841 
842  void
843  VectorBase::equ(const PetscScalar a, const VectorBase &v)
844  {
846  AssertIsFinite(a);
847 
848  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
849 
850  // there is no simple operation for this
851  // in PETSc. there are multiple ways to
852  // emulate it, we choose this one:
853  const PetscErrorCode ierr = VecCopy(v.vector, vector);
854  AssertThrow(ierr == 0, ExcPETScError(ierr));
855 
856  *this *= a;
857  }
858 
859 
860 
861  void
862  VectorBase::write_ascii(const PetscViewerFormat format)
863  {
864  // TODO[TH]:assert(is_compressed())
865  MPI_Comm comm = PetscObjectComm((PetscObject)vector);
866 
867  // Set options
868  PetscErrorCode ierr =
869  PetscViewerSetFormat(PETSC_VIEWER_STDOUT_(comm), format);
870  AssertThrow(ierr == 0, ExcPETScError(ierr));
871 
872  // Write to screen
873  ierr = VecView(vector, PETSC_VIEWER_STDOUT_(comm));
874  AssertThrow(ierr == 0, ExcPETScError(ierr));
875  }
876 
877 
878 
879  void
880  VectorBase::print(std::ostream & out,
881  const unsigned int precision,
882  const bool scientific,
883  const bool across) const
884  {
885  AssertThrow(out.fail() == false, ExcIO());
886 
887  // get a representation of the vector and
888  // loop over all the elements
889  const PetscScalar *val;
890  PetscErrorCode ierr = VecGetArrayRead(vector, &val);
891 
892  AssertThrow(ierr == 0, ExcPETScError(ierr));
893 
894  // save the state of out stream
895  const std::ios::fmtflags old_flags = out.flags();
896  const unsigned int old_precision = out.precision(precision);
897 
898  out.precision(precision);
899  if (scientific)
900  out.setf(std::ios::scientific, std::ios::floatfield);
901  else
902  out.setf(std::ios::fixed, std::ios::floatfield);
903 
904  if (across)
905  for (size_type i = 0; i < locally_owned_size(); ++i)
906  out << val[i] << ' ';
907  else
908  for (size_type i = 0; i < locally_owned_size(); ++i)
909  out << val[i] << std::endl;
910  out << std::endl;
911 
912  // reset output format
913  out.flags(old_flags);
914  out.precision(old_precision);
915 
916  // restore the representation of the
917  // vector
918  ierr = VecRestoreArrayRead(vector, &val);
919  AssertThrow(ierr == 0, ExcPETScError(ierr));
920 
921  AssertThrow(out.fail() == false, ExcIO());
922  }
923 
924 
925 
926  void
928  {
929  const PetscErrorCode ierr = VecSwap(vector, v.vector);
930  AssertThrow(ierr == 0, ExcPETScError(ierr));
931  }
932 
933 
934  VectorBase::operator const Vec &() const
935  {
936  return vector;
937  }
938 
939 
940  Vec &
942  {
943  return vector;
944  }
945 
946 
947  std::size_t
949  {
950  std::size_t mem = sizeof(Vec) + sizeof(last_action) +
953 
954  // TH: I am relatively sure that PETSc is
955  // storing the local data in a contiguous
956  // block without indices:
957  mem += locally_owned_size() * sizeof(PetscScalar);
958  // assume that PETSc is storing one index
959  // and one double per ghost element
960  if (ghosted)
961  mem += ghost_indices.n_elements() * (sizeof(PetscScalar) + sizeof(int));
962 
963  // TODO[TH]: size of constant memory for PETSc?
964  return mem;
965  }
966 
967 
968 
969  void
971  const size_type * indices,
972  const PetscScalar *values,
973  const bool add_values)
974  {
975  VectorOperation::values action =
978  internal::VectorReference::ExcWrongMode(action, last_action));
980  // VecSetValues complains if we
981  // come with an empty
982  // vector. however, it is not a
983  // collective operation, so we
984  // can skip the call if necessary
985  // (unlike the above calls)
986  if (n_elements != 0)
987  {
988  const PetscInt *petsc_indices =
989  reinterpret_cast<const PetscInt *>(indices);
990 
991  const InsertMode mode = (add_values ? ADD_VALUES : INSERT_VALUES);
992  const PetscErrorCode ierr =
993  VecSetValues(vector, n_elements, petsc_indices, values, mode);
994  AssertThrow(ierr == 0, ExcPETScError(ierr));
995  }
996 
997  // set the mode here, independent of whether we have actually
998  // written elements or whether the list was empty
999  last_action = action;
1000  }
1001 
1002 } // namespace PETScWrappers
1003 
1005 
1006 #endif // DEAL_II_WITH_PETSC
size_type n_elements() const
Definition: index_set.h:1796
void clear()
Definition: index_set.h:1624
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)
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)
const MPI_Comm & get_mpi_communicator() const
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
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1586
#define AssertIsFinite(number)
Definition: exceptions.h:1854
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1886
#define AssertNothrow(cond, exc)
Definition: exceptions.h:1628
#define AssertIndexRange(index, range)
Definition: exceptions.h:1827
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcGhostsPresent()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1675
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< T >::value, std::size_t > memory_consumption(const T &t)
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)
const MPI_Comm & comm