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 
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 
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 
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 
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 
261  const SparseMatrix<T> & A);
262 
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 
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 LIBMESH_BEST_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
452  NumericVector<T>(comm_in, ptype),
453  _array_is_present(false),
454  _first(0),
455  _last(0),
459  _destroy_vec_on_exit(true),
461  _values_read_only(false)
462 {
463  this->_type = ptype;
464 }
465 
466 
467 
468 template <typename T>
469 inline
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
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
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
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
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
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 
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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);
1191 
1192 #ifdef LIBMESH_HAVE_CXX11_THREAD
1193  // Only truly atomic for v... but swap() doesn't really need to be thread safe!
1195 #else
1197 #endif
1198 
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)
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:451
virtual void zero() libmesh_override
Definition: petsc_vector.h:859
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)
Definition: petsc_vector.C:277
virtual NumericVector< T > & operator-=(const NumericVector< T > &v) libmesh_override
Definition: petsc_vector.C:125
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:399
numeric_index_type _local_size
Definition: petsc_vector.h:361
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:625
uint8_t processor_id_type
Definition: id_types.h:99
numeric_index_type _last
Definition: petsc_vector.h:356
Threads::spin_mutex _petsc_vector_mutex
Definition: petsc_vector.h:398
const class libmesh_nullptr_t libmesh_nullptr
virtual T dot(const NumericVector< T > &v) const libmesh_override
Definition: petsc_vector.C:465
PetscScalar * _values
Definition: petsc_vector.h:388
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:339
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:396
const PetscScalar * _read_only_values
Definition: petsc_vector.h:378
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
Definition: petsc_vector.h:896
GlobalToLocalMap _global_to_local_map
Definition: petsc_vector.h:424
virtual NumericVector< T > & operator/=(NumericVector< T > &v) libmesh_override
Definition: petsc_vector.C:426
dof_id_type numeric_index_type
Definition: id_types.h:92
virtual NumericVector< T > & operator=(const T s) libmesh_override
Definition: petsc_vector.C:509
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
Definition: petsc_vector.h:907
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:380
virtual numeric_index_type size() const libmesh_override
Definition: petsc_vector.h:919
PetscErrorCode Vec Mat libmesh_dbg_var(j)
virtual NumericVector< T > & operator+=(const NumericVector< T > &v) libmesh_override
Definition: petsc_vector.C:111
T indefinite_dot(const NumericVector< T > &v) const
Definition: petsc_vector.C:486
virtual void add_vector_transpose(const NumericVector< T > &v, const SparseMatrix< T > &A) libmesh_override
Definition: petsc_vector.C:249
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:439
virtual void close() libmesh_override
Definition: petsc_vector.h:807
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:928
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
Definition: petsc_vector.h:939
ParallelType type() const
numeric_index_type _first
Definition: petsc_vector.h:349
LIBMESH_BEST_UNORDERED_MAP< numeric_index_type, numeric_index_type > GlobalToLocalMap
Definition: petsc_vector.h:418
virtual void swap(NumericVector< T > &v) libmesh_override
virtual numeric_index_type first_local_index() const libmesh_override
Definition: petsc_vector.h:956
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
Definition: petsc_vector.h:835
virtual T sum() const libmesh_override
Definition: petsc_vector.C:41
virtual numeric_index_type last_local_index() const libmesh_override
Definition: petsc_vector.h:980
PetscScalar * get_array()
virtual T operator()(const numeric_index_type i) const libmesh_override
virtual Real min() const libmesh_override