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