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