petsc_vector.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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 
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)
Definition: elem_quality.C:42
virtual Real l2_norm() const override
Definition: petsc_vector.C:74
PetscVector(const Parallel::Communicator &comm_in, const ParallelType type=AUTOMATIC)
Definition: petsc_vector.h:454
bool closed()
Definition: libmesh.C:265
virtual void localize(std::vector< T > &v_local) const override
Definition: petsc_vector.C:926
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 void add_vector_transpose(const NumericVector< T > &v, const SparseMatrix< T > &A) override
Definition: petsc_vector.C:249
virtual void zero() override
Definition: petsc_vector.h:862
virtual numeric_index_type size() const override
Definition: petsc_vector.h:922
virtual Real linfty_norm() const override
Definition: petsc_vector.C:92
numeric_index_type _local_size
Definition: petsc_vector.h:364
virtual numeric_index_type first_local_index() const override
Definition: petsc_vector.h:959
virtual bool initialized() const
void _get_array(bool read_only) const
virtual numeric_index_type local_size() const override
Definition: petsc_vector.h:942
virtual void add(const numeric_index_type i, const T value) override
Definition: petsc_vector.C:180
uint8_t processor_id_type
Definition: id_types.h:99
numeric_index_type _last
Definition: petsc_vector.h:359
Threads::spin_mutex _petsc_vector_mutex
Definition: petsc_vector.h:401
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
Definition: diff_context.h:40
PetscScalar * _values
Definition: petsc_vector.h:391
IterBase * end
virtual void get(const std::vector< numeric_index_type > &index, T *values) const override
const Number zero
Definition: libmesh.h:239
virtual void init(const numeric_index_type n, const numeric_index_type n_local, const bool fast=false, const ParallelType ptype=AUTOMATIC)=0
void _restore_array() const
std::atomic< bool > _array_is_present
Definition: petsc_vector.h:342
PetscInt * numeric_petsc_cast(const numeric_index_type *p)
std::mutex _petsc_vector_mutex
Definition: petsc_vector.h:399
virtual void swap(NumericVector< T > &v) override
virtual void insert(const T *v, const std::vector< numeric_index_type > &dof_indices) override
Definition: petsc_vector.C:380
const PetscScalar * _read_only_values
Definition: petsc_vector.h:381
virtual void pointwise_mult(const NumericVector< T > &vec1, const NumericVector< T > &vec2) override
virtual Real l1_norm() const override
Definition: petsc_vector.C:57
GlobalToLocalMap _global_to_local_map
Definition: petsc_vector.h:427
virtual void clear() override
Definition: petsc_vector.h:838
virtual void abs() override
Definition: petsc_vector.C:440
dof_id_type numeric_index_type
Definition: id_types.h:92
const PetscScalar * get_array_read() const
virtual NumericVector< T > & operator+=(const NumericVector< T > &v) override
Definition: petsc_vector.C:111
virtual void conjugate() override
Definition: petsc_vector.C:168
void init(triangulateio &t)
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
numeric_index_type map_global_to_local_index(const numeric_index_type i) const
virtual void create_subvector(NumericVector< T > &subvector, const std::vector< numeric_index_type > &rows) const override
virtual numeric_index_type last_local_index() const override
Definition: petsc_vector.h:983
virtual T sum() const override
Definition: petsc_vector.C:41
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
Definition: petsc_vector.C:487
std::unordered_map< numeric_index_type, numeric_index_type > GlobalToLocalMap
Definition: petsc_vector.h:421
virtual void swap(NumericVector< T > &v)
PetscErrorCode ierr
virtual void init(const numeric_index_type N, const numeric_index_type n_local, const bool fast=false, const ParallelType type=AUTOMATIC) override
Definition: petsc_vector.h:628
static PetscErrorCode Mat * A
void swap(Iterator &lhs, Iterator &rhs)
virtual Real max() const override
static const bool value
Definition: xdr_io.C:109
bool initialized()
Definition: libmesh.C:258
numeric_index_type _first
Definition: petsc_vector.h:352
PetscVector< T > & operator=(const PetscVector< T > &v)
Definition: petsc_vector.C:560
virtual void reciprocal() override
Definition: petsc_vector.C:156
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices) override
Definition: petsc_vector.C:198
virtual void print_matlab(const std::string &name="") const override
virtual void scale(const T factor) override
Definition: petsc_vector.C:400
PetscScalar * get_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
Definition: petsc_vector.h:810
virtual NumericVector< T > & operator/=(const NumericVector< T > &v) override
Definition: petsc_vector.C:427
virtual NumericVector< T > & operator-=(const NumericVector< T > &v) override
Definition: petsc_vector.C:125
virtual T dot(const NumericVector< T > &v) const override
Definition: petsc_vector.C:466