libMesh
petsc_vector.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 
21 #ifndef LIBMESH_PETSC_VECTOR_H
22 #define LIBMESH_PETSC_VECTOR_H
23 
24 
25 #include "libmesh/libmesh_config.h"
26 
27 #ifdef LIBMESH_HAVE_PETSC
28 
29 // Local includes
30 #include "libmesh/numeric_vector.h"
31 #include "libmesh/petsc_macro.h"
32 #include "libmesh/int_range.h"
33 #include "libmesh/libmesh_common.h"
34 #include "libmesh/petsc_solver_exception.h"
35 #include "libmesh/parallel_only.h"
36 #include "libmesh/enum_to_string.h"
37 
38 // PETSc include files.
39 #ifdef I
40 # define LIBMESH_SAW_I
41 #endif
42 #include <petscvec.h>
43 #ifndef LIBMESH_SAW_I
44 # undef I // Avoid complex.h contamination
45 #endif
46 
47 // C++ includes
48 #include <cstddef>
49 #include <cstring>
50 #include <vector>
51 #include <unordered_map>
52 #include <limits>
53 
54 #include <atomic>
55 #include <mutex>
56 #include <condition_variable>
57 
58 namespace libMesh
59 {
60 
61 // forward declarations
62 template <typename T> class SparseMatrix;
63 
72 template <typename T>
73 class PetscVector final : public NumericVector<T>
74 {
75 public:
76 
80  explicit
81  PetscVector (const Parallel::Communicator & comm_in,
82  const ParallelType type = AUTOMATIC);
83 
87  explicit
88  PetscVector (const Parallel::Communicator & comm_in,
89  const numeric_index_type n,
90  const ParallelType type = AUTOMATIC);
91 
96  PetscVector (const Parallel::Communicator & comm_in,
97  const numeric_index_type n,
98  const numeric_index_type n_local,
99  const ParallelType type = AUTOMATIC);
100 
106  PetscVector (const Parallel::Communicator & comm_in,
107  const numeric_index_type N,
108  const numeric_index_type n_local,
109  const std::vector<numeric_index_type> & ghost,
110  const ParallelType type = AUTOMATIC);
111 
119  PetscVector(Vec v,
120  const Parallel::Communicator & comm_in);
121 
128 
134  PetscVector (PetscVector &&) = delete;
135  PetscVector (const PetscVector &) = delete;
136  PetscVector & operator= (PetscVector &&) = delete;
137  virtual ~PetscVector ();
138 
139  virtual void close () override;
140 
144  virtual void clear () noexcept override;
145 
146  virtual void zero () override;
147 
148  virtual std::unique_ptr<NumericVector<T>> zero_clone () const override;
149 
150  virtual std::unique_ptr<NumericVector<T>> clone () const override;
151 
152  virtual void init (const numeric_index_type N,
153  const numeric_index_type n_local,
154  const bool fast=false,
155  const ParallelType type=AUTOMATIC) override;
156 
157  virtual void init (const numeric_index_type N,
158  const bool fast=false,
159  const ParallelType type=AUTOMATIC) override;
160 
161  virtual void init (const numeric_index_type N,
162  const numeric_index_type n_local,
163  const std::vector<numeric_index_type> & ghost,
164  const bool fast = false,
165  const ParallelType = AUTOMATIC) override;
166 
167  virtual void init (const NumericVector<T> & other,
168  const bool fast = false) override;
169 
170  virtual NumericVector<T> & operator= (const T s) override;
171 
172  virtual NumericVector<T> & operator= (const NumericVector<T> & v) override;
173 
174  virtual NumericVector<T> & operator= (const std::vector<T> & v) override;
175 
176  virtual Real min () const override;
177 
178  virtual Real max () const override;
179 
180  virtual T sum () const override;
181 
182  virtual Real l1_norm () const override;
183 
184  virtual Real l2_norm () const override;
185 
186  virtual Real linfty_norm () const override;
187 
188  virtual numeric_index_type size () const override;
189 
190  virtual numeric_index_type local_size() const override;
191 
192  virtual numeric_index_type first_local_index() const override;
193 
194  virtual numeric_index_type last_local_index() const override;
195 
204 
205  virtual T operator() (const numeric_index_type i) const override;
206 
207  virtual void get(const std::vector<numeric_index_type> & index,
208  T * values) const override;
209 
216  PetscScalar * get_array();
217 
224  const PetscScalar * get_array_read() const;
225 
232  void restore_array();
233 
234  virtual NumericVector<T> & operator += (const NumericVector<T> & v) override;
235 
236  virtual NumericVector<T> & operator -= (const NumericVector<T> & v) override;
237 
238  virtual void reciprocal() override;
239 
240  virtual void conjugate() override;
241 
242  virtual void set (const numeric_index_type i,
243  const T value) override;
244 
245  virtual void add (const numeric_index_type i,
246  const T value) override;
247 
248  virtual void add (const T s) override;
249 
250  virtual void add (const NumericVector<T> & v) override;
251 
252  virtual void add (const T a, const NumericVector<T> & v) override;
253 
259 
260  virtual void add_vector (const T * v,
261  const std::vector<numeric_index_type> & dof_indices) override;
262 
263  virtual void add_vector (const NumericVector<T> & v,
264  const SparseMatrix<T> & A) override;
265 
266  virtual void add_vector_transpose (const NumericVector<T> & v,
267  const SparseMatrix<T> & A) override;
268 
276  const SparseMatrix<T> & A);
277 
283 
284  virtual void insert (const T * v,
285  const std::vector<numeric_index_type> & dof_indices) override;
286 
287  virtual void scale (const T factor) override;
288 
289  virtual NumericVector<T> & operator *= (const NumericVector<T> & v) override;
290 
291  virtual NumericVector<T> & operator /= (const NumericVector<T> & v) override;
292 
293  virtual void abs() override;
294 
295  virtual T dot(const NumericVector<T> & v) const override;
296 
302  T indefinite_dot(const NumericVector<T> & v) const;
303 
304  virtual void localize (std::vector<T> & v_local) const override;
305 
306  virtual void localize (NumericVector<T> & v_local) const override;
307 
308  virtual void localize (NumericVector<T> & v_local,
309  const std::vector<numeric_index_type> & send_list) const override;
310 
311  virtual void localize (std::vector<T> & v_local,
312  const std::vector<numeric_index_type> & indices) const override;
313 
314  virtual void localize (const numeric_index_type first_local_idx,
315  const numeric_index_type last_local_idx,
316  const std::vector<numeric_index_type> & send_list) override;
317 
318  virtual void localize_to_one (std::vector<T> & v_local,
319  const processor_id_type proc_id=0) const override;
320 
321  virtual void pointwise_mult (const NumericVector<T> & vec1,
322  const NumericVector<T> & vec2) override;
323 
324  virtual void pointwise_divide (const NumericVector<T> & vec1,
325  const NumericVector<T> & vec2) override;
326 
327  virtual void print_matlab(const std::string & name = "") const override;
328 
329  virtual void create_subvector(NumericVector<T> & subvector,
330  const std::vector<numeric_index_type> & rows,
331  bool supplying_global_rows = true) const override;
332 
333  virtual void swap (NumericVector<T> & v) override;
334 
335  virtual std::size_t max_allowed_id() const override;
336 
345  Vec vec () { libmesh_assert (_vec); return _vec; }
346 
347  Vec vec () const { libmesh_assert (_vec); return _vec; }
348 
349  virtual std::unique_ptr<NumericVector<T>>
350  get_subvector(const std::vector<numeric_index_type> & rows) override;
351 
352  virtual void restore_subvector(std::unique_ptr<NumericVector<T>> subvector,
353  const std::vector<numeric_index_type> & rows) override;
354 
355 private:
356 
360  Vec _vec;
361 
367 #ifdef LIBMESH_HAVE_CXX11_THREAD
368  // We can't use std::atomic_flag here because we need load and store operations.
369  mutable std::atomic<bool> _array_is_present;
370 #else
371  mutable bool _array_is_present;
372 #endif
373 
380 
387 
392 
398  mutable Vec _local_form;
399 
408  mutable const PetscScalar * _read_only_values;
409 
418  mutable PetscScalar * _values;
419 
425  mutable std::mutex _petsc_get_restore_array_mutex;
426 
433  void _get_array(bool read_only) const;
434 
439  void _restore_array() const;
440 
444  typedef std::unordered_map<numeric_index_type,numeric_index_type> GlobalToLocalMap;
445 
451 
457 
463 
467  mutable bool _values_read_only;
468 
475  template <NormType N> Real norm () const;
476 };
477 
478 
479 /*----------------------- Inline functions ----------------------------------*/
480 
481 
482 
483 template <typename T>
484 inline
486  NumericVector<T>(comm_in, ptype),
487  _vec(nullptr),
488  _array_is_present(false),
489  _first(0),
490  _last(0),
491  _local_size(0),
492  _local_form(nullptr),
493  _read_only_values(nullptr),
494  _values(nullptr),
495  _global_to_local_map(),
496  _destroy_vec_on_exit(true),
497  _values_manually_retrieved(false),
498  _values_read_only(false)
499 {
500  this->_type = ptype;
501 }
502 
503 
504 
505 template <typename T>
506 inline
508  const numeric_index_type n,
509  const ParallelType ptype) :
510  NumericVector<T>(comm_in, ptype),
511  _vec(nullptr),
512  _array_is_present(false),
513  _first(0),
514  _last(0),
515  _local_size(0),
516  _local_form(nullptr),
517  _read_only_values(nullptr),
518  _values(nullptr),
519  _global_to_local_map(),
520  _destroy_vec_on_exit(true),
521  _values_manually_retrieved(false),
522  _values_read_only(false)
523 {
524  this->init(n, n, false, ptype);
525 }
526 
527 
528 
529 template <typename T>
530 inline
532  const numeric_index_type n,
533  const numeric_index_type n_local,
534  const ParallelType ptype) :
535  NumericVector<T>(comm_in, ptype),
536  _vec(nullptr),
537  _array_is_present(false),
538  _first(0),
539  _last(0),
540  _local_size(0),
541  _local_form(nullptr),
542  _read_only_values(nullptr),
543  _values(nullptr),
544  _global_to_local_map(),
545  _destroy_vec_on_exit(true),
546  _values_manually_retrieved(false),
547  _values_read_only(false)
548 {
549  this->init(n, n_local, false, ptype);
550 }
551 
552 
553 
554 template <typename T>
555 inline
557  const numeric_index_type n,
558  const numeric_index_type n_local,
559  const std::vector<numeric_index_type> & ghost,
560  const ParallelType ptype) :
561  NumericVector<T>(comm_in, ptype),
562  _vec(nullptr),
563  _array_is_present(false),
564  _first(0),
565  _last(0),
566  _local_size(0),
567  _local_form(nullptr),
568  _read_only_values(nullptr),
569  _values(nullptr),
570  _global_to_local_map(),
571  _destroy_vec_on_exit(true),
572  _values_manually_retrieved(false),
573  _values_read_only(false)
574 {
575  this->init(n, n_local, ghost, false, ptype);
576 }
577 
578 
579 
580 
581 
582 template <typename T>
583 inline
585  const Parallel::Communicator & comm_in) :
586  NumericVector<T>(comm_in, AUTOMATIC),
587  _vec(v),
588  _array_is_present(false),
589  _first(0),
590  _last(0),
591  _local_size(0),
592  _local_form(nullptr),
593  _read_only_values(nullptr),
594  _values(nullptr),
595  _global_to_local_map(),
596  _destroy_vec_on_exit(false),
597  _values_manually_retrieved(false),
598  _values_read_only(false)
599 {
600  this->_is_closed = true;
601  this->_is_initialized = true;
602 
603  /* We need to ask PETSc about the (local to global) ghost value
604  mapping and create the inverse mapping out of it. */
605  PetscInt petsc_local_size=0;
606  LibmeshPetscCall(VecGetLocalSize(_vec, &petsc_local_size));
607 
608  // Get the vector type from PETSc.
609  VecType ptype;
610  LibmeshPetscCall(VecGetType(_vec, &ptype));
611 
612 #if PETSC_RELEASE_GREATER_EQUALS(3, 21, 0)
613  // Fande only implemented VecGhostGetGhostIS for VECMPI
614  if (std::strcmp(ptype, VECMPI) == 0)
615  {
616  IS ghostis;
617  LibmeshPetscCall(VecGhostGetGhostIS(_vec, &ghostis));
618 
619  Vec localrep;
620  LibmeshPetscCall(VecGhostGetLocalForm(_vec, &localrep));
621 
622  // If is a sparsely stored vector, set up our new mapping
623  // Only checking mapping is not enough to determine if a vec is ghosted
624  // We need to check if vec has a local representation
625  if (ghostis && localrep)
626  {
627  PetscInt ghost_size;
628  LibmeshPetscCall(ISGetSize(ghostis, &ghost_size));
629 
630  const PetscInt * indices;
631  LibmeshPetscCall(ISGetIndices(ghostis, &indices));
632 
633  for (const auto i : make_range(ghost_size))
634  _global_to_local_map[indices[i]] = i;
635  this->_type = GHOSTED;
636  LibmeshPetscCall(ISRestoreIndices(ghostis, &indices));
637  }
638  }
639  else if (std::strcmp(ptype,VECSHARED) == 0)
640 #else
641  if ((std::strcmp(ptype,VECSHARED) == 0) || (std::strcmp(ptype,VECMPI) == 0))
642 #endif
643  {
644  ISLocalToGlobalMapping mapping;
645  LibmeshPetscCall(VecGetLocalToGlobalMapping(_vec, &mapping));
646 
647  Vec localrep;
648  LibmeshPetscCall(VecGhostGetLocalForm(_vec,&localrep));
649  // If is a sparsely stored vector, set up our new mapping
650  // Only checking mapping is not enough to determine if a vec is ghosted
651  // We need to check if vec has a local representation
652  if (mapping && localrep)
653  {
654  const numeric_index_type my_local_size = static_cast<numeric_index_type>(petsc_local_size);
655  const numeric_index_type ghost_begin = static_cast<numeric_index_type>(petsc_local_size);
656  PetscInt n;
657  LibmeshPetscCall(ISLocalToGlobalMappingGetSize(mapping, &n));
658 
659  const numeric_index_type ghost_end = static_cast<numeric_index_type>(n);
660  const PetscInt * indices;
661  LibmeshPetscCall(ISLocalToGlobalMappingGetIndices(mapping,&indices));
662 
663  for (numeric_index_type i=ghost_begin; i<ghost_end; i++)
664  _global_to_local_map[indices[i]] = i-my_local_size;
665  this->_type = GHOSTED;
666  LibmeshPetscCall(ISLocalToGlobalMappingRestoreIndices(mapping, &indices));
667  }
668  else
669  this->_type = PARALLEL;
670 
671  LibmeshPetscCall(VecGhostRestoreLocalForm(_vec,&localrep));
672  }
673  else
674  this->_type = SERIAL;
675 }
676 
677 
678 
679 
680 template <typename T>
681 inline
683 {
684  this->clear ();
685 }
686 
687 
688 
689 template <typename T>
690 inline
692  const numeric_index_type n_local,
693  const bool fast,
694  const ParallelType ptype)
695 {
696  parallel_object_only();
697 
698  PetscInt petsc_n=static_cast<PetscInt>(n);
699 
700  // Clear initialized vectors
701  if (this->initialized())
702  this->clear();
703 
704  if (ptype == AUTOMATIC)
705  {
706  if (n == n_local)
707  this->_type = SERIAL;
708  else
709  this->_type = PARALLEL;
710  }
711  else
712  this->_type = ptype;
713 
714  libmesh_assert ((this->_type==SERIAL && n==n_local) ||
715  this->_type==PARALLEL);
716 
717  // create a sequential vector if on only 1 processor
718  if (this->_type == SERIAL)
719  {
720  LibmeshPetscCallA(PETSC_COMM_SELF, VecCreate(PETSC_COMM_SELF, &_vec));
721  LibmeshPetscCallA(PETSC_COMM_SELF, VecSetSizes(_vec, petsc_n, petsc_n));
722  LibmeshPetscCallA(PETSC_COMM_SELF, VecSetFromOptions (_vec));
723  }
724  // otherwise create an MPI-enabled vector
725  else if (this->_type == PARALLEL)
726  {
727 #ifdef LIBMESH_HAVE_MPI
728  PetscInt petsc_n_local=cast_int<PetscInt>(n_local);
729  libmesh_assert_less_equal (n_local, n);
730  // Use more generic function instead of VecCreateSeq/MPI
731  LibmeshPetscCall(VecCreate(this->comm().get(), &_vec));
732  LibmeshPetscCall(VecSetSizes(_vec, petsc_n_local, petsc_n));
733 #else
734  libmesh_assert_equal_to (n_local, n);
735  LibmeshPetscCallA(PETSC_COMM_SELF, VecCreate(PETSC_COMM_SELF, &_vec));
736  LibmeshPetscCallA(PETSC_COMM_SELF, VecSetSizes(_vec, petsc_n, petsc_n));
737 #endif
738  LibmeshPetscCall(VecSetFromOptions(_vec));
739  }
740  else
741  libmesh_error_msg("Unsupported type " << Utility::enum_to_string(this->_type));
742 
743  this->_is_initialized = true;
744  this->_is_closed = true;
745 
746 
747  if (fast == false)
748  this->zero ();
749 }
750 
751 
752 
753 template <typename T>
754 inline
756  const bool fast,
757  const ParallelType ptype)
758 {
759  this->init(n,n,fast,ptype);
760 }
761 
762 
763 
764 template <typename T>
765 inline
767  const numeric_index_type n_local,
768  const std::vector<numeric_index_type> & ghost,
769  const bool fast,
770  const ParallelType libmesh_dbg_var(ptype))
771 {
772  parallel_object_only();
773 
774  PetscInt petsc_n=static_cast<PetscInt>(n);
775  PetscInt petsc_n_local=static_cast<PetscInt>(n_local);
776  PetscInt petsc_n_ghost=static_cast<PetscInt>(ghost.size());
777 
778  // If the mesh is not disjoint, every processor will either have
779  // all the dofs, none of the dofs, or some non-zero dofs at the
780  // boundary between processors.
781  //
782  // However we can't assert this, because someone might want to
783  // construct a GHOSTED vector which doesn't include neighbor element
784  // dofs. Boyce tried to do so in user code, and we're going to want
785  // to do so in System::project_vector().
786  //
787  // libmesh_assert(n_local == 0 || n_local == n || !ghost.empty());
788 
789  libmesh_assert(sizeof(PetscInt) == sizeof(numeric_index_type));
790 
791  PetscInt * petsc_ghost = ghost.empty() ? LIBMESH_PETSC_NULLPTR :
792  const_cast<PetscInt *>(reinterpret_cast<const PetscInt *>(ghost.data()));
793 
794  // Clear initialized vectors
795  if (this->initialized())
796  this->clear();
797 
798  libmesh_assert(ptype == AUTOMATIC || ptype == GHOSTED);
799  this->_type = GHOSTED;
800 
801  /* Make the global-to-local ghost cell map. */
802  for (auto i : index_range(ghost))
803  {
804  _global_to_local_map[ghost[i]] = i;
805  }
806 
807  /* Create vector. */
808  LibmeshPetscCall(VecCreateGhost (this->comm().get(), petsc_n_local, petsc_n,
809  petsc_n_ghost, petsc_ghost,
810  &_vec));
811 
812  // Add a prefix so that we can at least distinguish a ghosted vector from a
813  // nonghosted vector when using a petsc option.
814  // PETSc does not fully support VecGhost on GPU yet. This change allows us to
815  // trigger a nonghosted vector to use GPU without bothering the ghosted vectors.
816  LibmeshPetscCall(PetscObjectAppendOptionsPrefix((PetscObject)_vec,"ghost_"));
817 
818  LibmeshPetscCall(VecSetFromOptions (_vec));
819 
820  this->_is_initialized = true;
821  this->_is_closed = true;
822 
823  if (fast == false)
824  this->zero ();
825 }
826 
827 
828 
829 template <typename T>
830 inline
832  const bool fast)
833 {
834  parallel_object_only();
835 
836  // Clear initialized vectors
837  if (this->initialized())
838  this->clear();
839 
840  const PetscVector<T> & v = cast_ref<const PetscVector<T> &>(other);
841 
842  // Other vector should restore array.
843  if (v.initialized())
844  {
845  v._restore_array();
846  }
847 
848  this->_global_to_local_map = v._global_to_local_map;
849 
850  // Even if we're initializing sizes based on an uninitialized or
851  // unclosed vector, *this* vector is being initialized now and is
852  // initially closed.
853  this->_is_closed = true; // v._is_closed;
854  this->_is_initialized = true; // v._is_initialized;
855 
856  this->_type = v._type;
857 
858  // We want to have a valid Vec, even if it's initially of size zero
859  LibmeshPetscCall(VecDuplicate (v._vec, &this->_vec));
860 
861  if (fast == false)
862  this->zero ();
863 }
864 
865 
866 
867 template <typename T>
868 inline
870 {
871  parallel_object_only();
872 
873  this->_restore_array();
874 
875  VecAssemblyBeginEnd(this->comm(), _vec);
876 
877  if (this->type() == GHOSTED)
878  VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
879 
880  this->_is_closed = true;
881 }
882 
883 
884 
885 template <typename T>
886 inline
887 void PetscVector<T>::clear () noexcept
888 {
889  exceptionless_parallel_object_only();
890 
891  if (this->initialized())
892  this->_restore_array();
893 
894  if ((this->initialized()) && (this->_destroy_vec_on_exit))
895  {
896  // If we encounter an error here, print a warning but otherwise
897  // keep going since we may be recovering from an exception.
898  PetscErrorCode ierr = VecDestroy(&_vec);
899  if (ierr)
900  libmesh_warning("Warning: VecDestroy returned a non-zero error code which we ignored.");
901  }
902 
903  this->_is_closed = this->_is_initialized = false;
904 
905  _global_to_local_map.clear();
906 }
907 
908 
909 
910 template <typename T>
911 inline
913 {
914  parallel_object_only();
915 
916  libmesh_assert(this->closed());
917 
918  this->_restore_array();
919 
920  PetscScalar z=0.;
921 
922  if (this->type() != GHOSTED)
923  LibmeshPetscCall(VecSet (_vec, z));
924  else
925  {
926  /* Vectors that include ghost values require a special
927  handling. */
928  Vec loc_vec;
929  LibmeshPetscCall(VecGhostGetLocalForm (_vec,&loc_vec));
930  LibmeshPetscCall(VecSet (loc_vec, z));
931  LibmeshPetscCall(VecGhostRestoreLocalForm (_vec, &loc_vec));
932  }
933 }
934 
935 
936 
937 template <typename T>
938 inline
939 std::unique_ptr<NumericVector<T>> PetscVector<T>::zero_clone () const
940 {
941  std::unique_ptr<NumericVector<T>> cloned_vector =
942  std::make_unique<PetscVector<T>>(this->comm(), this->type());
943  cloned_vector->init(*this);
944  return cloned_vector;
945 }
946 
947 
948 
949 template <typename T>
950 inline
951 std::unique_ptr<NumericVector<T>> PetscVector<T>::clone () const
952 {
953  std::unique_ptr<NumericVector<T>> cloned_vector =
954  std::make_unique<PetscVector<T>>(this->comm(), this->type());
955  cloned_vector->init(*this, true);
956  *cloned_vector = *this;
957  return cloned_vector;
958 }
959 
960 
961 
962 template <typename T>
963 inline
965 {
966  libmesh_assert (this->initialized());
967 
968  PetscInt petsc_size=0;
969 
970  if (!this->initialized())
971  return 0;
972 
973  LibmeshPetscCall(VecGetSize(_vec, &petsc_size));
974 
975  return static_cast<numeric_index_type>(petsc_size);
976 }
977 
978 
979 
980 template <typename T>
981 inline
983 {
984  libmesh_assert (this->initialized());
985 
986  PetscInt petsc_size=0;
987 
988  LibmeshPetscCall(VecGetLocalSize(_vec, &petsc_size));
989 
990  return static_cast<numeric_index_type>(petsc_size);
991 }
992 
993 
994 
995 template <typename T>
996 inline
998 {
999  libmesh_assert (this->initialized());
1000 
1001  numeric_index_type first = 0;
1002 
1003  if (_array_is_present) // Can we use cached values?
1004  first = _first;
1005  else
1006  {
1007  PetscInt petsc_first=0, petsc_last=0;
1008  LibmeshPetscCall(VecGetOwnershipRange (_vec, &petsc_first, &petsc_last));
1009  first = static_cast<numeric_index_type>(petsc_first);
1010  }
1011 
1012  return first;
1013 }
1014 
1015 
1016 
1017 template <typename T>
1018 inline
1020 {
1021  libmesh_assert (this->initialized());
1022 
1023  numeric_index_type last = 0;
1024 
1025  if (_array_is_present) // Can we use cached values?
1026  last = _last;
1027  else
1028  {
1029  PetscInt petsc_first=0, petsc_last=0;
1030  LibmeshPetscCall(VecGetOwnershipRange (_vec, &petsc_first, &petsc_last));
1031  last = static_cast<numeric_index_type>(petsc_last);
1032  }
1033 
1034  return last;
1035 }
1036 
1037 
1038 
1039 template <typename T>
1040 inline
1042 {
1043  libmesh_assert (this->initialized());
1044 
1045  numeric_index_type first=0;
1046  numeric_index_type last=0;
1047 
1048  if (_array_is_present) // Can we use cached values?
1049  {
1050  first = _first;
1051  last = _last;
1052  }
1053  else
1054  {
1055  PetscInt petsc_first=0, petsc_last=0;
1056  LibmeshPetscCall(VecGetOwnershipRange (_vec, &petsc_first, &petsc_last));
1057  first = static_cast<numeric_index_type>(petsc_first);
1058  last = static_cast<numeric_index_type>(petsc_last);
1059  }
1060 
1061 
1062  if ((i>=first) && (i<last))
1063  {
1064  return i-first;
1065  }
1066 
1067  GlobalToLocalMap::const_iterator it = _global_to_local_map.find(i);
1068 #ifndef NDEBUG
1069  const GlobalToLocalMap::const_iterator end = _global_to_local_map.end();
1070  if (it == end)
1071  {
1072  std::ostringstream error_message;
1073  error_message << "No index " << i << " in ghosted vector.\n"
1074  << "Vector contains [" << first << ',' << last << ")\n";
1075  GlobalToLocalMap::const_iterator b = _global_to_local_map.begin();
1076  if (b == end)
1077  {
1078  error_message << "And empty ghost array.\n";
1079  }
1080  else
1081  {
1082  error_message << "And ghost array {" << b->first;
1083  for (++b; b != end; ++b)
1084  error_message << ',' << b->first;
1085  error_message << "}\n";
1086  }
1087 
1088  libmesh_error_msg(error_message.str());
1089  }
1090  libmesh_assert (it != _global_to_local_map.end());
1091 #endif
1092  return it->second+last-first;
1093 }
1094 
1095 
1096 
1097 template <typename T>
1098 inline
1100 {
1101  this->_get_array(true);
1102 
1103  const numeric_index_type local_index = this->map_global_to_local_index(i);
1104 
1105 #ifndef NDEBUG
1106  if (this->type() == GHOSTED)
1107  {
1108  libmesh_assert_less (local_index, _local_size);
1109  }
1110 #endif
1111 
1112  return static_cast<T>(_read_only_values[local_index]);
1113 }
1114 
1115 
1116 
1117 template <typename T>
1118 inline
1119 void PetscVector<T>::get(const std::vector<numeric_index_type> & index,
1120  T * values) const
1121 {
1122  this->_get_array(true);
1123 
1124  const std::size_t num = index.size();
1125 
1126  for (std::size_t i=0; i<num; i++)
1127  {
1128  const numeric_index_type local_index = this->map_global_to_local_index(index[i]);
1129 #ifndef NDEBUG
1130  if (this->type() == GHOSTED)
1131  {
1132  libmesh_assert_less (local_index, _local_size);
1133  }
1134 #endif
1135  values[i] = static_cast<T>(_read_only_values[local_index]);
1136  }
1137 }
1138 
1139 
1140 template <typename T>
1141 inline
1143 {
1144  _get_array(false);
1145  _values_manually_retrieved = true;
1146 
1147  return _values;
1148 }
1149 
1150 
1151 template <typename T>
1152 inline
1153 const PetscScalar * PetscVector<T>::get_array_read() const
1154 {
1155  _get_array(true);
1156  _values_manually_retrieved = true;
1157 
1158  return _read_only_values;
1159 }
1160 
1161 template <typename T>
1162 inline
1164 {
1165  // \note \p _values_manually_retrieved needs to be set to \p false
1166  // \e before calling \p _restore_array()!
1167  _values_manually_retrieved = false;
1168  _restore_array();
1169 }
1170 
1171 template <typename T>
1172 inline
1174 {
1175  parallel_object_only();
1176 
1177  this->_restore_array();
1178 
1179  PetscInt index=0;
1180  PetscReal returnval=0.;
1181 
1182  LibmeshPetscCall(VecMin (_vec, &index, &returnval));
1183 
1184  // this return value is correct: VecMin returns a PetscReal
1185  return static_cast<Real>(returnval);
1186 }
1187 
1188 
1189 
1190 template <typename T>
1191 inline
1193 {
1194  parallel_object_only();
1195 
1196  this->_restore_array();
1197 
1198  PetscInt index=0;
1199  PetscReal returnval=0.;
1200 
1201  LibmeshPetscCall(VecMax (_vec, &index, &returnval));
1202 
1203  // this return value is correct: VecMax returns a PetscReal
1204  return static_cast<Real>(returnval);
1205 }
1206 
1207 
1208 
1209 template <typename T>
1210 inline
1212 {
1213  parallel_object_only();
1214 
1215  NumericVector<T>::swap(other);
1216 
1217  PetscVector<T> & v = cast_ref<PetscVector<T> &>(other);
1218 
1219  std::swap(_vec, v._vec);
1220  std::swap(_destroy_vec_on_exit, v._destroy_vec_on_exit);
1221  std::swap(_global_to_local_map, v._global_to_local_map);
1222 
1223 #ifdef LIBMESH_HAVE_CXX11_THREAD
1224  // Only truly atomic for v... but swap() doesn't really need to be thread safe!
1225  _array_is_present = v._array_is_present.exchange(_array_is_present);
1226 #else
1227  std::swap(_array_is_present, v._array_is_present);
1228 #endif
1229 
1230  std::swap(_local_form, v._local_form);
1231  std::swap(_values, v._values);
1232 }
1233 
1234 
1235 
1236 template <typename T>
1237 inline
1239 {
1240  // The PetscInt type is used for indexing, it may be either a signed
1241  // 4-byte or 8-byte integer depending on how PETSc is configured.
1242  return std::numeric_limits<PetscInt>::max();
1243 }
1244 
1245 
1246 
1247 #ifdef LIBMESH_HAVE_CXX11
1248 static_assert(sizeof(PetscInt) == sizeof(numeric_index_type),
1249  "PETSc and libMesh integer sizes must match!");
1250 #endif
1251 
1252 
1253 inline
1255 {
1256  return reinterpret_cast<PetscInt *>(const_cast<numeric_index_type *>(p));
1257 }
1258 
1259 } // namespace libMesh
1260 
1261 
1262 #endif // #ifdef LIBMESH_HAVE_PETSC
1263 #endif // LIBMESH_PETSC_VECTOR_H
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
virtual Real l2_norm() const override
Definition: petsc_vector.C:79
PetscVector(const Parallel::Communicator &comm_in, const ParallelType type=AUTOMATIC)
Dummy-Constructor.
Definition: petsc_vector.h:485
bool closed()
Checks that the library has been closed.
Definition: libmesh.C:283
virtual void localize(std::vector< T > &v_local) const override
Creates a copy of the global vector in the local vector v_local.
Definition: petsc_vector.C:787
This class provides a nice interface to PETSc&#39;s Vec object.
Definition: petsc_vector.h:73
void add_vector_conjugate_transpose(const NumericVector< T > &v, const SparseMatrix< T > &A)
.
Definition: petsc_vector.C:260
virtual void add_vector_transpose(const NumericVector< T > &v, const SparseMatrix< T > &A) override
Computes , i.e.
Definition: petsc_vector.C:232
virtual void zero() override
Set all entries to zero.
Definition: petsc_vector.h:912
bool _destroy_vec_on_exit
This boolean value should only be set to false for the constructor which takes a PETSc Vec object...
Definition: petsc_vector.h:456
virtual numeric_index_type size() const override
Definition: petsc_vector.h:964
virtual Real linfty_norm() const override
Definition: petsc_vector.C:84
numeric_index_type _local_size
Size of the local values from _get_array()
Definition: petsc_vector.h:391
virtual numeric_index_type first_local_index() const override
Definition: petsc_vector.h:997
virtual bool initialized() const
void _get_array(bool read_only) const
Queries the array (and the local form if the vector is ghosted) from PETSc.
virtual numeric_index_type local_size() const override
Definition: petsc_vector.h:982
virtual void add(const numeric_index_type i, const T value) override
Adds value to the vector entry specified by i.
Definition: petsc_vector.C:165
Vec _vec
Actual PETSc vector datatype to hold vector entries.
Definition: petsc_vector.h:360
numeric_index_type _last
Last local index.
Definition: petsc_vector.h:386
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
Definition: vector_fe_ex5.C:44
PetscScalar * _values
Pointer to the actual PETSc array of the values of the vector.
Definition: petsc_vector.h:418
bool _values_read_only
Whether or not the data array is for read only access.
Definition: petsc_vector.h:467
The libMesh namespace provides an interface to certain functionality in the library.
const Number zero
.
Definition: libmesh.h:304
virtual void pointwise_divide(const NumericVector< T > &vec1, const NumericVector< T > &vec2) override
Computes (summation not implied) i.e.
void _restore_array() const
Restores the array (and the local form if the vector is ghosted) to PETSc.
std::atomic< bool > _array_is_present
If true, the actual PETSc array of the values of the vector is currently accessible.
Definition: petsc_vector.h:369
uint8_t processor_id_type
Definition: id_types.h:104
PetscInt * numeric_petsc_cast(const numeric_index_type *p)
virtual void swap(NumericVector< T > &v) override
Swaps the contents of this with v.
Generic sparse matrix.
Definition: vector_fe_ex5.C:46
virtual void insert(const T *v, const std::vector< numeric_index_type > &dof_indices) override
Inserts the entries of v in *this at the locations specified by v.
Definition: petsc_vector.C:349
const PetscScalar * _read_only_values
Pointer to the actual PETSc array of the values of the vector.
Definition: petsc_vector.h:408
virtual void pointwise_mult(const NumericVector< T > &vec1, const NumericVector< T > &vec2) override
Computes (summation not implied) i.e.
Definition: petsc_vector.C:989
virtual Real l1_norm() const override
Definition: petsc_vector.C:74
Vec _local_form
PETSc vector datatype to hold the local form of a ghosted vector.
Definition: petsc_vector.h:398
GlobalToLocalMap _global_to_local_map
Map that maps global to local ghost cells (will be empty if not in ghost cell mode).
Definition: petsc_vector.h:450
virtual void abs() override
Sets for each entry in the vector.
Definition: petsc_vector.C:411
dof_id_type numeric_index_type
Definition: id_types.h:99
bool _is_initialized
Flag that tells if init() has been called.
Definition: libmesh.C:257
const PetscScalar * get_array_read() const
Get read only access to the raw PETSc Vector data array.
virtual NumericVector< T > & operator+=(const NumericVector< T > &v) override
Adds v to *this, .
Definition: petsc_vector.C:93
virtual void conjugate() override
Negates the imaginary component of each entry in the vector.
Definition: petsc_vector.C:153
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
libmesh_assert(ctx)
virtual T operator()(const numeric_index_type i) const override
virtual void localize_to_one(std::vector< T > &v_local, const processor_id_type proc_id=0) const override
Creates a local copy of the global vector in v_local only on processor proc_id.
numeric_index_type map_global_to_local_index(const numeric_index_type i) const
virtual std::size_t max_allowed_id() const override
ParallelType _type
Type of vector.
virtual void create_subvector(NumericVector< T > &subvector, const std::vector< numeric_index_type > &rows, bool supplying_global_rows=true) const override
Fills in subvector from this vector using the indices in rows.
virtual numeric_index_type last_local_index() const override
std::string enum_to_string(const T e)
Real norm() const
Definition: petsc_vector.C:60
virtual T sum() const override
Definition: petsc_vector.C:44
std::mutex _petsc_get_restore_array_mutex
Mutex for _get_array and _restore_array.
Definition: petsc_vector.h:425
void restore_array()
Restore the data array.
ParallelType type() const
virtual std::unique_ptr< NumericVector< T > > clone() const override
Definition: petsc_vector.h:951
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
T indefinite_dot(const NumericVector< T > &v) const
Definition: petsc_vector.C:447
std::unordered_map< numeric_index_type, numeric_index_type > GlobalToLocalMap
Type for map that maps global to local ghost cells.
Definition: petsc_vector.h:444
virtual void restore_subvector(std::unique_ptr< NumericVector< T >> subvector, const std::vector< numeric_index_type > &rows) override
Restores a view into this vector using the indices in rows.
virtual void clear() noexcept override
clear() is called from the destructor, so it should not throw.
Definition: petsc_vector.h:887
virtual void init(const numeric_index_type N, const numeric_index_type n_local, const bool fast=false, const ParallelType type=AUTOMATIC) override
Change the dimension of the vector to n.
Definition: petsc_vector.h:691
virtual Real max() const override
static const bool value
Definition: xdr_io.C:54
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
bool initialized()
Checks that library initialization has been done.
Definition: libmesh.C:276
numeric_index_type _first
First local index.
Definition: petsc_vector.h:379
PetscVector< T > & operator=(const PetscVector< T > &v)
Copy assignment operator.
Definition: petsc_vector.C:512
virtual void reciprocal() override
Computes the component-wise reciprocal, .
Definition: petsc_vector.C:141
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices) override
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
Definition: petsc_vector.C:182
virtual NumericVector< T > & operator*=(const NumericVector< T > &v) override
Computes the component-wise multiplication of this vector&#39;s entries by another&#39;s, ...
Definition: petsc_vector.C:385
bool _values_manually_retrieved
Whether or not the data array has been manually retrieved using get_array() or get_array_read() ...
Definition: petsc_vector.h:462
virtual void print_matlab(const std::string &name="") const override
Print the contents of the vector in Matlab&#39;s sparse matrix format.
virtual std::unique_ptr< NumericVector< T > > get_subvector(const std::vector< numeric_index_type > &rows) override
Creates a view into this vector using the indices in rows.
virtual void scale(const T factor) override
Scale each element of the vector by the given factor.
Definition: petsc_vector.C:368
PetscScalar * get_array()
Get read/write access to the raw PETSc Vector data array.
virtual std::unique_ptr< NumericVector< T > > zero_clone() const override
Definition: petsc_vector.h:939
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117
const Elem & get(const ElemType type_in)
virtual Real min() const override
virtual void close() override
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
Definition: petsc_vector.h:869
virtual NumericVector< T > & operator/=(const NumericVector< T > &v) override
Computes the component-wise division of this vector&#39;s entries by another&#39;s, .
Definition: petsc_vector.C:398
ParallelType
Defines an enum for parallel data structure types.
virtual NumericVector< T > & operator-=(const NumericVector< T > &v) override
Subtracts v from *this, .
Definition: petsc_vector.C:109
virtual T dot(const NumericVector< T > &v) const override
Definition: petsc_vector.C:426