libMesh
petsc_vector.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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/libmesh_common.h"
33 #include "libmesh/petsc_solver_exception.h"
34 
35 // PETSc include files.
36 #include <petscvec.h>
37 
38 // C++ includes
39 #include <cstddef>
40 #include <cstring>
41 #include <vector>
42 #include <unordered_map>
43 
44 #ifdef LIBMESH_HAVE_CXX11_THREAD
45 #include <atomic>
46 #include <mutex>
47 #endif
48 
49 namespace libMesh
50 {
51 
52 // forward declarations
53 template <typename T> class SparseMatrix;
54 
63 template <typename T>
64 class PetscVector libmesh_final : public NumericVector<T>
65 {
66 public:
67 
71  explicit
72  PetscVector (const Parallel::Communicator & comm_in,
73  const ParallelType type = AUTOMATIC);
74 
78  explicit
79  PetscVector (const Parallel::Communicator & comm_in,
80  const numeric_index_type n,
81  const ParallelType type = AUTOMATIC);
82 
87  PetscVector (const Parallel::Communicator & comm_in,
88  const numeric_index_type n,
89  const numeric_index_type n_local,
90  const ParallelType type = AUTOMATIC);
91 
97  PetscVector (const Parallel::Communicator & comm_in,
98  const numeric_index_type N,
99  const numeric_index_type n_local,
100  const std::vector<numeric_index_type> & ghost,
101  const ParallelType type = AUTOMATIC);
102 
110  PetscVector(Vec v,
111  const Parallel::Communicator & comm_in
112  LIBMESH_CAN_DEFAULT_TO_COMMWORLD);
113 
118  ~PetscVector ();
119 
120  virtual void close () libmesh_override;
121 
122  virtual void clear () libmesh_override;
123 
124  virtual void zero () libmesh_override;
125 
126  virtual UniquePtr<NumericVector<T>> zero_clone () const libmesh_override;
127 
128  virtual UniquePtr<NumericVector<T>> clone () const libmesh_override;
129 
130  virtual void init (const numeric_index_type N,
131  const numeric_index_type n_local,
132  const bool fast=false,
133  const ParallelType type=AUTOMATIC) libmesh_override;
134 
135  virtual void init (const numeric_index_type N,
136  const bool fast=false,
137  const ParallelType type=AUTOMATIC) libmesh_override;
138 
139  virtual void init (const numeric_index_type N,
140  const numeric_index_type n_local,
141  const std::vector<numeric_index_type> & ghost,
142  const bool fast = false,
143  const ParallelType = AUTOMATIC) libmesh_override;
144 
145  virtual void init (const NumericVector<T> & other,
146  const bool fast = false) libmesh_override;
147 
148  virtual NumericVector<T> & operator= (const T s) libmesh_override;
149 
150  virtual NumericVector<T> & operator= (const NumericVector<T> & v) libmesh_override;
151 
152  virtual NumericVector<T> & operator= (const std::vector<T> & v) libmesh_override;
153 
159  PetscVector<T> & operator= (const PetscVector<T> & v);
160 
161  virtual Real min () const libmesh_override;
162 
163  virtual Real max () const libmesh_override;
164 
165  virtual T sum () const libmesh_override;
166 
167  virtual Real l1_norm () const libmesh_override;
168 
169  virtual Real l2_norm () const libmesh_override;
170 
171  virtual Real linfty_norm () const libmesh_override;
172 
173  virtual numeric_index_type size () const libmesh_override;
174 
175  virtual numeric_index_type local_size() const libmesh_override;
176 
177  virtual numeric_index_type first_local_index() const libmesh_override;
178 
179  virtual numeric_index_type last_local_index() const libmesh_override;
180 
188  numeric_index_type map_global_to_local_index(const numeric_index_type i) const;
189 
190  virtual T operator() (const numeric_index_type i) const libmesh_override;
191 
192  virtual void get(const std::vector<numeric_index_type> & index,
193  T * values) const libmesh_override;
194 
201  PetscScalar * get_array();
202 
209  const PetscScalar * get_array_read() const;
210 
217  void restore_array();
218 
219  virtual NumericVector<T> & operator += (const NumericVector<T> & v) libmesh_override;
220 
221  virtual NumericVector<T> & operator -= (const NumericVector<T> & v) libmesh_override;
222 
223  virtual void reciprocal() libmesh_override;
224 
225  virtual void conjugate() libmesh_override;
226 
227  virtual void set (const numeric_index_type i,
228  const T value) libmesh_override;
229 
230  virtual void add (const numeric_index_type i,
231  const T value) libmesh_override;
232 
233  virtual void add (const T s) libmesh_override;
234 
235  virtual void add (const NumericVector<T> & v) libmesh_override;
236 
237  virtual void add (const T a, const NumericVector<T> & v) libmesh_override;
238 
243  using NumericVector<T>::add_vector;
244 
245  virtual void add_vector (const T * v,
246  const std::vector<numeric_index_type> & dof_indices) libmesh_override;
247 
248  virtual void add_vector (const NumericVector<T> & v,
249  const SparseMatrix<T> & A) libmesh_override;
250 
251  virtual void add_vector_transpose (const NumericVector<T> & v,
252  const SparseMatrix<T> & A) libmesh_override;
253 
260  void add_vector_conjugate_transpose (const NumericVector<T> & v,
261  const SparseMatrix<T> & A);
262 
267  using NumericVector<T>::insert;
268 
269  virtual void insert (const T * v,
270  const std::vector<numeric_index_type> & dof_indices) libmesh_override;
271 
272  virtual void scale (const T factor) libmesh_override;
273 
274  virtual NumericVector<T> & operator /= (NumericVector<T> & v) libmesh_override;
275 
276  virtual void abs() libmesh_override;
277 
278  virtual T dot(const NumericVector<T> & v) const libmesh_override;
279 
285  T indefinite_dot(const NumericVector<T> & v) const;
286 
287  virtual void localize (std::vector<T> & v_local) const libmesh_override;
288 
289  virtual void localize (NumericVector<T> & v_local) const libmesh_override;
290 
291  virtual void localize (NumericVector<T> & v_local,
292  const std::vector<numeric_index_type> & send_list) const libmesh_override;
293 
294  virtual void localize (std::vector<T> & v_local,
295  const std::vector<numeric_index_type> & indices) const libmesh_override;
296 
297  virtual void localize (const numeric_index_type first_local_idx,
298  const numeric_index_type last_local_idx,
299  const std::vector<numeric_index_type> & send_list) libmesh_override;
300 
301  virtual void localize_to_one (std::vector<T> & v_local,
302  const processor_id_type proc_id=0) const libmesh_override;
303 
304  virtual void pointwise_mult (const NumericVector<T> & vec1,
305  const NumericVector<T> & vec2) libmesh_override;
306 
307  virtual void print_matlab(const std::string & name = "") const libmesh_override;
308 
309  virtual void create_subvector(NumericVector<T> & subvector,
310  const std::vector<numeric_index_type> & rows) const libmesh_override;
311 
312  virtual void swap (NumericVector<T> & v) libmesh_override;
313 
322  Vec vec () { libmesh_assert (_vec); return _vec; }
323 
324 
325 private:
326 
330  Vec _vec;
331 
337 #ifdef LIBMESH_HAVE_CXX11_THREAD
338  // We can't use std::atomic_flag here because we need load and store operations.
339  mutable std::atomic<bool> _array_is_present;
340 #else
341  mutable bool _array_is_present;
342 #endif
343 
350 
357 
361  mutable numeric_index_type _local_size;
362 
368  mutable Vec _local_form;
369 
378  mutable const PetscScalar * _read_only_values;
379 
388  mutable PetscScalar * _values;
389 
395 #ifdef LIBMESH_HAVE_CXX11_THREAD
396  mutable std::mutex _petsc_vector_mutex;
397 #else
399 #endif
400 
407  void _get_array(bool read_only) const;
408 
413  void _restore_array() const;
414 
418  typedef std::unordered_map<numeric_index_type,numeric_index_type> GlobalToLocalMap;
419 
424  GlobalToLocalMap _global_to_local_map;
425 
431 
437 
441  mutable bool _values_read_only;
442 };
443 
444 
445 /*----------------------- Inline functions ----------------------------------*/
446 
447 
448 
449 template <typename T>
450 inline
451 PetscVector<T>::PetscVector (const Parallel::Communicator & comm_in, const ParallelType ptype) :
452  NumericVector<T>(comm_in, ptype),
453  _array_is_present(false),
454  _first(0),
455  _last(0),
456  _local_form(libmesh_nullptr),
457  _values(libmesh_nullptr),
458  _global_to_local_map(),
459  _destroy_vec_on_exit(true),
460  _values_manually_retrieved(false),
461  _values_read_only(false)
462 {
463  this->_type = ptype;
464 }
465 
466 
467 
468 template <typename T>
469 inline
470 PetscVector<T>::PetscVector (const Parallel::Communicator & comm_in,
471  const numeric_index_type n,
472  const ParallelType ptype) :
473  NumericVector<T>(comm_in, ptype),
474  _array_is_present(false),
478  _destroy_vec_on_exit(true),
480  _values_read_only(false)
481 {
482  this->init(n, n, false, ptype);
483 }
484 
485 
486 
487 template <typename T>
488 inline
489 PetscVector<T>::PetscVector (const Parallel::Communicator & comm_in,
490  const numeric_index_type n,
491  const numeric_index_type n_local,
492  const ParallelType ptype) :
493  NumericVector<T>(comm_in, ptype),
494  _array_is_present(false),
498  _destroy_vec_on_exit(true),
500  _values_read_only(false)
501 {
502  this->init(n, n_local, false, ptype);
503 }
504 
505 
506 
507 template <typename T>
508 inline
509 PetscVector<T>::PetscVector (const Parallel::Communicator & comm_in,
510  const numeric_index_type n,
511  const numeric_index_type n_local,
512  const std::vector<numeric_index_type> & ghost,
513  const ParallelType ptype) :
514  NumericVector<T>(comm_in, ptype),
515  _array_is_present(false),
519  _destroy_vec_on_exit(true),
521  _values_read_only(false)
522 {
523  this->init(n, n_local, ghost, false, ptype);
524 }
525 
526 
527 
528 
529 
530 template <typename T>
531 inline
532 PetscVector<T>::PetscVector (Vec v,
533  const Parallel::Communicator & comm_in) :
534  NumericVector<T>(comm_in, AUTOMATIC),
535  _array_is_present(false),
539  _destroy_vec_on_exit(false),
541  _values_read_only(false)
542 {
543  this->_vec = v;
544  this->_is_closed = true;
545  this->_is_initialized = true;
546 
547  /* We need to ask PETSc about the (local to global) ghost value
548  mapping and create the inverse mapping out of it. */
549  PetscErrorCode ierr=0;
550  PetscInt petsc_local_size=0;
551  ierr = VecGetLocalSize(_vec, &petsc_local_size);
552  LIBMESH_CHKERR(ierr);
553 
554  // Get the vector type from PETSc.
555  // As of PETSc 3.0.0, the VecType #define lost its const-ness, so we
556  // need to have it in the code
557 #if PETSC_VERSION_LESS_THAN(3,0,0) || !PETSC_VERSION_LESS_THAN(3,4,0)
558  // Pre-3.0 and petsc-dev (as of October 2012) use non-const versions
559  VecType ptype;
560 #else
561  const VecType ptype;
562 #endif
563  ierr = VecGetType(_vec, &ptype);
564  LIBMESH_CHKERR(ierr);
565 
566  if ((std::strcmp(ptype,VECSHARED) == 0) || (std::strcmp(ptype,VECMPI) == 0))
567  {
568 #if PETSC_RELEASE_LESS_THAN(3,1,1)
569  ISLocalToGlobalMapping mapping = _vec->mapping;
570 #else
571  ISLocalToGlobalMapping mapping;
572  ierr = VecGetLocalToGlobalMapping(_vec, &mapping);
573  LIBMESH_CHKERR(ierr);
574 #endif
575 
576  // If is a sparsely stored vector, set up our new mapping
577  if (mapping)
578  {
579  const numeric_index_type my_local_size = static_cast<numeric_index_type>(petsc_local_size);
580  const numeric_index_type ghost_begin = static_cast<numeric_index_type>(petsc_local_size);
581 #if PETSC_RELEASE_LESS_THAN(3,4,0)
582  const numeric_index_type ghost_end = static_cast<numeric_index_type>(mapping->n);
583 #else
584  PetscInt n;
585  ierr = ISLocalToGlobalMappingGetSize(mapping, &n);
586  LIBMESH_CHKERR(ierr);
587  const numeric_index_type ghost_end = static_cast<numeric_index_type>(n);
588 #endif
589 #if PETSC_RELEASE_LESS_THAN(3,1,1)
590  const PetscInt * indices = mapping->indices;
591 #else
592  const PetscInt * indices;
593  ierr = ISLocalToGlobalMappingGetIndices(mapping,&indices);
594  LIBMESH_CHKERR(ierr);
595 #endif
596  for (numeric_index_type i=ghost_begin; i<ghost_end; i++)
597  _global_to_local_map[indices[i]] = i-my_local_size;
598  this->_type = GHOSTED;
599 #if !PETSC_RELEASE_LESS_THAN(3,1,1)
600  ierr = ISLocalToGlobalMappingRestoreIndices(mapping, &indices);
601  LIBMESH_CHKERR(ierr);
602 #endif
603  }
604  else
605  this->_type = PARALLEL;
606  }
607  else
608  this->_type = SERIAL;
609 }
610 
611 
612 
613 
614 template <typename T>
615 inline
616 PetscVector<T>::~PetscVector ()
617 {
618  this->clear ();
619 }
620 
621 
622 
623 template <typename T>
624 inline
626  const numeric_index_type n_local,
627  const bool fast,
628  const ParallelType ptype)
629 {
630  parallel_object_only();
631 
632  PetscErrorCode ierr=0;
633  PetscInt petsc_n=static_cast<PetscInt>(n);
634  PetscInt petsc_n_local=static_cast<PetscInt>(n_local);
635 
636 
637  // Clear initialized vectors
638  if (this->initialized())
639  this->clear();
640 
641  if (ptype == AUTOMATIC)
642  {
643  if (n == n_local)
644  this->_type = SERIAL;
645  else
646  this->_type = PARALLEL;
647  }
648  else
649  this->_type = ptype;
650 
651  libmesh_assert ((this->_type==SERIAL && n==n_local) ||
652  this->_type==PARALLEL);
653 
654  // create a sequential vector if on only 1 processor
655  if (this->_type == SERIAL)
656  {
657  ierr = VecCreateSeq (PETSC_COMM_SELF, petsc_n, &_vec);
658  CHKERRABORT(PETSC_COMM_SELF,ierr);
659 
660  ierr = VecSetFromOptions (_vec);
661  CHKERRABORT(PETSC_COMM_SELF,ierr);
662  }
663  // otherwise create an MPI-enabled vector
664  else if (this->_type == PARALLEL)
665  {
666 #ifdef LIBMESH_HAVE_MPI
667  libmesh_assert_less_equal (n_local, n);
668  ierr = VecCreateMPI (this->comm().get(), petsc_n_local, petsc_n,
669  &_vec);
670  LIBMESH_CHKERR(ierr);
671 #else
672  libmesh_assert_equal_to (n_local, n);
673  ierr = VecCreateSeq (PETSC_COMM_SELF, petsc_n, &_vec);
674  CHKERRABORT(PETSC_COMM_SELF,ierr);
675 #endif
676 
677  ierr = VecSetFromOptions (_vec);
678  LIBMESH_CHKERR(ierr);
679  }
680  else
681  libmesh_error_msg("Unsupported type " << this->_type);
682 
683  this->_is_initialized = true;
684  this->_is_closed = true;
685 
686 
687  if (fast == false)
688  this->zero ();
689 }
690 
691 
692 
693 template <typename T>
694 inline
696  const bool fast,
697  const ParallelType ptype)
698 {
699  this->init(n,n,fast,ptype);
700 }
701 
702 
703 
704 template <typename T>
705 inline
707  const numeric_index_type n_local,
708  const std::vector<numeric_index_type> & ghost,
709  const bool fast,
710  const ParallelType libmesh_dbg_var(ptype))
711 {
712  parallel_object_only();
713 
714  PetscErrorCode ierr=0;
715  PetscInt petsc_n=static_cast<PetscInt>(n);
716  PetscInt petsc_n_local=static_cast<PetscInt>(n_local);
717  PetscInt petsc_n_ghost=static_cast<PetscInt>(ghost.size());
718 
719  // If the mesh is not disjoint, every processor will either have
720  // all the dofs, none of the dofs, or some non-zero dofs at the
721  // boundary between processors.
722  //
723  // However we can't assert this, because someone might want to
724  // construct a GHOSTED vector which doesn't include neighbor element
725  // dofs. Boyce tried to do so in user code, and we're going to want
726  // to do so in System::project_vector().
727  //
728  // libmesh_assert(n_local == 0 || n_local == n || !ghost.empty());
729 
730  libmesh_assert(sizeof(PetscInt) == sizeof(numeric_index_type));
731 
732  PetscInt * petsc_ghost = ghost.empty() ? PETSC_NULL :
733  const_cast<PetscInt *>(reinterpret_cast<const PetscInt *>(&ghost[0]));
734 
735  // Clear initialized vectors
736  if (this->initialized())
737  this->clear();
738 
739  libmesh_assert(ptype == AUTOMATIC || ptype == GHOSTED);
740  this->_type = GHOSTED;
741 
742  /* Make the global-to-local ghost cell map. */
743  for (numeric_index_type i=0; i<ghost.size(); i++)
744  {
745  _global_to_local_map[ghost[i]] = i;
746  }
747 
748  /* Create vector. */
749  ierr = VecCreateGhost (this->comm().get(), petsc_n_local, petsc_n,
750  petsc_n_ghost, petsc_ghost,
751  &_vec);
752  LIBMESH_CHKERR(ierr);
753 
754  ierr = VecSetFromOptions (_vec);
755  LIBMESH_CHKERR(ierr);
756 
757  this->_is_initialized = true;
758  this->_is_closed = true;
759 
760  if (fast == false)
761  this->zero ();
762 }
763 
764 
765 
766 template <typename T>
767 inline
768 void PetscVector<T>::init (const NumericVector<T> & other,
769  const bool fast)
770 {
771  parallel_object_only();
772 
773  // Clear initialized vectors
774  if (this->initialized())
775  this->clear();
776 
777  const PetscVector<T> & v = cast_ref<const PetscVector<T> &>(other);
778 
779  // Other vector should restore array.
780  if (v.initialized())
781  {
782  v._restore_array();
783  }
784 
785  this->_global_to_local_map = v._global_to_local_map;
786 
787  // Even if we're initializing sizes based on an uninitialized or
788  // unclosed vector, *this* vector is being initialized now and is
789  // initially closed.
790  this->_is_closed = true; // v._is_closed;
791  this->_is_initialized = true; // v._is_initialized;
792 
793  this->_type = v._type;
794 
795  // We want to have a valid Vec, even if it's initially of size zero
796  PetscErrorCode ierr = VecDuplicate (v._vec, &this->_vec);
797  LIBMESH_CHKERR(ierr);
798 
799  if (fast == false)
800  this->zero ();
801 }
802 
803 
804 
805 template <typename T>
806 inline
807 void PetscVector<T>::close ()
808 {
809  parallel_object_only();
810 
811  this->_restore_array();
812 
813  PetscErrorCode ierr=0;
814 
815  ierr = VecAssemblyBegin(_vec);
816  LIBMESH_CHKERR(ierr);
817  ierr = VecAssemblyEnd(_vec);
818  LIBMESH_CHKERR(ierr);
819 
820  if (this->type() == GHOSTED)
821  {
822  ierr = VecGhostUpdateBegin(_vec,INSERT_VALUES,SCATTER_FORWARD);
823  LIBMESH_CHKERR(ierr);
824  ierr = VecGhostUpdateEnd(_vec,INSERT_VALUES,SCATTER_FORWARD);
825  LIBMESH_CHKERR(ierr);
826  }
827 
828  this->_is_closed = true;
829 }
830 
831 
832 
833 template <typename T>
834 inline
835 void PetscVector<T>::clear ()
836 {
837  parallel_object_only();
838 
839  if (this->initialized())
840  this->_restore_array();
841 
842  if ((this->initialized()) && (this->_destroy_vec_on_exit))
843  {
844  PetscErrorCode ierr=0;
845 
846  ierr = LibMeshVecDestroy(&_vec);
847  LIBMESH_CHKERR(ierr);
848  }
849 
850  this->_is_closed = this->_is_initialized = false;
851 
852  _global_to_local_map.clear();
853 }
854 
855 
856 
857 template <typename T>
858 inline
859 void PetscVector<T>::zero ()
860 {
861  parallel_object_only();
862 
863  libmesh_assert(this->closed());
864 
865  this->_restore_array();
866 
867  PetscErrorCode ierr=0;
868 
869  PetscScalar z=0.;
870 
871  if (this->type() != GHOSTED)
872  {
873  ierr = VecSet (_vec, z);
874  LIBMESH_CHKERR(ierr);
875  }
876  else
877  {
878  /* Vectors that include ghost values require a special
879  handling. */
880  Vec loc_vec;
881  ierr = VecGhostGetLocalForm (_vec,&loc_vec);
882  LIBMESH_CHKERR(ierr);
883 
884  ierr = VecSet (loc_vec, z);
885  LIBMESH_CHKERR(ierr);
886 
887  ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
888  LIBMESH_CHKERR(ierr);
889  }
890 }
891 
892 
893 
894 template <typename T>
895 inline
896 UniquePtr<NumericVector<T>> PetscVector<T>::zero_clone () const
897 {
898  NumericVector<T> * cloned_vector = new PetscVector<T>(this->comm(), this->type());
899  cloned_vector->init(*this);
900  return UniquePtr<NumericVector<T>>(cloned_vector);
901 }
902 
903 
904 
905 template <typename T>
906 inline
907 UniquePtr<NumericVector<T>> PetscVector<T>::clone () const
908 {
909  NumericVector<T> * cloned_vector = new PetscVector<T>(this->comm(), this->type());
910  cloned_vector->init(*this, true);
911  *cloned_vector = *this;
912  return UniquePtr<NumericVector<T>>(cloned_vector);
913 }
914 
915 
916 
917 template <typename T>
918 inline
919 numeric_index_type PetscVector<T>::size () const
920 {
921  libmesh_assert (this->initialized());
922 
923  PetscErrorCode ierr=0;
924  PetscInt petsc_size=0;
925 
926  if (!this->initialized())
927  return 0;
928 
929  ierr = VecGetSize(_vec, &petsc_size);
930  LIBMESH_CHKERR(ierr);
931 
932  return static_cast<numeric_index_type>(petsc_size);
933 }
934 
935 
936 
937 template <typename T>
938 inline
939 numeric_index_type PetscVector<T>::local_size () const
940 {
941  libmesh_assert (this->initialized());
942 
943  PetscErrorCode ierr=0;
944  PetscInt petsc_size=0;
945 
946  ierr = VecGetLocalSize(_vec, &petsc_size);
947  LIBMESH_CHKERR(ierr);
948 
949  return static_cast<numeric_index_type>(petsc_size);
950 }
951 
952 
953 
954 template <typename T>
955 inline
956 numeric_index_type PetscVector<T>::first_local_index () const
957 {
958  libmesh_assert (this->initialized());
959 
960  numeric_index_type first = 0;
961 
962  if (_array_is_present) // Can we use cached values?
963  first = _first;
964  else
965  {
966  PetscErrorCode ierr=0;
967  PetscInt petsc_first=0, petsc_last=0;
968  ierr = VecGetOwnershipRange (_vec, &petsc_first, &petsc_last);
969  LIBMESH_CHKERR(ierr);
970  first = static_cast<numeric_index_type>(petsc_first);
971  }
972 
973  return first;
974 }
975 
976 
977 
978 template <typename T>
979 inline
980 numeric_index_type PetscVector<T>::last_local_index () const
981 {
982  libmesh_assert (this->initialized());
983 
984  numeric_index_type last = 0;
985 
986  if (_array_is_present) // Can we use cached values?
987  last = _last;
988  else
989  {
990  PetscErrorCode ierr=0;
991  PetscInt petsc_first=0, petsc_last=0;
992  ierr = VecGetOwnershipRange (_vec, &petsc_first, &petsc_last);
993  LIBMESH_CHKERR(ierr);
994  last = static_cast<numeric_index_type>(petsc_last);
995  }
996 
997  return last;
998 }
999 
1000 
1001 
1002 template <typename T>
1003 inline
1004 numeric_index_type PetscVector<T>::map_global_to_local_index (const numeric_index_type i) const
1005 {
1006  libmesh_assert (this->initialized());
1007 
1008  numeric_index_type first=0;
1009  numeric_index_type last=0;
1010 
1011  if (_array_is_present) // Can we use cached values?
1012  {
1013  first = _first;
1014  last = _last;
1015  }
1016  else
1017  {
1018  PetscErrorCode ierr=0;
1019  PetscInt petsc_first=0, petsc_last=0;
1020  ierr = VecGetOwnershipRange (_vec, &petsc_first, &petsc_last);
1021  LIBMESH_CHKERR(ierr);
1022  first = static_cast<numeric_index_type>(petsc_first);
1023  last = static_cast<numeric_index_type>(petsc_last);
1024  }
1025 
1026 
1027  if ((i>=first) && (i<last))
1028  {
1029  return i-first;
1030  }
1031 
1032  GlobalToLocalMap::const_iterator it = _global_to_local_map.find(i);
1033 #ifndef NDEBUG
1034  const GlobalToLocalMap::const_iterator end = _global_to_local_map.end();
1035  if (it == end)
1036  {
1037  std::ostringstream error_message;
1038  error_message << "No index " << i << " in ghosted vector.\n"
1039  << "Vector contains [" << first << ',' << last << ")\n";
1040  GlobalToLocalMap::const_iterator b = _global_to_local_map.begin();
1041  if (b == end)
1042  {
1043  error_message << "And empty ghost array.\n";
1044  }
1045  else
1046  {
1047  error_message << "And ghost array {" << b->first;
1048  for (++b; b != end; ++b)
1049  error_message << ',' << b->first;
1050  error_message << "}\n";
1051  }
1052 
1053  libmesh_error_msg(error_message.str());
1054  }
1055  libmesh_assert (it != _global_to_local_map.end());
1056 #endif
1057  return it->second+last-first;
1058 }
1059 
1060 
1061 
1062 template <typename T>
1063 inline
1064 T PetscVector<T>::operator() (const numeric_index_type i) const
1065 {
1066  this->_get_array(true);
1067 
1068  const numeric_index_type local_index = this->map_global_to_local_index(i);
1069 
1070 #ifndef NDEBUG
1071  if (this->type() == GHOSTED)
1072  {
1073  libmesh_assert_less (local_index, _local_size);
1074  }
1075 #endif
1076 
1077  return static_cast<T>(_read_only_values[local_index]);
1078 }
1079 
1080 
1081 
1082 template <typename T>
1083 inline
1084 void PetscVector<T>::get(const std::vector<numeric_index_type> & index,
1085  T * values) const
1086 {
1087  this->_get_array(true);
1088 
1089  const std::size_t num = index.size();
1090 
1091  for (std::size_t i=0; i<num; i++)
1092  {
1093  const numeric_index_type local_index = this->map_global_to_local_index(index[i]);
1094 #ifndef NDEBUG
1095  if (this->type() == GHOSTED)
1096  {
1097  libmesh_assert_less (local_index, _local_size);
1098  }
1099 #endif
1100  values[i] = static_cast<T>(_read_only_values[local_index]);
1101  }
1102 }
1103 
1104 
1105 template <typename T>
1106 inline
1107 PetscScalar * PetscVector<T>::get_array()
1108 {
1110  _get_array(false);
1111 
1112  return _values;
1113 }
1114 
1115 
1116 template <typename T>
1117 inline
1118 const PetscScalar * PetscVector<T>::get_array_read() const
1119 {
1121  _get_array(true);
1122 
1123  return _read_only_values;
1124 }
1125 
1126 template <typename T>
1127 inline
1128 void PetscVector<T>::restore_array()
1129 {
1130  // \note \p _values_manually_retrieved needs to be set to \p false
1131  // \e before calling \p _restore_array()!
1133  _restore_array();
1134 }
1135 
1136 template <typename T>
1137 inline
1138 Real PetscVector<T>::min () const
1139 {
1140  parallel_object_only();
1141 
1142  this->_restore_array();
1143 
1144  PetscErrorCode ierr=0;
1145  PetscInt index=0;
1146  PetscReal returnval=0.;
1147 
1148  ierr = VecMin (_vec, &index, &returnval);
1149  LIBMESH_CHKERR(ierr);
1150 
1151  // this return value is correct: VecMin returns a PetscReal
1152  return static_cast<Real>(returnval);
1153 }
1154 
1155 
1156 
1157 template <typename T>
1158 inline
1159 Real PetscVector<T>::max() const
1160 {
1161  parallel_object_only();
1162 
1163  this->_restore_array();
1164 
1165  PetscErrorCode ierr=0;
1166  PetscInt index=0;
1167  PetscReal returnval=0.;
1168 
1169  ierr = VecMax (_vec, &index, &returnval);
1170  LIBMESH_CHKERR(ierr);
1171 
1172  // this return value is correct: VecMax returns a PetscReal
1173  return static_cast<Real>(returnval);
1174 }
1175 
1176 
1177 
1178 template <typename T>
1179 inline
1181 {
1182  parallel_object_only();
1183 
1184  NumericVector<T>::swap(other);
1185 
1186  PetscVector<T> & v = cast_ref<PetscVector<T> &>(other);
1187 
1188  std::swap(_vec, v._vec);
1189  std::swap(_destroy_vec_on_exit, v._destroy_vec_on_exit);
1190  std::swap(_global_to_local_map, v._global_to_local_map);
1191 
1192 #ifdef LIBMESH_HAVE_CXX11_THREAD
1193  // Only truly atomic for v... but swap() doesn't really need to be thread safe!
1194  _array_is_present = v._array_is_present.exchange(_array_is_present);
1195 #else
1196  std::swap(_array_is_present, v._array_is_present);
1197 #endif
1198 
1199  std::swap(_local_form, v._local_form);
1200  std::swap(_values, v._values);
1201 }
1202 
1203 
1204 
1205 
1206 
1207 #ifdef LIBMESH_HAVE_CXX11
1208 static_assert(sizeof(PetscInt) == sizeof(numeric_index_type),
1209  "PETSc and libMesh integer sizes must match!");
1210 #endif
1211 
1212 
1213 inline
1215 {
1216  return reinterpret_cast<PetscInt *>(const_cast<numeric_index_type *>(p));
1217 }
1218 
1219 } // namespace libMesh
1220 
1221 
1222 #endif // #ifdef LIBMESH_HAVE_PETSC
1223 #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:39
virtual numeric_index_type n() const libmesh_override
double abs(double a)
Encapsulates the MPI_Comm object.
Definition: parallel.h:657
virtual void init(const numeric_index_type, const numeric_index_type, const bool=false, const ParallelType=AUTOMATIC)=0
Change the dimension of the vector to N.
void scale(MeshBase &mesh, const Real xs, const Real ys=0., const Real zs=0.)
Scales the mesh.
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:436
Vec _local_form
PETSc vector datatype to hold the local form of a ghosted vector.
Definition: petsc_vector.h:368
uint8_t processor_id_type
Definition: id_types.h:99
const class libmesh_nullptr_t libmesh_nullptr
Numeric vector.
Definition: dof_map.h:66
numeric_index_type map_global_to_local_index(const numeric_index_type i) const
IterBase * end
Also have a polymorphic pointer to the end object, this prevents iterating past the end...
bool _is_initialized
true once init() has been called.
virtual bool initialized() const
The libMesh namespace provides an interface to certain functionality in the library.
void sum(T &r, const Communicator &comm=Communicator_World)
std::mutex _petsc_vector_mutex
Mutex for _get_array and _restore_array.
Definition: petsc_vector.h:396
const Number zero
.
Definition: libmesh.h:178
long double max(long double a, double b)
numeric_index_type _first
First local index.
Definition: petsc_vector.h:349
libmesh_assert(j)
virtual bool closed() const libmesh_override
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
PetscInt * numeric_petsc_cast(const numeric_index_type *p)
const PetscScalar * _read_only_values
Pointer to the actual PETSc array of the values of the vector.
Definition: petsc_vector.h:378
std::vector< T > _values
Actual vector datatype to hold vector entries.
DataType _vec
Actual Eigen::SparseVector<> we are wrapping.
ParallelType
Defines an enum for parallel data structure types.
dof_id_type numeric_index_type
Definition: id_types.h:92
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:430
void _get_array(bool read_only) const
Queries the array (and the local form if the vector is ghosted) from PETSc.
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
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:424
PetscErrorCode Vec Mat libmesh_dbg_var(j)
virtual ElemType type() const libmesh_override
Definition: cell_hex20.h:81
Threads::spin_mutex _petsc_vector_mutex
Definition: petsc_vector.h:398
virtual void init() libmesh_override
Initialize this matrix using the sparsity structure computed by dof_map.
ParallelType _type
Type of vector.
numeric_index_type _last
Last local index.
Definition: petsc_vector.h:356
void _restore_array() const
Restores the array (and the local form if the vector is ghosted) to PETSc.
std::unordered_map< numeric_index_type, numeric_index_type > GlobalToLocalMap
Type for map that maps global to local ghost cells.
Definition: petsc_vector.h:418
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void clear() libmesh_override
Restores the NumericVector<T> to a pristine state.
virtual void swap(NumericVector< T > &v)
Swaps the contents of this with v.
PetscErrorCode ierr
X_input swap(X_system)
static PetscErrorCode Mat * A
const Parallel::Communicator & comm() const
void swap(Iterator &lhs, Iterator &rhs)
swap, used to implement op=
numeric_index_type _local_size
The local vector size.
virtual void zero() libmesh_override
Set all entries to zero.
static const bool value
Definition: xdr_io.C:108
bool _values_read_only
Whether or not the data array is for read only access.
Definition: petsc_vector.h:441
PetscScalar * _values
Pointer to the actual PETSc array of the values of the vector.
Definition: petsc_vector.h:388
Vec _vec
Actual PETSc vector datatype to hold vector entries.
Definition: petsc_vector.h:330
long double min(long double a, double b)
const Elem & get(const ElemType type_in)
bool _is_closed
Flag which tracks whether the vector&#39;s values are consistent on all processors after insertion or add...
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:339