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"
34 #include LIBMESH_INCLUDE_UNORDERED_MAP
35 
36 // Petsc include files.
37 #include <petscvec.h>
38 
39 // C++ includes
40 #include <cstddef>
41 #include <cstring>
42 #include <vector>
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 : 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 
123  virtual void close () libmesh_override;
124 
128  virtual void clear () libmesh_override;
129 
134  virtual void zero () libmesh_override;
135 
141  virtual UniquePtr<NumericVector<T> > zero_clone () const libmesh_override;
142 
146  virtual UniquePtr<NumericVector<T> > clone () const libmesh_override;
147 
160  virtual void init (const numeric_index_type N,
161  const numeric_index_type n_local,
162  const bool fast=false,
163  const ParallelType type=AUTOMATIC) libmesh_override;
164 
168  virtual void init (const numeric_index_type N,
169  const bool fast=false,
170  const ParallelType type=AUTOMATIC) libmesh_override;
171 
176  virtual void init (const numeric_index_type /*N*/,
177  const numeric_index_type /*n_local*/,
178  const std::vector<numeric_index_type> & /*ghost*/,
179  const bool /*fast*/ = false,
180  const ParallelType = AUTOMATIC) libmesh_override;
181 
186  virtual void init (const NumericVector<T> & other,
187  const bool fast = false) libmesh_override;
188 
192  virtual NumericVector<T> & operator= (const T s) libmesh_override;
193 
197  virtual NumericVector<T> & operator= (const NumericVector<T> & V) libmesh_override;
198 
203 
207  virtual NumericVector<T> & operator= (const std::vector<T> & v) libmesh_override;
208 
214  virtual Real min () const libmesh_override;
215 
221  virtual Real max () const libmesh_override;
222 
226  virtual T sum () const libmesh_override;
227 
232  virtual Real l1_norm () const libmesh_override;
233 
239  virtual Real l2_norm () const libmesh_override;
240 
246  virtual Real linfty_norm () const libmesh_override;
247 
255  virtual numeric_index_type size () const libmesh_override;
256 
261  virtual numeric_index_type local_size() const libmesh_override;
262 
267  virtual numeric_index_type first_local_index() const libmesh_override;
268 
273  virtual numeric_index_type last_local_index() const libmesh_override;
274 
282 
286  virtual T operator() (const numeric_index_type i) const libmesh_override;
287 
294  virtual void get(const std::vector<numeric_index_type> & index,
295  T * values) const libmesh_override;
296 
303  PetscScalar * get_array();
304 
311  const PetscScalar * get_array_read() const;
312 
319  void restore_array();
320 
325  virtual NumericVector<T> & operator += (const NumericVector<T> & V) libmesh_override;
326 
331  virtual NumericVector<T> & operator -= (const NumericVector<T> & V) libmesh_override;
332 
336  virtual void reciprocal() libmesh_override;
337 
342  virtual void conjugate() libmesh_override;
343 
347  virtual void set (const numeric_index_type i,
348  const T value) libmesh_override;
349 
353  virtual void add (const numeric_index_type i,
354  const T value) libmesh_override;
355 
361  virtual void add (const T s) libmesh_override;
362 
368  virtual void add (const NumericVector<T> & V) libmesh_override;
369 
375  virtual void add (const T a, const NumericVector<T> & v) libmesh_override;
376 
382 
387  virtual void add_vector (const T * v,
388  const std::vector<numeric_index_type> & dof_indices) libmesh_override;
389 
394  virtual void add_vector (const NumericVector<T> & V,
395  const SparseMatrix<T> & A) libmesh_override;
396 
402  virtual void add_vector_transpose (const NumericVector<T> & V,
403  const SparseMatrix<T> & A_trans) libmesh_override;
404 
411  const SparseMatrix<T> & A_trans);
412 
418 
423  virtual void insert (const T * v,
424  const std::vector<numeric_index_type> & dof_indices) libmesh_override;
425 
430  virtual void scale (const T factor) libmesh_override;
431 
435  virtual NumericVector<T> & operator /= (NumericVector<T> & v) libmesh_override;
436 
441  virtual void abs() libmesh_override;
442 
447  virtual T dot(const NumericVector<T> &) const libmesh_override;
448 
453  T indefinite_dot(const NumericVector<T> &) const;
454 
459  virtual void localize (std::vector<T> & v_local) const libmesh_override;
460 
465  virtual void localize (NumericVector<T> & v_local) const libmesh_override;
466 
472  virtual void localize (NumericVector<T> & v_local,
473  const std::vector<numeric_index_type> & send_list) const libmesh_override;
474 
479  virtual void localize (std::vector<T> & v_local,
480  const std::vector<numeric_index_type> & indices) const libmesh_override;
481 
486  virtual void localize (const numeric_index_type first_local_idx,
487  const numeric_index_type last_local_idx,
488  const std::vector<numeric_index_type> & send_list) libmesh_override;
489 
496  virtual void localize_to_one (std::vector<T> & v_local,
497  const processor_id_type proc_id=0) const libmesh_override;
498 
503  virtual void pointwise_mult (const NumericVector<T> & vec1,
504  const NumericVector<T> & vec2) libmesh_override;
505 
512  virtual void print_matlab(const std::string & name = "") const libmesh_override;
513 
518  virtual void create_subvector(NumericVector<T> & subvector,
519  const std::vector<numeric_index_type> & rows) const libmesh_override;
520 
524  virtual void swap (NumericVector<T> & v) libmesh_override;
525 
531  Vec vec () { libmesh_assert (_vec); return _vec; }
532 
533 
534 private:
535 
540  Vec _vec;
541 
547 #ifdef LIBMESH_HAVE_CXX11_THREAD
548  // Note: we can't use std::atomic_flag here because we need load and store operations
549  mutable std::atomic<bool> _array_is_present;
550 #else
551  mutable bool _array_is_present;
552 #endif
553 
560 
567 
572 
578  mutable Vec _local_form;
579 
588  mutable const PetscScalar * _read_only_values;
589 
598  mutable PetscScalar * _values;
599 
605 #ifdef LIBMESH_HAVE_CXX11_THREAD
606  mutable std::mutex _petsc_vector_mutex;
607 #else
609 #endif
610 
617  void _get_array(bool read_only) const;
618 
623  void _restore_array() const;
624 
628  typedef LIBMESH_BEST_UNORDERED_MAP<numeric_index_type,numeric_index_type> GlobalToLocalMap;
629 
634  GlobalToLocalMap _global_to_local_map;
635 
641 
647 
651  mutable bool _values_read_only;
652 };
653 
654 
655 /*----------------------- Inline functions ----------------------------------*/
656 
657 
658 
659 template <typename T>
660 inline
662  NumericVector<T>(comm_in, ptype),
663  _array_is_present(false),
664  _first(0),
665  _last(0),
669  _destroy_vec_on_exit(true),
671  _values_read_only(false)
672 {
673  this->_type = ptype;
674 }
675 
676 
677 
678 template <typename T>
679 inline
681  const numeric_index_type n,
682  const ParallelType ptype) :
683  NumericVector<T>(comm_in, ptype),
684  _array_is_present(false),
688  _destroy_vec_on_exit(true),
690  _values_read_only(false)
691 {
692  this->init(n, n, false, ptype);
693 }
694 
695 
696 
697 template <typename T>
698 inline
700  const numeric_index_type n,
701  const numeric_index_type n_local,
702  const ParallelType ptype) :
703  NumericVector<T>(comm_in, ptype),
704  _array_is_present(false),
708  _destroy_vec_on_exit(true),
710  _values_read_only(false)
711 {
712  this->init(n, n_local, false, ptype);
713 }
714 
715 
716 
717 template <typename T>
718 inline
720  const numeric_index_type n,
721  const numeric_index_type n_local,
722  const std::vector<numeric_index_type> & ghost,
723  const ParallelType ptype) :
724  NumericVector<T>(comm_in, ptype),
725  _array_is_present(false),
729  _destroy_vec_on_exit(true),
731  _values_read_only(false)
732 {
733  this->init(n, n_local, ghost, false, ptype);
734 }
735 
736 
737 
738 
739 
740 template <typename T>
741 inline
743  const Parallel::Communicator & comm_in) :
744  NumericVector<T>(comm_in, AUTOMATIC),
745  _array_is_present(false),
749  _destroy_vec_on_exit(false),
751  _values_read_only(false)
752 {
753  this->_vec = v;
754  this->_is_closed = true;
755  this->_is_initialized = true;
756 
757  /* We need to ask PETSc about the (local to global) ghost value
758  mapping and create the inverse mapping out of it. */
759  PetscErrorCode ierr=0;
760  PetscInt petsc_local_size=0;
761  ierr = VecGetLocalSize(_vec, &petsc_local_size);
762  LIBMESH_CHKERR(ierr);
763 
764  // Get the vector type from PETSc.
765  // As of Petsc 3.0.0, the VecType #define lost its const-ness, so we
766  // need to have it in the code
767 #if PETSC_VERSION_LESS_THAN(3,0,0) || !PETSC_VERSION_LESS_THAN(3,4,0)
768  // Pre-3.0 and petsc-dev (as of October 2012) use non-const versions
769  VecType ptype;
770 #else
771  const VecType ptype;
772 #endif
773  ierr = VecGetType(_vec, &ptype);
774  LIBMESH_CHKERR(ierr);
775 
776  if((std::strcmp(ptype,VECSHARED) == 0) || (std::strcmp(ptype,VECMPI) == 0))
777  {
778 #if PETSC_RELEASE_LESS_THAN(3,1,1)
779  ISLocalToGlobalMapping mapping = _vec->mapping;
780 #else
781  ISLocalToGlobalMapping mapping;
782  ierr = VecGetLocalToGlobalMapping(_vec, &mapping);
783  LIBMESH_CHKERR(ierr);
784 #endif
785 
786  // If is a sparsely stored vector, set up our new mapping
787  if (mapping)
788  {
789  const numeric_index_type my_local_size = static_cast<numeric_index_type>(petsc_local_size);
790  const numeric_index_type ghost_begin = static_cast<numeric_index_type>(petsc_local_size);
791 #if PETSC_RELEASE_LESS_THAN(3,4,0)
792  const numeric_index_type ghost_end = static_cast<numeric_index_type>(mapping->n);
793 #else
794  PetscInt n;
795  ierr = ISLocalToGlobalMappingGetSize(mapping, &n);
796  LIBMESH_CHKERR(ierr);
797  const numeric_index_type ghost_end = static_cast<numeric_index_type>(n);
798 #endif
799 #if PETSC_RELEASE_LESS_THAN(3,1,1)
800  const PetscInt * indices = mapping->indices;
801 #else
802  const PetscInt * indices;
803  ierr = ISLocalToGlobalMappingGetIndices(mapping,&indices);
804  LIBMESH_CHKERR(ierr);
805 #endif
806  for(numeric_index_type i=ghost_begin; i<ghost_end; i++)
807  _global_to_local_map[indices[i]] = i-my_local_size;
808  this->_type = GHOSTED;
809 #if !PETSC_RELEASE_LESS_THAN(3,1,1)
810  ierr = ISLocalToGlobalMappingRestoreIndices(mapping, &indices);
811  LIBMESH_CHKERR(ierr);
812 #endif
813  }
814  else
815  this->_type = PARALLEL;
816  }
817  else
818  this->_type = SERIAL;
819 }
820 
821 
822 
823 
824 template <typename T>
825 inline
827 {
828  this->clear ();
829 }
830 
831 
832 
833 template <typename T>
834 inline
836  const numeric_index_type n_local,
837  const bool fast,
838  const ParallelType ptype)
839 {
840  parallel_object_only();
841 
842  PetscErrorCode ierr=0;
843  PetscInt petsc_n=static_cast<PetscInt>(n);
844  PetscInt petsc_n_local=static_cast<PetscInt>(n_local);
845 
846 
847  // Clear initialized vectors
848  if (this->initialized())
849  this->clear();
850 
851  if (ptype == AUTOMATIC)
852  {
853  if (n == n_local)
854  this->_type = SERIAL;
855  else
856  this->_type = PARALLEL;
857  }
858  else
859  this->_type = ptype;
860 
861  libmesh_assert ((this->_type==SERIAL && n==n_local) ||
862  this->_type==PARALLEL);
863 
864  // create a sequential vector if on only 1 processor
865  if (this->_type == SERIAL)
866  {
867  ierr = VecCreateSeq (PETSC_COMM_SELF, petsc_n, &_vec);
868  CHKERRABORT(PETSC_COMM_SELF,ierr);
869 
870  ierr = VecSetFromOptions (_vec);
871  CHKERRABORT(PETSC_COMM_SELF,ierr);
872  }
873  // otherwise create an MPI-enabled vector
874  else if (this->_type == PARALLEL)
875  {
876 #ifdef LIBMESH_HAVE_MPI
877  libmesh_assert_less_equal (n_local, n);
878  ierr = VecCreateMPI (this->comm().get(), petsc_n_local, petsc_n,
879  &_vec);
880  LIBMESH_CHKERR(ierr);
881 #else
882  libmesh_assert_equal_to (n_local, n);
883  ierr = VecCreateSeq (PETSC_COMM_SELF, petsc_n, &_vec);
884  CHKERRABORT(PETSC_COMM_SELF,ierr);
885 #endif
886 
887  ierr = VecSetFromOptions (_vec);
888  LIBMESH_CHKERR(ierr);
889  }
890  else
891  libmesh_error_msg("Unsupported type " << this->_type);
892 
893  this->_is_initialized = true;
894  this->_is_closed = true;
895 
896 
897  if (fast == false)
898  this->zero ();
899 }
900 
901 
902 
903 template <typename T>
904 inline
906  const bool fast,
907  const ParallelType ptype)
908 {
909  this->init(n,n,fast,ptype);
910 }
911 
912 
913 
914 template <typename T>
915 inline
917  const numeric_index_type n_local,
918  const std::vector<numeric_index_type> & ghost,
919  const bool fast,
920  const ParallelType libmesh_dbg_var(ptype))
921 {
922  parallel_object_only();
923 
924  PetscErrorCode ierr=0;
925  PetscInt petsc_n=static_cast<PetscInt>(n);
926  PetscInt petsc_n_local=static_cast<PetscInt>(n_local);
927  PetscInt petsc_n_ghost=static_cast<PetscInt>(ghost.size());
928 
929  // If the mesh is not disjoint, every processor will either have
930  // all the dofs, none of the dofs, or some non-zero dofs at the
931  // boundary between processors.
932  //
933  // However we can't assert this, because someone might want to
934  // construct a GHOSTED vector which doesn't include neighbor element
935  // dofs. Boyce tried to do so in user code, and we're going to want
936  // to do so in System::project_vector().
937  //
938  // libmesh_assert(n_local == 0 || n_local == n || !ghost.empty());
939 
940  libmesh_assert(sizeof(PetscInt) == sizeof(numeric_index_type));
941 
942  PetscInt * petsc_ghost = ghost.empty() ? PETSC_NULL :
943  const_cast<PetscInt *>(reinterpret_cast<const PetscInt *>(&ghost[0]));
944 
945  // Clear initialized vectors
946  if (this->initialized())
947  this->clear();
948 
949  libmesh_assert(ptype == AUTOMATIC || ptype == GHOSTED);
950  this->_type = GHOSTED;
951 
952  /* Make the global-to-local ghost cell map. */
953  for(numeric_index_type i=0; i<ghost.size(); i++)
954  {
955  _global_to_local_map[ghost[i]] = i;
956  }
957 
958  /* Create vector. */
959  ierr = VecCreateGhost (this->comm().get(), petsc_n_local, petsc_n,
960  petsc_n_ghost, petsc_ghost,
961  &_vec);
962  LIBMESH_CHKERR(ierr);
963 
964  ierr = VecSetFromOptions (_vec);
965  LIBMESH_CHKERR(ierr);
966 
967  this->_is_initialized = true;
968  this->_is_closed = true;
969 
970  if (fast == false)
971  this->zero ();
972 }
973 
974 
975 
976 template <typename T>
977 inline
979  const bool fast)
980 {
981  parallel_object_only();
982 
983  // Clear initialized vectors
984  if (this->initialized())
985  this->clear();
986 
987  const PetscVector<T> & v = cast_ref<const PetscVector<T> &>(other);
988 
989  // Other vector should restore array.
990  if(v.initialized())
991  {
992  v._restore_array();
993  }
994 
996 
997  // Even if we're initializeing sizes based on an uninitialized or
998  // unclosed vector, *this* vector is being initialized now and is
999  // initially closed.
1000  this->_is_closed = true; // v._is_closed;
1001  this->_is_initialized = true; // v._is_initialized;
1002 
1003  this->_type = v._type;
1004 
1005  // We want to have a valid Vec, even if it's initially of size zero
1006  PetscErrorCode ierr = VecDuplicate (v._vec, &this->_vec);
1007  LIBMESH_CHKERR(ierr);
1008 
1009  if (fast == false)
1010  this->zero ();
1011 }
1012 
1013 
1014 
1015 template <typename T>
1016 inline
1018 {
1019  parallel_object_only();
1020 
1021  this->_restore_array();
1022 
1023  PetscErrorCode ierr=0;
1024 
1025  ierr = VecAssemblyBegin(_vec);
1026  LIBMESH_CHKERR(ierr);
1027  ierr = VecAssemblyEnd(_vec);
1028  LIBMESH_CHKERR(ierr);
1029 
1030  if(this->type() == GHOSTED)
1031  {
1032  ierr = VecGhostUpdateBegin(_vec,INSERT_VALUES,SCATTER_FORWARD);
1033  LIBMESH_CHKERR(ierr);
1034  ierr = VecGhostUpdateEnd(_vec,INSERT_VALUES,SCATTER_FORWARD);
1035  LIBMESH_CHKERR(ierr);
1036  }
1037 
1038  this->_is_closed = true;
1039 }
1040 
1041 
1042 
1043 template <typename T>
1044 inline
1046 {
1047  parallel_object_only();
1048 
1049  if (this->initialized())
1050  this->_restore_array();
1051 
1052  if ((this->initialized()) && (this->_destroy_vec_on_exit))
1053  {
1054  PetscErrorCode ierr=0;
1055 
1056  ierr = LibMeshVecDestroy(&_vec);
1057  LIBMESH_CHKERR(ierr);
1058  }
1059 
1060  this->_is_closed = this->_is_initialized = false;
1061 
1062  _global_to_local_map.clear();
1063 }
1064 
1065 
1066 
1067 template <typename T>
1068 inline
1070 {
1071  parallel_object_only();
1072 
1073  libmesh_assert(this->closed());
1074 
1075  this->_restore_array();
1076 
1077  PetscErrorCode ierr=0;
1078 
1079  PetscScalar z=0.;
1080 
1081  if(this->type() != GHOSTED)
1082  {
1083  ierr = VecSet (_vec, z);
1084  LIBMESH_CHKERR(ierr);
1085  }
1086  else
1087  {
1088  /* Vectors that include ghost values require a special
1089  handling. */
1090  Vec loc_vec;
1091  ierr = VecGhostGetLocalForm (_vec,&loc_vec);
1092  LIBMESH_CHKERR(ierr);
1093 
1094  ierr = VecSet (loc_vec, z);
1095  LIBMESH_CHKERR(ierr);
1096 
1097  ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
1098  LIBMESH_CHKERR(ierr);
1099  }
1100 }
1101 
1102 
1103 
1104 template <typename T>
1105 inline
1107 {
1108  NumericVector<T> * cloned_vector = new PetscVector<T>(this->comm(), this->type());
1109  cloned_vector->init(*this);
1110  return UniquePtr<NumericVector<T> >(cloned_vector);
1111 }
1112 
1113 
1114 
1115 template <typename T>
1116 inline
1118 {
1119  NumericVector<T> * cloned_vector = new PetscVector<T>(this->comm(), this->type());
1120  cloned_vector->init(*this, true);
1121  *cloned_vector = *this;
1122  return UniquePtr<NumericVector<T> >(cloned_vector);
1123 }
1124 
1125 
1126 
1127 template <typename T>
1128 inline
1130 {
1131  libmesh_assert (this->initialized());
1132 
1133  PetscErrorCode ierr=0;
1134  PetscInt petsc_size=0;
1135 
1136  if (!this->initialized())
1137  return 0;
1138 
1139  ierr = VecGetSize(_vec, &petsc_size);
1140  LIBMESH_CHKERR(ierr);
1141 
1142  return static_cast<numeric_index_type>(petsc_size);
1143 }
1144 
1145 
1146 
1147 template <typename T>
1148 inline
1150 {
1151  libmesh_assert (this->initialized());
1152 
1153  PetscErrorCode ierr=0;
1154  PetscInt petsc_size=0;
1155 
1156  ierr = VecGetLocalSize(_vec, &petsc_size);
1157  LIBMESH_CHKERR(ierr);
1158 
1159  return static_cast<numeric_index_type>(petsc_size);
1160 }
1161 
1162 
1163 
1164 template <typename T>
1165 inline
1167 {
1168  libmesh_assert (this->initialized());
1169 
1170  numeric_index_type first = 0;
1171 
1172  if(_array_is_present) // Can we use cached values?
1173  first = _first;
1174  else
1175  {
1176  PetscErrorCode ierr=0;
1177  PetscInt petsc_first=0, petsc_last=0;
1178  ierr = VecGetOwnershipRange (_vec, &petsc_first, &petsc_last);
1179  LIBMESH_CHKERR(ierr);
1180  first = static_cast<numeric_index_type>(petsc_first);
1181  }
1182 
1183  return first;
1184 }
1185 
1186 
1187 
1188 template <typename T>
1189 inline
1191 {
1192  libmesh_assert (this->initialized());
1193 
1194  numeric_index_type last = 0;
1195 
1196  if(_array_is_present) // Can we use cached values?
1197  last = _last;
1198  else
1199  {
1200  PetscErrorCode ierr=0;
1201  PetscInt petsc_first=0, petsc_last=0;
1202  ierr = VecGetOwnershipRange (_vec, &petsc_first, &petsc_last);
1203  LIBMESH_CHKERR(ierr);
1204  last = static_cast<numeric_index_type>(petsc_last);
1205  }
1206 
1207  return last;
1208 }
1209 
1210 
1211 
1212 template <typename T>
1213 inline
1215 {
1216  libmesh_assert (this->initialized());
1217 
1218  numeric_index_type first=0;
1219  numeric_index_type last=0;
1220 
1221  if(_array_is_present) // Can we use cached values?
1222  {
1223  first = _first;
1224  last = _last;
1225  }
1226  else
1227  {
1228  PetscErrorCode ierr=0;
1229  PetscInt petsc_first=0, petsc_last=0;
1230  ierr = VecGetOwnershipRange (_vec, &petsc_first, &petsc_last);
1231  LIBMESH_CHKERR(ierr);
1232  first = static_cast<numeric_index_type>(petsc_first);
1233  last = static_cast<numeric_index_type>(petsc_last);
1234  }
1235 
1236 
1237  if((i>=first) && (i<last))
1238  {
1239  return i-first;
1240  }
1241 
1242  GlobalToLocalMap::const_iterator it = _global_to_local_map.find(i);
1243 #ifndef NDEBUG
1244  const GlobalToLocalMap::const_iterator end = _global_to_local_map.end();
1245  if (it == end)
1246  {
1247  std::ostringstream error_message;
1248  error_message << "No index " << i << " in ghosted vector.\n"
1249  << "Vector contains [" << first << ',' << last << ")\n";
1250  GlobalToLocalMap::const_iterator b = _global_to_local_map.begin();
1251  if (b == end)
1252  {
1253  error_message << "And empty ghost array.\n";
1254  }
1255  else
1256  {
1257  error_message << "And ghost array {" << b->first;
1258  for (++b; b != end; ++b)
1259  error_message << ',' << b->first;
1260  error_message << "}\n";
1261  }
1262 
1263  libmesh_error_msg(error_message.str());
1264  }
1265  libmesh_assert (it != _global_to_local_map.end());
1266 #endif
1267  return it->second+last-first;
1268 }
1269 
1270 
1271 
1272 template <typename T>
1273 inline
1275 {
1276  this->_get_array(true);
1277 
1278  const numeric_index_type local_index = this->map_global_to_local_index(i);
1279 
1280 #ifndef NDEBUG
1281  if(this->type() == GHOSTED)
1282  {
1283  libmesh_assert_less (local_index, _local_size);
1284  }
1285 #endif
1286 
1287  return static_cast<T>(_read_only_values[local_index]);
1288 }
1289 
1290 
1291 
1292 template <typename T>
1293 inline
1294 void PetscVector<T>::get(const std::vector<numeric_index_type> & index,
1295  T * values) const
1296 {
1297  this->_get_array(true);
1298 
1299  const std::size_t num = index.size();
1300 
1301  for(std::size_t i=0; i<num; i++)
1302  {
1303  const numeric_index_type local_index = this->map_global_to_local_index(index[i]);
1304 #ifndef NDEBUG
1305  if(this->type() == GHOSTED)
1306  {
1307  libmesh_assert_less (local_index, _local_size);
1308  }
1309 #endif
1310  values[i] = static_cast<T>(_read_only_values[local_index]);
1311  }
1312 }
1313 
1314 
1315 template <typename T>
1316 inline
1318 {
1320  _get_array(false);
1321 
1322  return _values;
1323 }
1324 
1325 
1326 template <typename T>
1327 inline
1328 const PetscScalar * PetscVector<T>::get_array_read() const
1329 {
1331  _get_array(true);
1332 
1333  return _read_only_values;
1334 }
1335 
1336 template <typename T>
1337 inline
1339 {
1340  // Note _values_manually_retrieved needs to be set to false
1341  // BEFORE calling _restore_array()!
1343  _restore_array();
1344 }
1345 
1346 template <typename T>
1347 inline
1349 {
1350  parallel_object_only();
1351 
1352  this->_restore_array();
1353 
1354  PetscErrorCode ierr=0;
1355  PetscInt index=0;
1356  PetscReal returnval=0.;
1357 
1358  ierr = VecMin (_vec, &index, &returnval);
1359  LIBMESH_CHKERR(ierr);
1360 
1361  // this return value is correct: VecMin returns a PetscReal
1362  return static_cast<Real>(returnval);
1363 }
1364 
1365 
1366 
1367 template <typename T>
1368 inline
1370 {
1371  parallel_object_only();
1372 
1373  this->_restore_array();
1374 
1375  PetscErrorCode ierr=0;
1376  PetscInt index=0;
1377  PetscReal returnval=0.;
1378 
1379  ierr = VecMax (_vec, &index, &returnval);
1380  LIBMESH_CHKERR(ierr);
1381 
1382  // this return value is correct: VecMax returns a PetscReal
1383  return static_cast<Real>(returnval);
1384 }
1385 
1386 
1387 
1388 template <typename T>
1389 inline
1391 {
1392  parallel_object_only();
1393 
1394  NumericVector<T>::swap(other);
1395 
1396  PetscVector<T> & v = cast_ref<PetscVector<T> &>(other);
1397 
1398  std::swap(_vec, v._vec);
1401 
1402 #ifdef LIBMESH_HAVE_CXX11_THREAD
1403  // Only truly atomic for v... but swap() doesn't really need to be thread safe!
1405 #else
1407 #endif
1408 
1411 }
1412 
1413 
1414 
1415 
1416 
1417 #ifdef LIBMESH_HAVE_CXX11
1418 static_assert(sizeof(PetscInt) == sizeof(numeric_index_type),
1419  "PETSc and libMesh integer sizes must match!");
1420 #endif
1421 
1422 
1423 inline
1425 {
1426  return reinterpret_cast<PetscInt *>(const_cast<numeric_index_type *>(p));
1427 }
1428 
1429 } // namespace libMesh
1430 
1431 
1432 #endif // #ifdef LIBMESH_HAVE_PETSC
1433 #endif // LIBMESH_PETSC_VECTOR_H
T indefinite_dot(const NumericVector< T > &) const
Definition: petsc_vector.C:465
std::string name(const ElemQuality q)
Definition: elem_quality.C:39
void _get_array(bool read_only) const
virtual bool closed() const
PetscVector(const Parallel::Communicator &comm_in, const ParallelType type=AUTOMATIC)
Definition: petsc_vector.h:661
virtual void zero() libmesh_override
virtual void pointwise_mult(const NumericVector< T > &vec1, const NumericVector< T > &vec2) libmesh_override
NumericVector interface to PETSc Vec.
Definition: petsc_vector.h:64
void add_vector_conjugate_transpose(const NumericVector< T > &V, const SparseMatrix< T > &A_trans)
Definition: petsc_vector.C:263
virtual void init(const numeric_index_type, const numeric_index_type, const bool=false, const ParallelType=AUTOMATIC)=0
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices) libmesh_override
Definition: petsc_vector.C:198
virtual void scale(const T factor) libmesh_override
Definition: petsc_vector.C:378
numeric_index_type _local_size
Definition: petsc_vector.h:571
virtual void init(const numeric_index_type N, const numeric_index_type n_local, const bool fast=false, const ParallelType type=AUTOMATIC) libmesh_override
Definition: petsc_vector.h:835
virtual NumericVector< T > & operator+=(const NumericVector< T > &V) libmesh_override
Definition: petsc_vector.C:111
uint8_t processor_id_type
Definition: id_types.h:99
numeric_index_type _last
Definition: petsc_vector.h:566
Threads::spin_mutex _petsc_vector_mutex
Definition: petsc_vector.h:608
const class libmesh_nullptr_t libmesh_nullptr
PetscScalar * _values
Definition: petsc_vector.h:598
IterBase * end
virtual bool initialized() const
virtual Real l1_norm() const libmesh_override
Definition: petsc_vector.C:57
std::atomic< bool > _array_is_present
Definition: petsc_vector.h:549
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
PetscInt * numeric_petsc_cast(const numeric_index_type *p)
std::mutex _petsc_vector_mutex
Definition: petsc_vector.h:606
const PetscScalar * _read_only_values
Definition: petsc_vector.h:588
virtual Real linfty_norm() const libmesh_override
Definition: petsc_vector.C:92
numeric_index_type map_global_to_local_index(const numeric_index_type i) const
virtual UniquePtr< NumericVector< T > > zero_clone() const libmesh_override
virtual T dot(const NumericVector< T > &) const libmesh_override
Definition: petsc_vector.C:444
GlobalToLocalMap _global_to_local_map
Definition: petsc_vector.h:634
virtual NumericVector< T > & operator/=(NumericVector< T > &v) libmesh_override
Definition: petsc_vector.C:405
dof_id_type numeric_index_type
Definition: id_types.h:92
virtual NumericVector< T > & operator=(const T s) libmesh_override
Definition: petsc_vector.C:488
virtual void localize_to_one(std::vector< T > &v_local, const processor_id_type proc_id=0) const libmesh_override
void _restore_array() const
virtual UniquePtr< NumericVector< T > > clone() const libmesh_override
const PetscScalar * get_array_read() const
virtual void insert(const T *v, const std::vector< numeric_index_type > &dof_indices) libmesh_override
Definition: petsc_vector.C:359
virtual numeric_index_type size() const libmesh_override
PetscErrorCode Vec Mat libmesh_dbg_var(j)
virtual void create_subvector(NumericVector< T > &subvector, const std::vector< numeric_index_type > &rows) const libmesh_override
virtual void get(const std::vector< numeric_index_type > &index, T *values) const libmesh_override
virtual void abs() libmesh_override
Definition: petsc_vector.C:418
virtual void close() libmesh_override
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void swap(NumericVector< T > &v)
PetscErrorCode ierr
virtual void localize(std::vector< T > &v_local) const libmesh_override
Definition: petsc_vector.C:907
static PetscErrorCode Mat * A
virtual void add(const numeric_index_type i, const T value) libmesh_override
Definition: petsc_vector.C:180
virtual void print_matlab(const std::string &name="") const libmesh_override
const Parallel::Communicator & comm() const
void swap(Iterator &lhs, Iterator &rhs)
virtual Real l2_norm() const libmesh_override
Definition: petsc_vector.C:74
virtual numeric_index_type local_size() const libmesh_override
ParallelType type() const
virtual NumericVector< T > & operator-=(const NumericVector< T > &V) libmesh_override
Definition: petsc_vector.C:125
numeric_index_type _first
Definition: petsc_vector.h:559
virtual void add_vector_transpose(const NumericVector< T > &V, const SparseMatrix< T > &A_trans) libmesh_override
Definition: petsc_vector.C:242
LIBMESH_BEST_UNORDERED_MAP< numeric_index_type, numeric_index_type > GlobalToLocalMap
Definition: petsc_vector.h:628
virtual void swap(NumericVector< T > &v) libmesh_override
virtual numeric_index_type first_local_index() const libmesh_override
virtual Real max() const libmesh_override
virtual void reciprocal() libmesh_override
Definition: petsc_vector.C:156
virtual void conjugate() libmesh_override
Definition: petsc_vector.C:168
virtual void clear() libmesh_override
virtual T sum() const libmesh_override
Definition: petsc_vector.C:41
virtual numeric_index_type last_local_index() const libmesh_override
PetscScalar * get_array()
virtual T operator()(const numeric_index_type i) const libmesh_override
virtual Real min() const libmesh_override