21 #ifndef LIBMESH_PETSC_VECTOR_H 22 #define LIBMESH_PETSC_VECTOR_H 25 #include "libmesh/libmesh_config.h" 27 #ifdef LIBMESH_HAVE_PETSC 30 #include "libmesh/numeric_vector.h" 31 #include "libmesh/petsc_macro.h" 32 #include "libmesh/int_range.h" 33 #include "libmesh/libmesh_common.h" 34 #include "libmesh/petsc_solver_exception.h" 35 #include "libmesh/parallel_only.h" 36 #include "libmesh/enum_to_string.h" 40 # define LIBMESH_SAW_I 44 # undef I // Avoid complex.h contamination 51 #include <unordered_map> 56 #include <condition_variable> 62 template <
typename T>
class SparseMatrix;
109 const std::vector<numeric_index_type> & ghost,
139 virtual void close ()
override;
144 virtual void clear () noexcept
override;
146 virtual void zero ()
override;
148 virtual std::unique_ptr<NumericVector<T>>
zero_clone ()
const override;
150 virtual std::unique_ptr<NumericVector<T>>
clone ()
const override;
154 const bool fast=
false,
158 const bool fast=
false,
163 const std::vector<numeric_index_type> & ghost,
164 const bool fast =
false,
168 const bool fast =
false)
override;
176 virtual Real min ()
const override;
178 virtual Real max ()
const override;
180 virtual T
sum ()
const override;
207 virtual void get(
const std::vector<numeric_index_type> & index,
208 T * values)
const override;
243 const T
value)
override;
246 const T
value)
override;
248 virtual void add (
const T s)
override;
261 const std::vector<numeric_index_type> & dof_indices)
override;
284 virtual void insert (
const T * v,
285 const std::vector<numeric_index_type> & dof_indices)
override;
287 virtual void scale (
const T factor)
override;
293 virtual void abs()
override;
304 virtual void localize (std::vector<T> & v_local)
const override;
309 const std::vector<numeric_index_type> & send_list)
const override;
311 virtual void localize (std::vector<T> & v_local,
312 const std::vector<numeric_index_type> & indices)
const override;
316 const std::vector<numeric_index_type> & send_list)
override;
330 const std::vector<numeric_index_type> & rows,
331 bool supplying_global_rows =
true)
const override;
349 virtual std::unique_ptr<NumericVector<T>>
350 get_subvector(
const std::vector<numeric_index_type> & rows)
override;
353 const std::vector<numeric_index_type> & rows)
override;
367 #ifdef LIBMESH_HAVE_CXX11_THREAD 475 template <NormType N>
Real norm ()
const;
483 template <
typename T>
488 _array_is_present(false),
492 _local_form(nullptr),
493 _read_only_values(nullptr),
495 _global_to_local_map(),
496 _destroy_vec_on_exit(true),
497 _values_manually_retrieved(false),
498 _values_read_only(false)
505 template <
typename T>
512 _array_is_present(false),
516 _local_form(nullptr),
517 _read_only_values(nullptr),
519 _global_to_local_map(),
520 _destroy_vec_on_exit(true),
521 _values_manually_retrieved(false),
522 _values_read_only(false)
524 this->
init(n, n,
false, ptype);
529 template <
typename T>
537 _array_is_present(false),
541 _local_form(nullptr),
542 _read_only_values(nullptr),
544 _global_to_local_map(),
545 _destroy_vec_on_exit(true),
546 _values_manually_retrieved(false),
547 _values_read_only(false)
549 this->
init(n, n_local,
false, ptype);
554 template <
typename T>
559 const std::vector<numeric_index_type> & ghost,
563 _array_is_present(false),
567 _local_form(nullptr),
568 _read_only_values(nullptr),
570 _global_to_local_map(),
571 _destroy_vec_on_exit(true),
572 _values_manually_retrieved(false),
573 _values_read_only(false)
575 this->
init(n, n_local, ghost,
false, ptype);
582 template <
typename T>
588 _array_is_present(false),
592 _local_form(nullptr),
593 _read_only_values(nullptr),
595 _global_to_local_map(),
596 _destroy_vec_on_exit(false),
597 _values_manually_retrieved(false),
598 _values_read_only(false)
600 this->_is_closed =
true;
605 PetscInt petsc_local_size=0;
606 LibmeshPetscCall(VecGetLocalSize(_vec, &petsc_local_size));
610 LibmeshPetscCall(VecGetType(_vec, &ptype));
612 #if PETSC_RELEASE_GREATER_EQUALS(3, 21, 0) 614 if (std::strcmp(ptype, VECMPI) == 0)
617 LibmeshPetscCall(VecGhostGetGhostIS(_vec, &ghostis));
620 LibmeshPetscCall(VecGhostGetLocalForm(_vec, &localrep));
625 if (ghostis && localrep)
628 LibmeshPetscCall(ISGetSize(ghostis, &ghost_size));
630 const PetscInt * indices;
631 LibmeshPetscCall(ISGetIndices(ghostis, &indices));
634 _global_to_local_map[indices[i]] = i;
636 LibmeshPetscCall(ISRestoreIndices(ghostis, &indices));
639 else if (std::strcmp(ptype,VECSHARED) == 0)
641 if ((std::strcmp(ptype,VECSHARED) == 0) || (std::strcmp(ptype,VECMPI) == 0))
644 ISLocalToGlobalMapping mapping;
645 LibmeshPetscCall(VecGetLocalToGlobalMapping(_vec, &mapping));
648 LibmeshPetscCall(VecGhostGetLocalForm(_vec,&localrep));
652 if (mapping && localrep)
657 LibmeshPetscCall(ISLocalToGlobalMappingGetSize(mapping, &n));
660 const PetscInt * indices;
661 LibmeshPetscCall(ISLocalToGlobalMappingGetIndices(mapping,&indices));
664 _global_to_local_map[indices[i]] = i-my_local_size;
666 LibmeshPetscCall(ISLocalToGlobalMappingRestoreIndices(mapping, &indices));
671 LibmeshPetscCall(VecGhostRestoreLocalForm(_vec,&localrep));
680 template <
typename T>
689 template <
typename T>
696 parallel_object_only();
698 PetscInt petsc_n=
static_cast<PetscInt
>(n);
718 if (this->_type ==
SERIAL)
720 LibmeshPetscCallA(PETSC_COMM_SELF, VecCreate(PETSC_COMM_SELF, &_vec));
721 LibmeshPetscCallA(PETSC_COMM_SELF, VecSetSizes(_vec, petsc_n, petsc_n));
722 LibmeshPetscCallA(PETSC_COMM_SELF, VecSetFromOptions (_vec));
727 #ifdef LIBMESH_HAVE_MPI 728 PetscInt petsc_n_local=cast_int<PetscInt>(n_local);
729 libmesh_assert_less_equal (n_local, n);
731 LibmeshPetscCall(VecCreate(this->comm().
get(), &_vec));
732 LibmeshPetscCall(VecSetSizes(_vec, petsc_n_local, petsc_n));
734 libmesh_assert_equal_to (n_local, n);
735 LibmeshPetscCallA(PETSC_COMM_SELF, VecCreate(PETSC_COMM_SELF, &_vec));
736 LibmeshPetscCallA(PETSC_COMM_SELF, VecSetSizes(_vec, petsc_n, petsc_n));
738 LibmeshPetscCall(VecSetFromOptions(_vec));
744 this->_is_closed =
true;
753 template <
typename T>
759 this->
init(n,n,fast,ptype);
764 template <
typename T>
768 const std::vector<numeric_index_type> & ghost,
772 parallel_object_only();
774 PetscInt petsc_n=
static_cast<PetscInt
>(n);
775 PetscInt petsc_n_local=
static_cast<PetscInt
>(n_local);
776 PetscInt petsc_n_ghost=
static_cast<PetscInt
>(ghost.size());
791 PetscInt * petsc_ghost = ghost.empty() ? LIBMESH_PETSC_NULLPTR :
792 const_cast<PetscInt *
>(
reinterpret_cast<const PetscInt *
>(ghost.data()));
804 _global_to_local_map[ghost[i]] = i;
808 LibmeshPetscCall(VecCreateGhost (this->comm().
get(), petsc_n_local, petsc_n,
809 petsc_n_ghost, petsc_ghost,
816 LibmeshPetscCall(PetscObjectAppendOptionsPrefix((PetscObject)_vec,
"ghost_"));
818 LibmeshPetscCall(VecSetFromOptions (_vec));
821 this->_is_closed =
true;
829 template <
typename T>
834 parallel_object_only();
840 const PetscVector<T> & v = cast_ref<const PetscVector<T> &>(other);
853 this->_is_closed =
true;
856 this->_type = v.
_type;
859 LibmeshPetscCall(VecDuplicate (v.
_vec, &this->_vec));
867 template <
typename T>
871 parallel_object_only();
873 this->_restore_array();
875 VecAssemblyBeginEnd(this->comm(), _vec);
878 VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
880 this->_is_closed =
true;
885 template <
typename T>
889 exceptionless_parallel_object_only();
892 this->_restore_array();
894 if ((this->
initialized()) && (this->_destroy_vec_on_exit))
898 PetscErrorCode ierr = VecDestroy(&_vec);
900 libmesh_warning(
"Warning: VecDestroy returned a non-zero error code which we ignored.");
905 _global_to_local_map.clear();
910 template <
typename T>
914 parallel_object_only();
918 this->_restore_array();
923 LibmeshPetscCall(VecSet (_vec, z));
929 LibmeshPetscCall(VecGhostGetLocalForm (_vec,&loc_vec));
930 LibmeshPetscCall(VecSet (loc_vec, z));
931 LibmeshPetscCall(VecGhostRestoreLocalForm (_vec, &loc_vec));
937 template <
typename T>
941 std::unique_ptr<NumericVector<T>> cloned_vector =
942 std::make_unique<PetscVector<T>>(this->comm(), this->type());
943 cloned_vector->
init(*
this);
944 return cloned_vector;
949 template <
typename T>
953 std::unique_ptr<NumericVector<T>> cloned_vector =
954 std::make_unique<PetscVector<T>>(this->comm(), this->type());
955 cloned_vector->
init(*
this,
true);
956 *cloned_vector = *
this;
957 return cloned_vector;
962 template <
typename T>
968 PetscInt petsc_size=0;
973 LibmeshPetscCall(VecGetSize(_vec, &petsc_size));
980 template <
typename T>
986 PetscInt petsc_size=0;
988 LibmeshPetscCall(VecGetLocalSize(_vec, &petsc_size));
995 template <
typename T>
1003 if (_array_is_present)
1007 PetscInt petsc_first=0, petsc_last=0;
1008 LibmeshPetscCall(VecGetOwnershipRange (_vec, &petsc_first, &petsc_last));
1017 template <
typename T>
1025 if (_array_is_present)
1029 PetscInt petsc_first=0, petsc_last=0;
1030 LibmeshPetscCall(VecGetOwnershipRange (_vec, &petsc_first, &petsc_last));
1039 template <
typename T>
1048 if (_array_is_present)
1055 PetscInt petsc_first=0, petsc_last=0;
1056 LibmeshPetscCall(VecGetOwnershipRange (_vec, &petsc_first, &petsc_last));
1062 if ((i>=first) && (i<last))
1067 GlobalToLocalMap::const_iterator it = _global_to_local_map.find(i);
1069 const GlobalToLocalMap::const_iterator end = _global_to_local_map.end();
1072 std::ostringstream error_message;
1073 error_message <<
"No index " << i <<
" in ghosted vector.\n" 1074 <<
"Vector contains [" << first <<
',' << last <<
")\n";
1075 GlobalToLocalMap::const_iterator b = _global_to_local_map.begin();
1078 error_message <<
"And empty ghost array.\n";
1082 error_message <<
"And ghost array {" << b->first;
1083 for (++b; b != end; ++b)
1084 error_message <<
',' << b->first;
1085 error_message <<
"}\n";
1088 libmesh_error_msg(error_message.str());
1092 return it->second+last-first;
1097 template <
typename T>
1101 this->_get_array(
true);
1108 libmesh_assert_less (local_index, _local_size);
1112 return static_cast<T
>(_read_only_values[local_index]);
1117 template <
typename T>
1122 this->_get_array(
true);
1124 const std::size_t num = index.size();
1126 for (std::size_t i=0; i<num; i++)
1132 libmesh_assert_less (local_index, _local_size);
1135 values[i] =
static_cast<T
>(_read_only_values[local_index]);
1140 template <
typename T>
1145 _values_manually_retrieved =
true;
1151 template <
typename T>
1156 _values_manually_retrieved =
true;
1158 return _read_only_values;
1161 template <
typename T>
1167 _values_manually_retrieved =
false;
1171 template <
typename T>
1175 parallel_object_only();
1177 this->_restore_array();
1180 PetscReal returnval=0.;
1182 LibmeshPetscCall(VecMin (_vec, &index, &returnval));
1185 return static_cast<Real>(returnval);
1190 template <
typename T>
1194 parallel_object_only();
1196 this->_restore_array();
1199 PetscReal returnval=0.;
1201 LibmeshPetscCall(VecMax (_vec, &index, &returnval));
1204 return static_cast<Real>(returnval);
1209 template <
typename T>
1213 parallel_object_only();
1219 std::swap(_vec, v.
_vec);
1223 #ifdef LIBMESH_HAVE_CXX11_THREAD 1231 std::swap(_values, v.
_values);
1236 template <
typename T>
1242 return std::numeric_limits<PetscInt>::max();
1247 #ifdef LIBMESH_HAVE_CXX11 1249 "PETSc and libMesh integer sizes must match!");
1262 #endif // #ifdef LIBMESH_HAVE_PETSC 1263 #endif // LIBMESH_PETSC_VECTOR_H std::string name(const ElemQuality q)
This function returns a string containing some name for q.
virtual Real l2_norm() const override
PetscVector(const Parallel::Communicator &comm_in, const ParallelType type=AUTOMATIC)
Dummy-Constructor.
bool closed()
Checks that the library has been closed.
virtual void localize(std::vector< T > &v_local) const override
Creates a copy of the global vector in the local vector v_local.
This class provides a nice interface to PETSc's Vec object.
void add_vector_conjugate_transpose(const NumericVector< T > &v, const SparseMatrix< T > &A)
.
virtual void add_vector_transpose(const NumericVector< T > &v, const SparseMatrix< T > &A) override
Computes , i.e.
virtual void zero() override
Set all entries to zero.
bool _destroy_vec_on_exit
This boolean value should only be set to false for the constructor which takes a PETSc Vec object...
virtual numeric_index_type size() const override
virtual Real linfty_norm() const override
numeric_index_type _local_size
Size of the local values from _get_array()
virtual numeric_index_type first_local_index() const override
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
virtual void add(const numeric_index_type i, const T value) override
Adds value to the vector entry specified by i.
Vec _vec
Actual PETSc vector datatype to hold vector entries.
numeric_index_type _last
Last local index.
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
PetscScalar * _values
Pointer to the actual PETSc array of the values of the vector.
bool _values_read_only
Whether or not the data array is for read only access.
The libMesh namespace provides an interface to certain functionality in the library.
virtual void pointwise_divide(const NumericVector< T > &vec1, const NumericVector< T > &vec2) override
Computes (summation not implied) i.e.
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.
uint8_t processor_id_type
PetscInt * numeric_petsc_cast(const numeric_index_type *p)
virtual void swap(NumericVector< T > &v) override
Swaps the contents of this with v.
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.
const PetscScalar * _read_only_values
Pointer to the actual PETSc array of the values of the vector.
virtual void pointwise_mult(const NumericVector< T > &vec1, const NumericVector< T > &vec2) override
Computes (summation not implied) i.e.
virtual Real l1_norm() const override
Vec _local_form
PETSc vector datatype to hold the local form of a ghosted vector.
GlobalToLocalMap _global_to_local_map
Map that maps global to local ghost cells (will be empty if not in ghost cell mode).
virtual void abs() override
Sets for each entry in the vector.
dof_id_type numeric_index_type
bool _is_initialized
Flag that tells if init() has been called.
const PetscScalar * get_array_read() const
Get read only access to the raw PETSc Vector data array.
virtual NumericVector< T > & operator+=(const NumericVector< T > &v) override
Adds v to *this, .
virtual void conjugate() override
Negates the imaginary component of each entry in the vector.
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 std::size_t max_allowed_id() const override
ParallelType _type
Type of vector.
virtual void create_subvector(NumericVector< T > &subvector, const std::vector< numeric_index_type > &rows, bool supplying_global_rows=true) const override
Fills in subvector from this vector using the indices in rows.
virtual numeric_index_type last_local_index() const override
std::string enum_to_string(const T e)
virtual T sum() const override
std::mutex _petsc_get_restore_array_mutex
Mutex for _get_array and _restore_array.
void restore_array()
Restore the data array.
ParallelType type() const
virtual std::unique_ptr< NumericVector< T > > clone() const override
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.
virtual void restore_subvector(std::unique_ptr< NumericVector< T >> subvector, const std::vector< numeric_index_type > &rows) override
Restores a view into this vector using the indices in rows.
virtual void clear() noexcept override
clear() is called from the destructor, so it should not throw.
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.
virtual Real max() const override
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
bool initialized()
Checks that library initialization has been done.
numeric_index_type _first
First local index.
PetscVector< T > & operator=(const PetscVector< T > &v)
Copy assignment operator.
virtual void reciprocal() override
Computes the component-wise reciprocal, .
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
Computes the component-wise multiplication of this vector's entries by another's, ...
bool _values_manually_retrieved
Whether or not the data array has been manually retrieved using get_array() or get_array_read() ...
virtual void print_matlab(const std::string &name="") const override
Print the contents of the vector in Matlab's sparse matrix format.
virtual std::unique_ptr< NumericVector< T > > get_subvector(const std::vector< numeric_index_type > &rows) override
Creates a view into this vector using the indices in rows.
virtual void scale(const T factor) override
Scale each element of the vector by the given factor.
PetscScalar * get_array()
Get read/write access to the raw PETSc Vector data array.
virtual std::unique_ptr< NumericVector< T > > zero_clone() const override
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
const Elem & get(const ElemType type_in)
virtual Real min() const override
virtual void close() override
Calls the NumericVector's internal assembly routines, ensuring that the values are consistent across ...
virtual NumericVector< T > & operator/=(const NumericVector< T > &v) override
Computes the component-wise division of this vector's entries by another's, .
ParallelType
Defines an enum for parallel data structure types.
virtual NumericVector< T > & operator-=(const NumericVector< T > &v) override
Subtracts v from *this, .
virtual T dot(const NumericVector< T > &v) const override