Reference documentation for deal.II version 8.5.1
petsc_vector_base.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2017 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 at
12 // the top level of the deal.II distribution.
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/lac/exceptions.h>
22 # include <deal.II/lac/petsc_compatibility.h>
23 # include <deal.II/lac/petsc_vector.h>
24 # include <deal.II/lac/petsc_parallel_vector.h>
25 # include <cmath>
26 # include <deal.II/base/multithread_info.h>
27 
28 DEAL_II_NAMESPACE_OPEN
29 
30 namespace PETScWrappers
31 {
32  namespace internal
33  {
34  VectorReference::operator PetscScalar () const
35  {
36  Assert (index < vector.size(),
37  ExcIndexRange (index, 0, vector.size()));
38 
39  // if the vector is local, then
40  // simply access the element we
41  // are interested in
42  if (dynamic_cast<const PETScWrappers::Vector *>(&vector) != 0)
43  {
44  PetscInt idx = index;
45  PetscScalar value;
46  PetscErrorCode ierr = VecGetValues(vector.vector, 1, &idx, &value);
47  AssertThrow (ierr == 0, ExcPETScError(ierr));
48  return value;
49  }
50  // else see if we are dealing
51  // with a parallel vector
52  else if (dynamic_cast<const PETScWrappers::MPI::Vector *>(&vector) != 0)
53  {
54  // there is the possibility
55  // that the vector has
56  // ghost elements. in that
57  // case, we first need to
58  // figure out which
59  // elements we own locally,
60  // then get a pointer to
61  // the elements that are
62  // stored here (both the
63  // ones we own as well as
64  // the ghost elements). in
65  // this array, the locally
66  // owned elements come
67  // first followed by the
68  // ghost elements whose
69  // position we can get from
70  // an index set
71  if (vector.ghosted)
72  {
73  PetscInt begin, end;
74  PetscErrorCode ierr = VecGetOwnershipRange (vector.vector, &begin,
75  &end);
76  AssertThrow (ierr == 0, ExcPETScError(ierr));
77 
78  Vec locally_stored_elements = PETSC_NULL;
79  ierr = VecGhostGetLocalForm(vector.vector, &locally_stored_elements);
80  AssertThrow (ierr == 0, ExcPETScError(ierr));
81 
82  PetscInt lsize;
83  ierr = VecGetSize(locally_stored_elements, &lsize);
84  AssertThrow (ierr == 0, ExcPETScError(ierr));
85 
86  PetscScalar *ptr;
87  ierr = VecGetArray(locally_stored_elements, &ptr);
88  AssertThrow (ierr == 0, ExcPETScError(ierr));
89 
90  PetscScalar value;
91 
92  if ( index>=static_cast<size_type>(begin)
93  && index<static_cast<size_type>(end) )
94  {
95  //local entry
96  value = *(ptr+index-begin);
97  }
98  else
99  {
100  //ghost entry
101  const size_type ghostidx
102  = vector.ghost_indices.index_within_set(index);
103 
104  Assert(ghostidx+end-begin<(size_type)lsize, ExcInternalError());
105  value = *(ptr+ghostidx+end-begin);
106  }
107 
108  ierr = VecRestoreArray(locally_stored_elements, &ptr);
109  AssertThrow (ierr == 0, ExcPETScError(ierr));
110 
111  ierr = VecGhostRestoreLocalForm(vector.vector, &locally_stored_elements);
112  AssertThrow (ierr == 0, ExcPETScError(ierr));
113 
114  return value;
115  }
116 
117 
118  // first verify that the requested
119  // element is actually locally
120  // available
121  PetscInt begin, end;
122 
123  PetscErrorCode ierr = VecGetOwnershipRange (vector.vector, &begin, &end);
124  AssertThrow (ierr == 0, ExcPETScError(ierr));
125 
126 
127 
128  AssertThrow ((index >= static_cast<size_type>(begin)) &&
129  (index < static_cast<size_type>(end)),
130  ExcAccessToNonlocalElement (index, begin, end-1));
131 
132  // old version which only work with
133  // VecGetArray()...
134  PetscInt idx = index;
135  PetscScalar value;
136  ierr = VecGetValues(vector.vector, 1, &idx, &value);
137  AssertThrow (ierr == 0, ExcPETScError(ierr));
138 
139  return value;
140  }
141  else
142  // what? what other kind of vector
143  // exists there?
144  Assert (false, ExcInternalError());
145 
146  return -1e20;
147  }
148  }
149 
151  :
152  vector (NULL),
153  ghosted(false),
154  last_action (::VectorOperation::unknown),
155  attained_ownership(true)
156  {
158  ExcMessage("PETSc does not support multi-threaded access, set "
159  "the thread limit to 1 in MPI_InitFinalize()."));
160  }
161 
162 
163 
165  :
166  Subscriptor (),
167  ghosted(v.ghosted),
168  ghost_indices(v.ghost_indices),
169  last_action (::VectorOperation::unknown),
170  attained_ownership(true)
171  {
173  ExcMessage("PETSc does not support multi-threaded access, set "
174  "the thread limit to 1 in MPI_InitFinalize()."));
175 
176  PetscErrorCode ierr = VecDuplicate (v.vector, &vector);
177  AssertThrow (ierr == 0, ExcPETScError(ierr));
178 
179  ierr = VecCopy (v.vector, vector);
180  AssertThrow (ierr == 0, ExcPETScError(ierr));
181  }
182 
183 
184 
185  VectorBase::VectorBase (const Vec &v)
186  :
187  Subscriptor (),
188  vector(v),
189  ghosted(false),
190  last_action (::VectorOperation::unknown),
191  attained_ownership(false)
192  {
194  ExcMessage("PETSc does not support multi-threaded access, set "
195  "the thread limit to 1 in MPI_InitFinalize()."));
196  }
197 
198 
199 
201  {
202  if (attained_ownership)
203  {
204 #if DEAL_II_PETSC_VERSION_LT(3,2,0)
205  const PetscErrorCode ierr = VecDestroy (vector);
206 #else
207  const PetscErrorCode ierr = VecDestroy (&vector);
208 #endif
209  AssertNothrow (ierr == 0, ExcPETScError(ierr));
210  (void) ierr;
211  }
212  }
213 
214 
215 
216  void
218  {
219  if (attained_ownership)
220  {
221 #if DEAL_II_PETSC_VERSION_LT(3,2,0)
222  const PetscErrorCode ierr = VecDestroy (vector);
223 #else
224  const PetscErrorCode ierr = VecDestroy (&vector);
225 #endif
226  AssertThrow (ierr == 0, ExcPETScError(ierr));
227  }
228 
229  ghosted = false;
230  ghost_indices.clear ();
232  attained_ownership = true;
233  }
234 
235 
236 
237  VectorBase &
238  VectorBase::operator = (const PetscScalar s)
239  {
240  AssertIsFinite(s);
241 
242  //TODO[TH]: assert(is_compressed())
243 
244  PetscErrorCode ierr = VecSet (vector, s);
245  AssertThrow (ierr == 0, ExcPETScError(ierr));
246 
247  if (has_ghost_elements())
248  {
249  Vec ghost = PETSC_NULL;
250  ierr = VecGhostGetLocalForm(vector, &ghost);
251  AssertThrow (ierr == 0, ExcPETScError(ierr));
252 
253  ierr = VecSet (ghost, s);
254  AssertThrow (ierr == 0, ExcPETScError(ierr));
255 
256  ierr = VecGhostRestoreLocalForm(vector, &ghost);
257  AssertThrow (ierr == 0, ExcPETScError(ierr));
258  }
259 
260  return *this;
261  }
262 
263 
264 
265  bool
267  {
268  Assert (size() == v.size(),
269  ExcDimensionMismatch(size(), v.size()));
270 
271  PetscBooleanType flag;
272  const PetscErrorCode ierr = VecEqual (vector, v.vector, &flag);
273  AssertThrow (ierr == 0, ExcPETScError(ierr));
274 
275  return (flag == PETSC_TRUE);
276  }
277 
278 
279 
280  bool
282  {
283  Assert (size() == v.size(),
284  ExcDimensionMismatch(size(), v.size()));
285 
286  PetscBooleanType flag;
287  const PetscErrorCode ierr = VecEqual (vector, v.vector, &flag);
288  AssertThrow (ierr == 0, ExcPETScError(ierr));
289 
290  return (flag == PETSC_FALSE);
291  }
292 
293 
294 
295  VectorBase::size_type
297  {
298  PetscInt sz;
299  const PetscErrorCode ierr = VecGetSize (vector, &sz);
300  AssertThrow (ierr == 0, ExcPETScError(ierr));
301 
302  return sz;
303  }
304 
305 
306 
307  VectorBase::size_type
309  {
310  PetscInt sz;
311  const PetscErrorCode ierr = VecGetLocalSize (vector, &sz);
312  AssertThrow (ierr == 0, ExcPETScError(ierr));
313 
314  return sz;
315  }
316 
317 
318 
319  std::pair<VectorBase::size_type, VectorBase::size_type>
321  {
322  PetscInt begin, end;
323  const PetscErrorCode ierr = VecGetOwnershipRange (static_cast<const Vec &>(vector),
324  &begin, &end);
325  AssertThrow (ierr == 0, ExcPETScError(ierr));
326 
327  return std::make_pair (begin, end);
328  }
329 
330 
331 
332  void
333  VectorBase::set (const std::vector<size_type> &indices,
334  const std::vector<PetscScalar> &values)
335  {
336  Assert (indices.size() == values.size(),
337  ExcMessage ("Function called with arguments of different sizes"));
338  do_set_add_operation(indices.size(), &indices[0], &values[0], false);
339  }
340 
341 
342 
343  void
344  VectorBase::add (const std::vector<size_type> &indices,
345  const std::vector<PetscScalar> &values)
346  {
347  Assert (indices.size() == values.size(),
348  ExcMessage ("Function called with arguments of different sizes"));
349  do_set_add_operation(indices.size(), &indices[0], &values[0], true);
350  }
351 
352 
353 
354  void
355  VectorBase::add (const std::vector<size_type> &indices,
356  const ::Vector<PetscScalar> &values)
357  {
358  Assert (indices.size() == values.size(),
359  ExcMessage ("Function called with arguments of different sizes"));
360  do_set_add_operation(indices.size(), &indices[0], values.begin(), true);
361  }
362 
363 
364 
365  void
366  VectorBase::add (const size_type n_elements,
367  const size_type *indices,
368  const PetscScalar *values)
369  {
370  do_set_add_operation(n_elements, indices, values, true);
371  }
372 
373 
374 
375  PetscScalar
377  {
378  Assert (size() == vec.size(),
379  ExcDimensionMismatch(size(), vec.size()));
380 
381  PetscScalar result;
382 
383  //For complex vectors, VecDot() computes
384  // val = (x,y) = y^H x,
385  //where y^H denotes the conjugate transpose of y.
386  //Note that this corresponds to the usual "mathematicians" complex inner product where the SECOND argument gets the complex conjugate.
387  const PetscErrorCode ierr = VecDot (vec.vector, vector, &result);
388  AssertThrow (ierr == 0, ExcPETScError(ierr));
389 
390  return result;
391  }
392 
393 
394 
395  PetscScalar
396  VectorBase::add_and_dot (const PetscScalar a,
397  const VectorBase &V,
398  const VectorBase &W)
399  {
400  this->add(a, V);
401  return *this * W;
402  }
403 
404 
405 
406  void
408  {
409  {
410 #ifdef DEBUG
411 #ifdef DEAL_II_WITH_MPI
412  // Check that all processors agree that last_action is the same (or none!)
413 
414  int my_int_last_action = last_action;
415  int all_int_last_action;
416 
417  const int ierr = MPI_Allreduce
418  (&my_int_last_action, &all_int_last_action, 1, MPI_INT, MPI_BOR,
420  AssertThrowMPI(ierr);
421 
422  AssertThrow(all_int_last_action !=
424  ExcMessage("Error: not all processors agree on the last "
425  "VectorOperation before this compress() call."));
426 #endif
427 #endif
428  }
429 
431  || last_action == operation,
432  ExcMessage("Missing compress() or calling with wrong VectorOperation argument."));
433 
434  // note that one may think that
435  // we only need to do something
436  // if in fact the state is
437  // anything but
438  // last_action::unknown. but
439  // that's not true: one
440  // frequently gets into
441  // situations where only one
442  // processor (or a subset of
443  // processors) actually writes
444  // something into a vector, but
445  // we still need to call
446  // VecAssemblyBegin/End on all
447  // processors.
448  PetscErrorCode ierr = VecAssemblyBegin(vector);
449  AssertThrow (ierr == 0, ExcPETScError(ierr));
450  ierr = VecAssemblyEnd(vector);
451  AssertThrow (ierr == 0, ExcPETScError(ierr));
452 
453  // reset the last action field to
454  // indicate that we're back to a
455  // pristine state
457  }
458 
459 
460 
461  VectorBase::real_type
463  {
464  const real_type d = l2_norm();
465  return d*d;
466  }
467 
468 
469 
470  PetscScalar
472  {
473  // We can only use our more efficient
474  // routine in the serial case.
475  if (dynamic_cast<const PETScWrappers::MPI::Vector *>(this) != 0)
476  {
477  PetscScalar sum;
478  const PetscErrorCode ierr = VecSum(vector, &sum);
479  AssertThrow (ierr == 0, ExcPETScError(ierr));
480  return sum/static_cast<PetscReal>(size());
481  }
482 
483  // get a representation of the vector and
484  // loop over all the elements
485  PetscScalar *start_ptr;
486  PetscErrorCode ierr = VecGetArray (vector, &start_ptr);
487  AssertThrow (ierr == 0, ExcPETScError(ierr));
488 
489  PetscScalar mean = 0;
490  {
491  PetscScalar sum0 = 0,
492  sum1 = 0,
493  sum2 = 0,
494  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 = VecRestoreArray (vector, &start_ptr);
518  AssertThrow (ierr == 0, ExcPETScError(ierr));
519 
520  return mean;
521  }
522 
523 
524  VectorBase::real_type
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 
537  VectorBase::real_type
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 
550  VectorBase::real_type
551  VectorBase::lp_norm (const real_type p) const
552  {
553  // get a representation of the vector and
554  // loop over all the elements
555  PetscScalar *start_ptr;
556  PetscErrorCode ierr = VecGetArray (vector, &start_ptr);
557  AssertThrow (ierr == 0, ExcPETScError(ierr));
558 
559  real_type norm = 0;
560  {
561  real_type sum0 = 0,
562  sum1 = 0,
563  sum2 = 0,
564  sum3 = 0;
565 
566  // use modern processors better by
567  // allowing pipelined commands to be
568  // executed in parallel
569  const PetscScalar *ptr = start_ptr;
570  const PetscScalar *eptr = ptr + (size()/4)*4;
571  while (ptr!=eptr)
572  {
573  sum0 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
574  sum1 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
575  sum2 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
576  sum3 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
577  }
578  // add up remaining elements
579  while (ptr != start_ptr+size())
580  sum0 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
581 
582  norm = std::pow(sum0+sum1+sum2+sum3,
583  1./p);
584  }
585 
586  // restore the representation of the
587  // vector
588  ierr = VecRestoreArray (vector, &start_ptr);
589  AssertThrow (ierr == 0, ExcPETScError(ierr));
590 
591  return norm;
592  }
593 
594 
595 
596  VectorBase::real_type
598  {
599  real_type d;
600 
601  const PetscErrorCode ierr = VecNorm (vector, NORM_INFINITY, &d);
602  AssertThrow (ierr == 0, ExcPETScError(ierr));
603 
604  return d;
605  }
606 
607 
608 
609  VectorBase::real_type
611  {
612  real_type d;
614  const PetscErrorCode ierr = VecNormalize (vector, &d);
615  AssertThrow (ierr == 0, ExcPETScError(ierr));
616 
617  return d;
618  }
619 
620 
621  VectorBase::real_type
623  {
624  PetscInt p;
625  real_type d;
626 
627  const PetscErrorCode ierr = VecMin (vector, &p, &d);
628  AssertThrow (ierr == 0, ExcPETScError(ierr));
629 
630  return d;
631  }
632 
633 
634  VectorBase::real_type
636  {
637  PetscInt p;
638  real_type d;
639 
640  const PetscErrorCode ierr = VecMax (vector, &p, &d);
641  AssertThrow (ierr == 0, ExcPETScError(ierr));
642 
643  return d;
644  }
645 
646 
647  VectorBase &
649  {
651 
652  const PetscErrorCode ierr = VecAbs (vector);
653  AssertThrow (ierr == 0, ExcPETScError(ierr));
654 
655  return *this;
656  }
657 
658 
659 
660  VectorBase &
662  {
664 
665  const PetscErrorCode ierr = VecConjugate (vector);
666  AssertThrow (ierr == 0, ExcPETScError(ierr));
667 
668  return *this;
669  }
670 
671 
672 
673  VectorBase &
675  {
677 
678  const PetscErrorCode ierr = VecPointwiseMult (vector,vector,vector);
679  AssertThrow (ierr == 0, ExcPETScError(ierr));
680 
681  return *this;
682  }
683 
684 
685  VectorBase &
687  {
689  const PetscErrorCode ierr = VecPointwiseMult (vector,vector,v);
690  AssertThrow (ierr == 0, ExcPETScError(ierr));
691 
692  return *this;
693  }
694 
695 
696  VectorBase &
698  const VectorBase &v)
699  {
701  const PetscErrorCode ierr = VecPointwiseMult (vector,u,v);
702  AssertThrow (ierr == 0, ExcPETScError(ierr));
703 
704  return *this;
705  }
706 
707  bool
709  {
710  // get a representation of the vector and
711  // loop over all the elements
712  PetscScalar *start_ptr;
713  PetscErrorCode ierr = VecGetArray (vector, &start_ptr);
714  AssertThrow (ierr == 0, ExcPETScError(ierr));
715 
716  const PetscScalar *ptr = start_ptr,
717  *eptr = start_ptr + local_size();
718  bool flag = true;
719  while (ptr != eptr)
720  {
721  if (*ptr != value_type())
722  {
723  flag = false;
724  break;
725  }
726  ++ptr;
727  }
728 
729  // restore the representation of the
730  // vector
731  ierr = VecRestoreArray (vector, &start_ptr);
732  AssertThrow (ierr == 0, ExcPETScError(ierr));
733 
734  return flag;
735  }
736 
737 
738  namespace internal
739  {
740  template <typename T>
741  bool is_non_negative (const T &t)
742  {
743  return t >= 0;
744  }
745 
746 
747 
748  template <typename T>
749  bool is_non_negative (const std::complex<T> &)
750  {
751  Assert (false,
752  ExcMessage ("You can't ask a complex value "
753  "whether it is non-negative."))
754  return true;
755  }
756  }
757 
758 
759 
760  bool
761  VectorBase::is_non_negative () const
762  {
763  // get a representation of the vector and
764  // loop over all the elements
765  PetscScalar *start_ptr;
766  PetscErrorCode ierr = VecGetArray (vector, &start_ptr);
767  AssertThrow (ierr == 0, ExcPETScError(ierr));
768 
769  const PetscScalar *ptr = start_ptr,
770  *eptr = start_ptr + local_size();
771  bool flag = true;
772  while (ptr != eptr)
773  {
774  if (! internal::is_non_negative(*ptr))
775  {
776  flag = false;
777  break;
778  }
779  ++ptr;
780  }
781 
782  // restore the representation of the
783  // vector
784  ierr = VecRestoreArray (vector, &start_ptr);
785  AssertThrow (ierr == 0, ExcPETScError(ierr));
786 
787  return flag;
788  }
789 
790 
791 
792  VectorBase &
793  VectorBase::operator *= (const PetscScalar a)
794  {
795  Assert (!has_ghost_elements(), ExcGhostsPresent());
796  AssertIsFinite(a);
797 
798  const PetscErrorCode ierr = VecScale (vector, a);
799  AssertThrow (ierr == 0, ExcPETScError(ierr));
800 
801  return *this;
802  }
803 
804 
805 
806  VectorBase &
807  VectorBase::operator /= (const PetscScalar a)
808  {
809  Assert (!has_ghost_elements(), ExcGhostsPresent());
810  AssertIsFinite(a);
811 
812  const PetscScalar factor = 1./a;
813  AssertIsFinite(factor);
814 
815  const PetscErrorCode ierr = VecScale (vector, factor);
816  AssertThrow (ierr == 0, ExcPETScError(ierr));
817 
818  return *this;
819  }
820 
821 
822 
823  VectorBase &
825  {
826  Assert (!has_ghost_elements(), ExcGhostsPresent());
827  const PetscErrorCode ierr = VecAXPY (vector, 1, v);
828  AssertThrow (ierr == 0, ExcPETScError(ierr));
829 
830  return *this;
831  }
832 
833 
834 
835  VectorBase &
837  {
838  Assert (!has_ghost_elements(), ExcGhostsPresent());
839  const PetscErrorCode ierr = VecAXPY (vector, -1, v);
840  AssertThrow (ierr == 0, ExcPETScError(ierr));
841 
842  return *this;
843  }
844 
845 
846 
847  void
848  VectorBase::add (const PetscScalar s)
849  {
850  Assert (!has_ghost_elements(), ExcGhostsPresent());
851  AssertIsFinite(s);
852 
853  const PetscErrorCode ierr = VecShift (vector, s);
854  AssertThrow (ierr == 0, ExcPETScError(ierr));
855  }
856 
857 
858 
859  void
861  {
862  *this += v;
863  }
864 
865 
866 
867  void
868  VectorBase::add (const PetscScalar a,
869  const VectorBase &v)
870  {
871  Assert (!has_ghost_elements(), ExcGhostsPresent());
872  AssertIsFinite(a);
873 
874  const PetscErrorCode ierr = VecAXPY (vector, a, v);
875  AssertThrow (ierr == 0, ExcPETScError(ierr));
876  }
877 
878 
879 
880  void
881  VectorBase::add (const PetscScalar a,
882  const VectorBase &v,
883  const PetscScalar b,
884  const VectorBase &w)
885  {
886  Assert (!has_ghost_elements(), ExcGhostsPresent());
887  AssertIsFinite(a);
888  AssertIsFinite(b);
889 
890  const PetscScalar weights[2] = {a,b};
891  Vec addends[2] = {v.vector, w.vector};
892 
893  const PetscErrorCode ierr = VecMAXPY (vector, 2, weights, addends);
894  AssertThrow (ierr == 0, ExcPETScError(ierr));
895  }
896 
897 
898 
899  void
900  VectorBase::sadd (const PetscScalar s,
901  const VectorBase &v)
902  {
903  Assert (!has_ghost_elements(), ExcGhostsPresent());
904  AssertIsFinite(s);
905 
906  const PetscErrorCode ierr = VecAYPX (vector, s, v);
907  AssertThrow (ierr == 0, ExcPETScError(ierr));
908  }
909 
910 
911 
912  void
913  VectorBase::sadd (const PetscScalar s,
914  const PetscScalar a,
915  const VectorBase &v)
916  {
917  Assert (!has_ghost_elements(), ExcGhostsPresent());
918  AssertIsFinite(s);
919  AssertIsFinite(a);
920 
921  // there is nothing like a AXPAY
922  // operation in Petsc, so do it in two
923  // steps
924  *this *= s;
925  add (a,v);
926  }
927 
928 
929 
930  void
931  VectorBase::sadd (const PetscScalar s,
932  const PetscScalar a,
933  const VectorBase &v,
934  const PetscScalar b,
935  const VectorBase &w)
936  {
937  Assert (!has_ghost_elements(), ExcGhostsPresent());
938  AssertIsFinite(s);
939  AssertIsFinite(a);
940  AssertIsFinite(b);
941 
942  // there is no operation like MAXPAY, so
943  // do it in two steps
944  *this *= s;
945 
946  const PetscScalar weights[2] = {a,b};
947  Vec addends[2] = {v.vector,w.vector};
948 
949  const PetscErrorCode ierr = VecMAXPY (vector, 2, weights, addends);
950  AssertThrow (ierr == 0, ExcPETScError(ierr));
951  }
952 
953 
954 
955  void
956  VectorBase::sadd (const PetscScalar s,
957  const PetscScalar a,
958  const VectorBase &v,
959  const PetscScalar b,
960  const VectorBase &w,
961  const PetscScalar c,
962  const VectorBase &x)
963  {
964  Assert (!has_ghost_elements(), ExcGhostsPresent());
965  AssertIsFinite(s);
966  AssertIsFinite(a);
967  AssertIsFinite(b);
968  AssertIsFinite(c);
969 
970  // there is no operation like MAXPAY, so
971  // do it in two steps
972  *this *= s;
973 
974  const PetscScalar weights[3] = {a,b,c};
975  Vec addends[3] = {v.vector, w.vector, x.vector};
976 
977  const PetscErrorCode ierr = VecMAXPY (vector, 3, weights, addends);
978  AssertThrow (ierr == 0, ExcPETScError(ierr));
979  }
980 
981 
982 
983  void
984  VectorBase::scale (const VectorBase &factors)
985  {
986  Assert (!has_ghost_elements(), ExcGhostsPresent());
987  const PetscErrorCode ierr
988  = VecPointwiseMult (vector, factors, vector);
989  AssertThrow (ierr == 0, ExcPETScError(ierr));
990  }
991 
992 
993 
994  void
995  VectorBase::equ (const PetscScalar a,
996  const VectorBase &v)
997  {
998  Assert (!has_ghost_elements(), ExcGhostsPresent());
999  AssertIsFinite(a);
1000 
1001  Assert (size() == v.size(),
1002  ExcDimensionMismatch (size(), v.size()));
1003 
1004  // there is no simple operation for this
1005  // in PETSc. there are multiple ways to
1006  // emulate it, we choose this one:
1007  const PetscErrorCode ierr = VecCopy (v.vector, vector);
1008  AssertThrow (ierr == 0, ExcPETScError(ierr));
1009 
1010  *this *= a;
1011  }
1012 
1013 
1014 
1015  void
1016  VectorBase::equ (const PetscScalar a,
1017  const VectorBase &v,
1018  const PetscScalar b,
1019  const VectorBase &w)
1020  {
1021  Assert (!has_ghost_elements(), ExcGhostsPresent());
1022  AssertIsFinite(a);
1023  AssertIsFinite(b);
1024 
1025  Assert (size() == v.size(),
1026  ExcDimensionMismatch (size(), v.size()));
1027 
1028  // there is no simple operation for this
1029  // in PETSc. there are multiple ways to
1030  // emulate it, we choose this one:
1031  const PetscErrorCode ierr = VecCopy (v.vector, vector);
1032  AssertThrow (ierr == 0, ExcPETScError(ierr));
1033 
1034  sadd (a, b, w);
1035  }
1036 
1037 
1038 
1039  void
1041  const VectorBase &b)
1042  {
1043  Assert (!has_ghost_elements(), ExcGhostsPresent());
1044  const PetscErrorCode ierr = VecPointwiseDivide (vector, a, b);
1045  AssertThrow (ierr == 0, ExcPETScError(ierr));
1046  }
1047 
1048 
1049 
1050  void
1051  VectorBase::write_ascii (const PetscViewerFormat format)
1052  {
1053  //TODO[TH]:assert(is_compressed())
1054 
1055  // Set options
1056  PetscErrorCode ierr = PetscViewerSetFormat (PETSC_VIEWER_STDOUT_WORLD,
1057  format);
1058  AssertThrow (ierr == 0, ExcPETScError(ierr));
1059 
1060  // Write to screen
1061  ierr = VecView (vector, PETSC_VIEWER_STDOUT_WORLD);
1062  AssertThrow (ierr == 0, ExcPETScError(ierr));
1063  }
1064 
1065 
1066 
1067  void
1068  VectorBase::print (std::ostream &out,
1069  const unsigned int precision,
1070  const bool scientific,
1071  const bool across) const
1072  {
1073  AssertThrow (out, ExcIO());
1074 
1075  // get a representation of the vector and
1076  // loop over all the elements
1077  PetscScalar *val;
1078  PetscErrorCode ierr = VecGetArray (vector, &val);
1079 
1080  AssertThrow (ierr == 0, ExcPETScError(ierr));
1081 
1082  // save the state of out stream
1083  const std::ios::fmtflags old_flags = out.flags();
1084  const unsigned int old_precision = out.precision (precision);
1085 
1086  out.precision (precision);
1087  if (scientific)
1088  out.setf (std::ios::scientific, std::ios::floatfield);
1089  else
1090  out.setf (std::ios::fixed, std::ios::floatfield);
1091 
1092  if (across)
1093  for (size_type i=0; i<local_size(); ++i)
1094  out << val[i] << ' ';
1095  else
1096  for (size_type i=0; i<local_size(); ++i)
1097  out << val[i] << std::endl;
1098  out << std::endl;
1099 
1100  // reset output format
1101  out.flags (old_flags);
1102  out.precision(old_precision);
1103 
1104  // restore the representation of the
1105  // vector
1106  ierr = VecRestoreArray (vector, &val);
1107  AssertThrow (ierr == 0, ExcPETScError(ierr));
1108 
1109  AssertThrow (out, ExcIO());
1110  }
1111 
1112 
1113 
1114  void
1116  {
1117  const PetscErrorCode ierr = VecSwap (vector, v.vector);
1118  AssertThrow (ierr == 0, ExcPETScError(ierr));
1119  }
1120 
1121 
1122 
1123  VectorBase::operator const Vec &() const
1124  {
1125  return vector;
1126  }
1127 
1128 
1129  std::size_t
1131  {
1132  std::size_t mem = sizeof(Vec)+sizeof(last_action)
1134  +MemoryConsumption::memory_consumption(ghost_indices);
1135 
1136  // TH: I am relatively sure that PETSc is
1137  // storing the local data in a contiguous
1138  // block without indices:
1139  mem += local_size()*sizeof(PetscScalar);
1140  // assume that PETSc is storing one index
1141  // and one double per ghost element
1142  if (ghosted)
1143  mem += ghost_indices.n_elements()*(sizeof(PetscScalar)+sizeof(int));
1144 
1145  //TODO[TH]: size of constant memory for PETSc?
1146  return mem;
1147  }
1148 
1149 
1150 
1151  void
1152  VectorBase::do_set_add_operation (const size_type n_elements,
1153  const size_type *indices,
1154  const PetscScalar *values,
1155  const bool add_values)
1156  {
1157  ::VectorOperation::values action = (add_values ?
1160  Assert ((last_action == action)
1161  ||
1162  (last_action == ::VectorOperation::unknown),
1163  internal::VectorReference::ExcWrongMode (action,
1164  last_action));
1165  Assert (!has_ghost_elements(), ExcGhostsPresent());
1166  // VecSetValues complains if we
1167  // come with an empty
1168  // vector. however, it is not a
1169  // collective operation, so we
1170  // can skip the call if necessary
1171  // (unlike the above calls)
1172  if (n_elements != 0)
1173  {
1174 #ifdef PETSC_USE_64BIT_INDICES
1175  std::vector<PetscInt> petsc_ind (n_elements);
1176  for (size_type i=0; i<n_elements; ++i)
1177  petsc_ind[i] = indices[i];
1178  const PetscInt *petsc_indices = &petsc_ind[0];
1179 #else
1180  const int *petsc_indices = (const int *)indices;
1181 #endif
1182 
1183  const InsertMode mode = (add_values ? ADD_VALUES : INSERT_VALUES);
1184  const PetscErrorCode ierr = VecSetValues (vector, n_elements, petsc_indices,
1185  values, mode);
1186  AssertThrow (ierr == 0, ExcPETScError(ierr));
1187  }
1188 
1189  // set the mode here, independent of whether we have actually
1190  // written elements or whether the list was empty
1191  last_action = action;
1192  }
1193 
1194 }
1195 
1196 DEAL_II_NAMESPACE_CLOSE
1197 
1198 #endif // DEAL_II_WITH_PETSC
real_type normalize() const 1
#define AssertNothrow(cond, exc)
Definition: exceptions.h:342
VectorBase & operator-=(const VectorBase &V)
void ratio(const VectorBase &a, const VectorBase &b) 1
static ::ExceptionBase & ExcIO()
void sadd(const PetscScalar s, const VectorBase &V)
VectorBase & operator/=(const PetscScalar factor)
#define AssertThrow(cond, exc)
Definition: exceptions.h:369
VectorBase & operator=(const PetscScalar s)
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
void scale(const VectorBase &scaling_factors)
std::pair< size_type, size_type > local_range() const
real_type lp_norm(const real_type p) const
void clear()
Definition: index_set.h:1393
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:313
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()
std_cxx11::enable_if< std_cxx11::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
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:1211
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
#define AssertIsFinite(number)
Definition: exceptions.h:1197
static ::ExceptionBase & ExcInternalError()
bool operator!=(const VectorBase &v) const
void set(const std::vector< size_type > &indices, const std::vector< PetscScalar > &values)