16 #include <deal.II/lac/petsc_vector_base.h> 18 #ifdef DEAL_II_WITH_PETSC 20 # include <deal.II/base/memory_consumption.h> 21 # include <deal.II/lac/exceptions.h> 22 # include <deal.II/lac/petsc_compatibility.h> 23 # include <deal.II/lac/petsc_vector.h> 24 # include <deal.II/lac/petsc_parallel_vector.h> 26 # include <deal.II/base/multithread_info.h> 28 DEAL_II_NAMESPACE_OPEN
34 VectorReference::operator PetscScalar ()
const 36 Assert (index < vector.size(),
42 if (dynamic_cast<const PETScWrappers::Vector *>(&vector) != 0)
46 PetscErrorCode ierr = VecGetValues(vector.vector, 1, &idx, &value);
52 else if (dynamic_cast<const PETScWrappers::MPI::Vector *>(&vector) != 0)
74 PetscErrorCode ierr = VecGetOwnershipRange (vector.vector, &begin,
78 Vec locally_stored_elements = PETSC_NULL;
79 ierr = VecGhostGetLocalForm(vector.vector, &locally_stored_elements);
83 ierr = VecGetSize(locally_stored_elements, &lsize);
87 ierr = VecGetArray(locally_stored_elements, &ptr);
92 if ( index>=static_cast<size_type>(begin)
93 && index<static_cast<size_type>(end) )
96 value = *(ptr+index-begin);
102 = vector.ghost_indices.index_within_set(index);
105 value = *(ptr+ghostidx+end-begin);
108 ierr = VecRestoreArray(locally_stored_elements, &ptr);
111 ierr = VecGhostRestoreLocalForm(vector.vector, &locally_stored_elements);
123 PetscErrorCode ierr = VecGetOwnershipRange (vector.vector, &begin, &end);
128 AssertThrow ((index >= static_cast<size_type>(begin)) &&
129 (index < static_cast<size_type>(end)),
130 ExcAccessToNonlocalElement (index, begin, end-1));
134 PetscInt idx = index;
136 ierr = VecGetValues(vector.vector, 1, &idx, &value);
155 attained_ownership(true)
158 ExcMessage(
"PETSc does not support multi-threaded access, set " 159 "the thread limit to 1 in MPI_InitFinalize()."));
168 ghost_indices(v.ghost_indices),
170 attained_ownership(true)
173 ExcMessage(
"PETSc does not support multi-threaded access, set " 174 "the thread limit to 1 in MPI_InitFinalize()."));
191 attained_ownership(false)
194 ExcMessage(
"PETSc does not support multi-threaded access, set " 195 "the thread limit to 1 in MPI_InitFinalize()."));
204 #if DEAL_II_PETSC_VERSION_LT(3,2,0) 205 const PetscErrorCode ierr = VecDestroy (
vector);
207 const PetscErrorCode ierr = VecDestroy (&
vector);
221 #if DEAL_II_PETSC_VERSION_LT(3,2,0) 222 const PetscErrorCode ierr = VecDestroy (
vector);
224 const PetscErrorCode ierr = VecDestroy (&
vector);
244 PetscErrorCode ierr = VecSet (
vector, s);
249 Vec ghost = PETSC_NULL;
250 ierr = VecGhostGetLocalForm(
vector, &ghost);
253 ierr = VecSet (ghost, s);
256 ierr = VecGhostRestoreLocalForm(
vector, &ghost);
271 PetscBooleanType flag;
272 const PetscErrorCode ierr = VecEqual (
vector, v.
vector, &flag);
275 return (flag == PETSC_TRUE);
286 PetscBooleanType flag;
287 const PetscErrorCode ierr = VecEqual (
vector, v.
vector, &flag);
290 return (flag == PETSC_FALSE);
295 VectorBase::size_type
299 const PetscErrorCode ierr = VecGetSize (
vector, &sz);
307 VectorBase::size_type
311 const PetscErrorCode ierr = VecGetLocalSize (
vector, &sz);
319 std::pair<VectorBase::size_type, VectorBase::size_type>
323 const PetscErrorCode ierr = VecGetOwnershipRange (static_cast<const Vec &>(
vector),
327 return std::make_pair (begin, end);
334 const std::vector<PetscScalar> &values)
336 Assert (indices.size() == values.size(),
337 ExcMessage (
"Function called with arguments of different sizes"));
345 const std::vector<PetscScalar> &values)
347 Assert (indices.size() == values.size(),
348 ExcMessage (
"Function called with arguments of different sizes"));
356 const ::Vector<PetscScalar> &values)
358 Assert (indices.size() == values.size(),
359 ExcMessage (
"Function called with arguments of different sizes"));
367 const size_type *indices,
368 const PetscScalar *values)
387 const PetscErrorCode ierr = VecDot (vec.
vector,
vector, &result);
411 #ifdef DEAL_II_WITH_MPI 415 int all_int_last_action;
417 const int ierr = MPI_Allreduce
418 (&my_int_last_action, &all_int_last_action, 1, MPI_INT, MPI_BOR,
424 ExcMessage(
"Error: not all processors agree on the last " 425 "VectorOperation before this compress() call."));
432 ExcMessage(
"Missing compress() or calling with wrong VectorOperation argument."));
448 PetscErrorCode ierr = VecAssemblyBegin(
vector);
450 ierr = VecAssemblyEnd(
vector);
461 VectorBase::real_type
475 if (dynamic_cast<const PETScWrappers::MPI::Vector *>(
this) != 0)
478 const PetscErrorCode ierr = VecSum(
vector, &sum);
480 return sum/
static_cast<PetscReal
>(
size());
485 PetscScalar *start_ptr;
486 PetscErrorCode ierr = VecGetArray (
vector, &start_ptr);
489 PetscScalar mean = 0;
491 PetscScalar sum0 = 0,
499 const PetscScalar *ptr = start_ptr;
500 const PetscScalar *eptr = ptr + (
size()/4)*4;
509 while (ptr != start_ptr+
size())
512 mean = (sum0+sum1+sum2+sum3)/static_cast<PetscReal>(
size());
517 ierr = VecRestoreArray (
vector, &start_ptr);
524 VectorBase::real_type
529 const PetscErrorCode ierr = VecNorm (
vector, NORM_1, &d);
537 VectorBase::real_type
542 const PetscErrorCode ierr = VecNorm (
vector, NORM_2, &d);
550 VectorBase::real_type
555 PetscScalar *start_ptr;
556 PetscErrorCode ierr = VecGetArray (
vector, &start_ptr);
569 const PetscScalar *ptr = start_ptr;
570 const PetscScalar *eptr = ptr + (
size()/4)*4;
579 while (ptr != start_ptr+
size())
582 norm = std::pow(sum0+sum1+sum2+sum3,
588 ierr = VecRestoreArray (
vector, &start_ptr);
596 VectorBase::real_type
601 const PetscErrorCode ierr = VecNorm (
vector, NORM_INFINITY, &d);
609 VectorBase::real_type
614 const PetscErrorCode ierr = VecNormalize (
vector, &d);
621 VectorBase::real_type
627 const PetscErrorCode ierr = VecMin (
vector, &p, &d);
634 VectorBase::real_type
640 const PetscErrorCode ierr = VecMax (
vector, &p, &d);
652 const PetscErrorCode ierr = VecAbs (
vector);
665 const PetscErrorCode ierr = VecConjugate (
vector);
689 const PetscErrorCode ierr = VecPointwiseMult (
vector,
vector,v);
701 const PetscErrorCode ierr = VecPointwiseMult (
vector,u,v);
712 PetscScalar *start_ptr;
713 PetscErrorCode ierr = VecGetArray (
vector, &start_ptr);
716 const PetscScalar *ptr = start_ptr,
731 ierr = VecRestoreArray (
vector, &start_ptr);
740 template <
typename T>
741 bool is_non_negative (
const T &t)
748 template <
typename T>
749 bool is_non_negative (
const std::complex<T> &)
753 "whether it is non-negative."))
765 PetscScalar *start_ptr;
766 PetscErrorCode ierr = VecGetArray (vector, &start_ptr);
769 const PetscScalar *ptr = start_ptr,
770 *eptr = start_ptr + local_size();
774 if (! internal::is_non_negative(*ptr))
784 ierr = VecRestoreArray (vector, &start_ptr);
798 const PetscErrorCode ierr = VecScale (vector, a);
812 const PetscScalar factor = 1./a;
815 const PetscErrorCode ierr = VecScale (vector, factor);
827 const PetscErrorCode ierr = VecAXPY (vector, 1, v);
839 const PetscErrorCode ierr = VecAXPY (vector, -1, v);
853 const PetscErrorCode ierr = VecShift (vector, s);
874 const PetscErrorCode ierr = VecAXPY (vector, a, v);
890 const PetscScalar weights[2] = {a,b};
891 Vec addends[2] = {v.
vector, w.vector};
893 const PetscErrorCode ierr = VecMAXPY (vector, 2, weights, addends);
906 const PetscErrorCode ierr = VecAYPX (vector, s, v);
946 const PetscScalar weights[2] = {a,b};
947 Vec addends[2] = {v.
vector,w.vector};
949 const PetscErrorCode ierr = VecMAXPY (vector, 2, weights, addends);
974 const PetscScalar weights[3] = {a,b,c};
977 const PetscErrorCode ierr = VecMAXPY (vector, 3, weights, addends);
987 const PetscErrorCode ierr
988 = VecPointwiseMult (vector, factors, vector);
1007 const PetscErrorCode ierr = VecCopy (v.
vector, vector);
1018 const PetscScalar b,
1031 const PetscErrorCode ierr = VecCopy (v.
vector, vector);
1044 const PetscErrorCode ierr = VecPointwiseDivide (vector, a, b);
1056 PetscErrorCode ierr = PetscViewerSetFormat (PETSC_VIEWER_STDOUT_WORLD,
1061 ierr = VecView (vector, PETSC_VIEWER_STDOUT_WORLD);
1069 const unsigned int precision,
1070 const bool scientific,
1071 const bool across)
const 1078 PetscErrorCode ierr = VecGetArray (vector, &val);
1083 const std::ios::fmtflags old_flags = out.flags();
1084 const unsigned int old_precision = out.precision (precision);
1086 out.precision (precision);
1088 out.setf (std::ios::scientific, std::ios::floatfield);
1090 out.setf (std::ios::fixed, std::ios::floatfield);
1093 for (size_type i=0; i<local_size(); ++i)
1094 out << val[i] <<
' ';
1096 for (size_type i=0; i<local_size(); ++i)
1097 out << val[i] << std::endl;
1101 out.flags (old_flags);
1102 out.precision(old_precision);
1106 ierr = VecRestoreArray (vector, &val);
1117 const PetscErrorCode ierr = VecSwap (vector, v.
vector);
1123 VectorBase::operator
const Vec &()
const 1132 std::size_t mem =
sizeof(Vec)+
sizeof(last_action)
1139 mem += local_size()*
sizeof(PetscScalar);
1143 mem += ghost_indices.n_elements()*(
sizeof(PetscScalar)+
sizeof(
int));
1153 const size_type *indices,
1154 const PetscScalar *values,
1155 const bool add_values)
1160 Assert ((last_action == action)
1163 internal::VectorReference::ExcWrongMode (action,
1172 if (n_elements != 0)
1174 #ifdef PETSC_USE_64BIT_INDICES 1175 std::vector<PetscInt> petsc_ind (n_elements);
1176 for (size_type i=0; i<n_elements; ++i)
1177 petsc_ind[i] = indices[i];
1178 const PetscInt *petsc_indices = &petsc_ind[0];
1180 const int *petsc_indices = (
const int *)indices;
1183 const InsertMode mode = (add_values ? ADD_VALUES : INSERT_VALUES);
1184 const PetscErrorCode ierr = VecSetValues (vector, n_elements, petsc_indices,
1191 last_action = action;
1196 DEAL_II_NAMESPACE_CLOSE
1198 #endif // DEAL_II_WITH_PETSC types::global_dof_index size_type
real_type normalize() const 1
#define AssertNothrow(cond, exc)
VectorBase & operator-=(const VectorBase &V)
void ratio(const VectorBase &a, const VectorBase &b) 1
static ::ExceptionBase & ExcIO()
void sadd(const PetscScalar s, const VectorBase &V)
real_type l2_norm() const
VectorBase & operator/=(const PetscScalar factor)
#define AssertThrow(cond, exc)
real_type norm_sqr() const
VectorBase & operator=(const PetscScalar s)
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
void scale(const VectorBase &scaling_factors)
real_type linfty_norm() const
std::pair< size_type, size_type > local_range() const
real_type lp_norm(const real_type p) const
static ::ExceptionBase & ExcMessage(std::string arg1)
size_type local_size() const
void compress(const VectorOperation::values operation)
void do_set_add_operation(const size_type n_elements, const size_type *indices, const PetscScalar *values, const bool add_values)
VectorOperation::values last_action
#define Assert(cond, exc)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
VectorBase & conjugate() 1
static bool is_running_single_threaded()
PetscScalar add_and_dot(const PetscScalar a, const VectorBase &V, const VectorBase &W)
std::size_t memory_consumption() const
void equ(const PetscScalar a, const VectorBase &V)
VectorBase & operator*=(const PetscScalar factor)
real_type l1_norm() const
static ::ExceptionBase & ExcGhostsPresent()
std_cxx11::enable_if< std_cxx11::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
VectorBase & operator+=(const VectorBase &V)
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
#define AssertThrowMPI(error_code)
PetscScalar operator*(const VectorBase &vec) const
bool has_ghost_elements() const
PetscScalar mean_value() const
bool operator==(const VectorBase &v) const
void add(const std::vector< size_type > &indices, const std::vector< PetscScalar > &values)
void write_ascii(const PetscViewerFormat format=PETSC_VIEWER_DEFAULT)
virtual const MPI_Comm & get_mpi_communicator() const
#define AssertIsFinite(number)
static ::ExceptionBase & ExcInternalError()
bool operator!=(const VectorBase &v) const
void set(const std::vector< size_type > &indices, const std::vector< PetscScalar > &values)