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  else
639  this->_type = PARALLEL;
640  }
641  else if (std::strcmp(ptype,VECSHARED) == 0)
642 #else // PETSc < 3.21.0
643  if ((std::strcmp(ptype,VECSHARED) == 0) || (std::strcmp(ptype,VECMPI) == 0))
644 #endif
645  {
646  ISLocalToGlobalMapping mapping;
647  LibmeshPetscCall(VecGetLocalToGlobalMapping(_vec, &mapping));
648 
649  Vec localrep;
650  LibmeshPetscCall(VecGhostGetLocalForm(_vec,&localrep));
651  // If is a sparsely stored vector, set up our new mapping
652  // Only checking mapping is not enough to determine if a vec is ghosted
653  // We need to check if vec has a local representation
654  if (mapping && localrep)
655  {
656  const numeric_index_type my_local_size = static_cast<numeric_index_type>(petsc_local_size);
657  const numeric_index_type ghost_begin = static_cast<numeric_index_type>(petsc_local_size);
658  PetscInt n;
659  LibmeshPetscCall(ISLocalToGlobalMappingGetSize(mapping, &n));
660 
661  const numeric_index_type ghost_end = static_cast<numeric_index_type>(n);
662  const PetscInt * indices;
663  LibmeshPetscCall(ISLocalToGlobalMappingGetIndices(mapping,&indices));
664 
665  for (numeric_index_type i=ghost_begin; i<ghost_end; i++)
666  _global_to_local_map[indices[i]] = i-my_local_size;
667  this->_type = GHOSTED;
668  LibmeshPetscCall(ISLocalToGlobalMappingRestoreIndices(mapping, &indices));
669  }
670  else
671  this->_type = PARALLEL;
672 
673  LibmeshPetscCall(VecGhostRestoreLocalForm(_vec,&localrep));
674  }
675  else
676  this->_type = SERIAL;
677 }
678 
679 
680 
681 
682 template <typename T>
683 inline
685 {
686  this->clear ();
687 }
688 
689 
690 
691 template <typename T>
692 inline
694  const numeric_index_type n_local,
695  const bool fast,
696  const ParallelType ptype)
697 {
698  parallel_object_only();
699 
700  PetscInt petsc_n=static_cast<PetscInt>(n);
701 
702  // Clear initialized vectors
703  if (this->initialized())
704  this->clear();
705 
706  if (ptype == AUTOMATIC)
707  {
708  if (n == n_local)
709  this->_type = SERIAL;
710  else
711  this->_type = PARALLEL;
712  }
713  else
714  this->_type = ptype;
715 
716  libmesh_assert ((this->_type==SERIAL && n==n_local) ||
717  this->_type==PARALLEL);
718 
719  // create a sequential vector if on only 1 processor
720  if (this->_type == SERIAL)
721  {
722  LibmeshPetscCallA(PETSC_COMM_SELF, VecCreate(PETSC_COMM_SELF, &_vec));
723  LibmeshPetscCallA(PETSC_COMM_SELF, VecSetSizes(_vec, petsc_n, petsc_n));
724  LibmeshPetscCallA(PETSC_COMM_SELF, VecSetFromOptions (_vec));
725  }
726  // otherwise create an MPI-enabled vector
727  else if (this->_type == PARALLEL)
728  {
729 #ifdef LIBMESH_HAVE_MPI
730  PetscInt petsc_n_local=cast_int<PetscInt>(n_local);
731  libmesh_assert_less_equal (n_local, n);
732  // Use more generic function instead of VecCreateSeq/MPI
733  LibmeshPetscCall(VecCreate(this->comm().get(), &_vec));
734  LibmeshPetscCall(VecSetSizes(_vec, petsc_n_local, petsc_n));
735 #else
736  libmesh_assert_equal_to (n_local, n);
737  LibmeshPetscCallA(PETSC_COMM_SELF, VecCreate(PETSC_COMM_SELF, &_vec));
738  LibmeshPetscCallA(PETSC_COMM_SELF, VecSetSizes(_vec, petsc_n, petsc_n));
739 #endif
740  LibmeshPetscCall(VecSetFromOptions(_vec));
741  }
742  else
743  libmesh_error_msg("Unsupported type " << Utility::enum_to_string(this->_type));
744 
745  this->_is_initialized = true;
746  this->_is_closed = true;
747 
748 
749  if (fast == false)
750  this->zero ();
751 }
752 
753 
754 
755 template <typename T>
756 inline
758  const bool fast,
759  const ParallelType ptype)
760 {
761  this->init(n,n,fast,ptype);
762 }
763 
764 
765 
766 template <typename T>
767 inline
769  const numeric_index_type n_local,
770  const std::vector<numeric_index_type> & ghost,
771  const bool fast,
772  const ParallelType libmesh_dbg_var(ptype))
773 {
774  parallel_object_only();
775 
776  PetscInt petsc_n=static_cast<PetscInt>(n);
777  PetscInt petsc_n_local=static_cast<PetscInt>(n_local);
778  PetscInt petsc_n_ghost=static_cast<PetscInt>(ghost.size());
779 
780  // If the mesh is not disjoint, every processor will either have
781  // all the dofs, none of the dofs, or some non-zero dofs at the
782  // boundary between processors.
783  //
784  // However we can't assert this, because someone might want to
785  // construct a GHOSTED vector which doesn't include neighbor element
786  // dofs. Boyce tried to do so in user code, and we're going to want
787  // to do so in System::project_vector().
788  //
789  // libmesh_assert(n_local == 0 || n_local == n || !ghost.empty());
790 
791  libmesh_assert(sizeof(PetscInt) == sizeof(numeric_index_type));
792 
793  PetscInt * petsc_ghost = ghost.empty() ? LIBMESH_PETSC_NULLPTR :
794  const_cast<PetscInt *>(reinterpret_cast<const PetscInt *>(ghost.data()));
795 
796  // Clear initialized vectors
797  if (this->initialized())
798  this->clear();
799 
800  libmesh_assert(ptype == AUTOMATIC || ptype == GHOSTED);
801  this->_type = GHOSTED;
802 
803  /* Make the global-to-local ghost cell map. */
804  for (auto i : index_range(ghost))
805  {
806  _global_to_local_map[ghost[i]] = i;
807  }
808 
809  /* Create vector. */
810  LibmeshPetscCall(VecCreateGhost (this->comm().get(), petsc_n_local, petsc_n,
811  petsc_n_ghost, petsc_ghost,
812  &_vec));
813 
814  // Add a prefix so that we can at least distinguish a ghosted vector from a
815  // nonghosted vector when using a petsc option.
816  // PETSc does not fully support VecGhost on GPU yet. This change allows us to
817  // trigger a nonghosted vector to use GPU without bothering the ghosted vectors.
818  LibmeshPetscCall(PetscObjectAppendOptionsPrefix((PetscObject)_vec,"ghost_"));
819 
820  LibmeshPetscCall(VecSetFromOptions (_vec));
821 
822  this->_is_initialized = true;
823  this->_is_closed = true;
824 
825  if (fast == false)
826  this->zero ();
827 }
828 
829 
830 
831 template <typename T>
832 inline
834  const bool fast)
835 {
836  parallel_object_only();
837 
838  // Clear initialized vectors
839  if (this->initialized())
840  this->clear();
841 
842  const PetscVector<T> & v = cast_ref<const PetscVector<T> &>(other);
843 
844  // Other vector should restore array.
845  if (v.initialized())
846  {
847  v._restore_array();
848  }
849 
850  this->_global_to_local_map = v._global_to_local_map;
851 
852  // Even if we're initializing sizes based on an uninitialized or
853  // unclosed vector, *this* vector is being initialized now and is
854  // initially closed.
855  this->_is_closed = true; // v._is_closed;
856  this->_is_initialized = true; // v._is_initialized;
857 
858  this->_type = v._type;
859 
860  // We want to have a valid Vec, even if it's initially of size zero
861  LibmeshPetscCall(VecDuplicate (v._vec, &this->_vec));
862 
863  if (fast == false)
864  this->zero ();
865 }
866 
867 
868 
869 template <typename T>
870 inline
872 {
873  parallel_object_only();
874 
875  this->_restore_array();
876 
877  VecAssemblyBeginEnd(this->comm(), _vec);
878 
879  if (this->type() == GHOSTED)
880  VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
881 
882  this->_is_closed = true;
883 }
884 
885 
886 
887 template <typename T>
888 inline
889 void PetscVector<T>::clear () noexcept
890 {
891  exceptionless_parallel_object_only();
892 
893  if (this->initialized())
894  this->_restore_array();
895 
896  if ((this->initialized()) && (this->_destroy_vec_on_exit))
897  {
898  // If we encounter an error here, print a warning but otherwise
899  // keep going since we may be recovering from an exception.
900  PetscErrorCode ierr = VecDestroy(&_vec);
901  if (ierr)
902  libmesh_warning("Warning: VecDestroy returned a non-zero error code which we ignored.");
903  }
904 
905  this->_is_closed = this->_is_initialized = false;
906 
907  _global_to_local_map.clear();
908 }
909 
910 
911 
912 template <typename T>
913 inline
915 {
916  parallel_object_only();
917 
918  libmesh_assert(this->closed());
919 
920  this->_restore_array();
921 
922  PetscScalar z=0.;
923 
924  if (this->type() != GHOSTED)
925  LibmeshPetscCall(VecSet (_vec, z));
926  else
927  {
928  /* Vectors that include ghost values require a special
929  handling. */
930  Vec loc_vec;
931  LibmeshPetscCall(VecGhostGetLocalForm (_vec,&loc_vec));
932  LibmeshPetscCall(VecSet (loc_vec, z));
933  LibmeshPetscCall(VecGhostRestoreLocalForm (_vec, &loc_vec));
934  }
935 }
936 
937 
938 
939 template <typename T>
940 inline
941 std::unique_ptr<NumericVector<T>> PetscVector<T>::zero_clone () const
942 {
943  std::unique_ptr<NumericVector<T>> cloned_vector =
944  std::make_unique<PetscVector<T>>(this->comm(), this->type());
945  cloned_vector->init(*this);
946  return cloned_vector;
947 }
948 
949 
950 
951 template <typename T>
952 inline
953 std::unique_ptr<NumericVector<T>> PetscVector<T>::clone () const
954 {
955  std::unique_ptr<NumericVector<T>> cloned_vector =
956  std::make_unique<PetscVector<T>>(this->comm(), this->type());
957  cloned_vector->init(*this, true);
958  *cloned_vector = *this;
959  return cloned_vector;
960 }
961 
962 
963 
964 template <typename T>
965 inline
967 {
968  libmesh_assert (this->initialized());
969 
970  PetscInt petsc_size=0;
971 
972  if (!this->initialized())
973  return 0;
974 
975  LibmeshPetscCall(VecGetSize(_vec, &petsc_size));
976 
977  return static_cast<numeric_index_type>(petsc_size);
978 }
979 
980 
981 
982 template <typename T>
983 inline
985 {
986  libmesh_assert (this->initialized());
987 
988  PetscInt petsc_size=0;
989 
990  LibmeshPetscCall(VecGetLocalSize(_vec, &petsc_size));
991 
992  return static_cast<numeric_index_type>(petsc_size);
993 }
994 
995 
996 
997 template <typename T>
998 inline
1000 {
1001  libmesh_assert (this->initialized());
1002 
1003  numeric_index_type first = 0;
1004 
1005  if (_array_is_present) // Can we use cached values?
1006  first = _first;
1007  else
1008  {
1009  PetscInt petsc_first=0, petsc_last=0;
1010  LibmeshPetscCall(VecGetOwnershipRange (_vec, &petsc_first, &petsc_last));
1011  first = static_cast<numeric_index_type>(petsc_first);
1012  }
1013 
1014  return first;
1015 }
1016 
1017 
1018 
1019 template <typename T>
1020 inline
1022 {
1023  libmesh_assert (this->initialized());
1024 
1025  numeric_index_type last = 0;
1026 
1027  if (_array_is_present) // Can we use cached values?
1028  last = _last;
1029  else
1030  {
1031  PetscInt petsc_first=0, petsc_last=0;
1032  LibmeshPetscCall(VecGetOwnershipRange (_vec, &petsc_first, &petsc_last));
1033  last = static_cast<numeric_index_type>(petsc_last);
1034  }
1035 
1036  return last;
1037 }
1038 
1039 
1040 
1041 template <typename T>
1042 inline
1044 {
1045  libmesh_assert (this->initialized());
1046 
1047  numeric_index_type first=0;
1048  numeric_index_type last=0;
1049 
1050  if (_array_is_present) // Can we use cached values?
1051  {
1052  first = _first;
1053  last = _last;
1054  }
1055  else
1056  {
1057  PetscInt petsc_first=0, petsc_last=0;
1058  LibmeshPetscCall(VecGetOwnershipRange (_vec, &petsc_first, &petsc_last));
1059  first = static_cast<numeric_index_type>(petsc_first);
1060  last = static_cast<numeric_index_type>(petsc_last);
1061  }
1062 
1063 
1064  if ((i>=first) && (i<last))
1065  {
1066  return i-first;
1067  }
1068 
1069  GlobalToLocalMap::const_iterator it = _global_to_local_map.find(i);
1070 #ifndef NDEBUG
1071  const GlobalToLocalMap::const_iterator end = _global_to_local_map.end();
1072  if (it == end)
1073  {
1074  std::ostringstream error_message;
1075  error_message << "No index " << i << " in ghosted vector.\n"
1076  << "Vector contains [" << first << ',' << last << ")\n";
1077  GlobalToLocalMap::const_iterator b = _global_to_local_map.begin();
1078  if (b == end)
1079  {
1080  error_message << "And empty ghost array.\n";
1081  }
1082  else
1083  {
1084  error_message << "And ghost array {" << b->first;
1085  for (++b; b != end; ++b)
1086  error_message << ',' << b->first;
1087  error_message << "}\n";
1088  }
1089 
1090  libmesh_error_msg(error_message.str());
1091  }
1092  libmesh_assert (it != _global_to_local_map.end());
1093 #endif
1094  return it->second+last-first;
1095 }
1096 
1097 
1098 
1099 template <typename T>
1100 inline
1102 {
1103  this->_get_array(true);
1104 
1105  const numeric_index_type local_index = this->map_global_to_local_index(i);
1106 
1107 #ifndef NDEBUG
1108  if (this->type() == GHOSTED)
1109  {
1110  libmesh_assert_less (local_index, _local_size);
1111  }
1112 #endif
1113 
1114  return static_cast<T>(_read_only_values[local_index]);
1115 }
1116 
1117 
1118 
1119 template <typename T>
1120 inline
1121 void PetscVector<T>::get(const std::vector<numeric_index_type> & index,
1122  T * values) const
1123 {
1124  this->_get_array(true);
1125 
1126  const std::size_t num = index.size();
1127 
1128  for (std::size_t i=0; i<num; i++)
1129  {
1130  const numeric_index_type local_index = this->map_global_to_local_index(index[i]);
1131 #ifndef NDEBUG
1132  if (this->type() == GHOSTED)
1133  {
1134  libmesh_assert_less (local_index, _local_size);
1135  }
1136 #endif
1137  values[i] = static_cast<T>(_read_only_values[local_index]);
1138  }
1139 }
1140 
1141 
1142 template <typename T>
1143 inline
1145 {
1146  _get_array(false);
1147  _values_manually_retrieved = true;
1148 
1149  return _values;
1150 }
1151 
1152 
1153 template <typename T>
1154 inline
1155 const PetscScalar * PetscVector<T>::get_array_read() const
1156 {
1157  _get_array(true);
1158  _values_manually_retrieved = true;
1159 
1160  return _read_only_values;
1161 }
1162 
1163 template <typename T>
1164 inline
1166 {
1167  // \note \p _values_manually_retrieved needs to be set to \p false
1168  // \e before calling \p _restore_array()!
1169  _values_manually_retrieved = false;
1170  _restore_array();
1171 }
1172 
1173 template <typename T>
1174 inline
1176 {
1177  parallel_object_only();
1178 
1179  this->_restore_array();
1180 
1181  PetscInt index=0;
1182  PetscReal returnval=0.;
1183 
1184  LibmeshPetscCall(VecMin (_vec, &index, &returnval));
1185 
1186  // this return value is correct: VecMin returns a PetscReal
1187  return static_cast<Real>(returnval);
1188 }
1189 
1190 
1191 
1192 template <typename T>
1193 inline
1195 {
1196  parallel_object_only();
1197 
1198  this->_restore_array();
1199 
1200  PetscInt index=0;
1201  PetscReal returnval=0.;
1202 
1203  LibmeshPetscCall(VecMax (_vec, &index, &returnval));
1204 
1205  // this return value is correct: VecMax returns a PetscReal
1206  return static_cast<Real>(returnval);
1207 }
1208 
1209 
1210 
1211 template <typename T>
1212 inline
1214 {
1215  parallel_object_only();
1216 
1217  NumericVector<T>::swap(other);
1218 
1219  PetscVector<T> & v = cast_ref<PetscVector<T> &>(other);
1220 
1221  std::swap(_vec, v._vec);
1222  std::swap(_destroy_vec_on_exit, v._destroy_vec_on_exit);
1223  std::swap(_global_to_local_map, v._global_to_local_map);
1224 
1225 #ifdef LIBMESH_HAVE_CXX11_THREAD
1226  // Only truly atomic for v... but swap() doesn't really need to be thread safe!
1227  _array_is_present = v._array_is_present.exchange(_array_is_present);
1228 #else
1229  std::swap(_array_is_present, v._array_is_present);
1230 #endif
1231 
1232  std::swap(_local_form, v._local_form);
1233  std::swap(_values, v._values);
1234 }
1235 
1236 
1237 
1238 template <typename T>
1239 inline
1241 {
1242  // The PetscInt type is used for indexing, it may be either a signed
1243  // 4-byte or 8-byte integer depending on how PETSc is configured.
1244  return std::numeric_limits<PetscInt>::max();
1245 }
1246 
1247 
1248 
1249 #ifdef LIBMESH_HAVE_CXX11
1250 static_assert(sizeof(PetscInt) == sizeof(numeric_index_type),
1251  "PETSc and libMesh integer sizes must match!");
1252 #endif
1253 
1254 
1255 inline
1257 {
1258  return reinterpret_cast<PetscInt *>(const_cast<numeric_index_type *>(p));
1259 }
1260 
1261 } // namespace libMesh
1262 
1263 
1264 #endif // #ifdef LIBMESH_HAVE_PETSC
1265 #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:341
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:914
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:966
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:999
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:984
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:315
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:953
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:889
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:693
virtual Real max() const override
static const bool value
Definition: xdr_io.C:55
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:334
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:941
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:871
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