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