From: Wolfgang Bangerth Date: Mon, 2 Sep 2013 22:03:06 +0000 (+0000) Subject: Patch by Fahad Alrasched: Add extract_subvector_to() function to all vectors. Use... X-Git-Tag: v8.1.0~912 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=40753babd3433a206f7f500d8f89cc1fd3693bad;p=dealii.git Patch by Fahad Alrasched: Add extract_subvector_to() function to all vectors. Use it in DoFAccessor::get_dof_values(). git-svn-id: https://svn.dealii.org/trunk@30555 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index c99ce95f9d..c918de3ea1 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -68,6 +68,13 @@ inconvenience this causes.

Specific improvements

    +
  1. + New: All vector classes now have functions extract_subvector_to() + that allow extracting not just a single value but a whole set. +
    + (Fahad Alrasched, 2013/09/02) +
    2. +
  2. Fixed: common/Make.global_options now exports enable-threads correctly, furthermore, lib-suffix, shared-lib-suffix diff --git a/deal.II/include/deal.II/dofs/dof_accessor.templates.h b/deal.II/include/deal.II/dofs/dof_accessor.templates.h index 3b17b0fa52..4cef5607b2 100644 --- a/deal.II/include/deal.II/dofs/dof_accessor.templates.h +++ b/deal.II/include/deal.II/dofs/dof_accessor.templates.h @@ -2803,8 +2803,10 @@ namespace internal types::global_dof_index *cache = &accessor.dof_handler->levels[accessor.level()] ->cell_dof_indices_cache[accessor.present_index * accessor.get_fe().dofs_per_cell]; - for ( ; local_values_begin != local_values_end; ++local_values_begin, ++cache) - *local_values_begin = values(*cache); + + values.extract_subvector_to (cache, + cache + accessor.get_fe().dofs_per_cell, + local_values_begin); } /** @@ -2834,8 +2836,10 @@ namespace internal std::vector local_dof_indices (dofs_per_cell); get_dof_indices (accessor, local_dof_indices); - for (unsigned int i=0; i + void extract_subvector_to (const std::vector &indices, + std::vector &values) const; + + /** + * Just as the above, but with pointers. + * Useful in minimizing copying of data around. + */ + template + void extract_subvector_to (ForwardIterator indices_begin, + const ForwardIterator indices_end, + OutputIterator values_begin) const; + /** * Copy operator: fill all components of * the vector with the given scalar @@ -2524,6 +2547,33 @@ BlockVectorBase::operator[] (const size_type i) return operator()(i); } + + +template +template +inline +void BlockVectorBase::extract_subvector_to (const std::vector &indices, + std::vector &values) const +{ + for (size_type i = 0; i < indices.size(); ++i) + values[i] = operator()(indices[i]); +} + + + +template +template +inline +void BlockVectorBase::extract_subvector_to (ForwardIterator indices_begin, + const ForwardIterator indices_end, + OutputIterator values_begin) const +{ + while (indices_begin != indices_end) { + *values_begin = operator()(*indices_begin); + indices_begin++; values_begin++; + } +} + #endif // DOXYGEN DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/include/deal.II/lac/parallel_vector.h b/deal.II/include/deal.II/lac/parallel_vector.h index 613fb61ab4..6d966154a0 100644 --- a/deal.II/include/deal.II/lac/parallel_vector.h +++ b/deal.II/include/deal.II/lac/parallel_vector.h @@ -571,6 +571,29 @@ namespace parallel */ Number &operator [] (const size_type global_index); + /** + * A collective get operation: instead + * of getting individual elements of a + * vector, this function allows to get + * a whole set of elements at once. The + * indices of the elements to be read + * are stated in the first argument, + * the corresponding values are returned in the + * second. + */ + template + void extract_subvector_to (const std::vector &indices, + std::vector &values) const; + + /** + * Just as the above, but with pointers. + * Useful in minimizing copying of data around. + */ + template + void extract_subvector_to (ForwardIterator indices_begin, + const ForwardIterator indices_end, + OutputIterator values_begin) const; + /** * Read access to the data field specified by @p local_index. Locally * owned indices can be accessed with indices @@ -1585,6 +1608,33 @@ namespace parallel + template + template + inline + void Vector::extract_subvector_to (const std::vector &indices, + std::vector &values) const + { + for (size_type i = 0; i < indices.size(); ++i) + values[i] = operator()(indices[i]); + } + + + + template + template + inline + void Vector::extract_subvector_to (ForwardIterator indices_begin, + const ForwardIterator indices_end, + OutputIterator values_begin) const + { + while (indices_begin != indices_end) { + *values_begin = operator()(*indices_begin); + indices_begin++; values_begin++; + } + } + + + template inline Number diff --git a/deal.II/include/deal.II/lac/petsc_vector_base.h b/deal.II/include/deal.II/lac/petsc_vector_base.h index 8d0f7e056f..709e1e6185 100644 --- a/deal.II/include/deal.II/lac/petsc_vector_base.h +++ b/deal.II/include/deal.II/lac/petsc_vector_base.h @@ -473,6 +473,28 @@ namespace PETScWrappers void set (const std::vector &indices, const std::vector &values); + /** + * A collective get operation: instead + * of getting individual elements of a + * vector, this function allows to get + * a whole set of elements at once. The + * indices of the elements to be read + * are stated in the first argument, + * the corresponding values are returned in the + * second. + */ + void extract_subvector_to (const std::vector &indices, + std::vector &values) const; + + /** + * Just as the above, but with pointers. + * Useful in minimizing copying of data around. + */ + template + void extract_subvector_to (const ForwardIterator indices_begin, + const ForwardIterator indices_end, + OutputIterator values_begin) const; + /** * A collective add operation: This * function adds a whole set of values @@ -1222,6 +1244,119 @@ namespace PETScWrappers return comm; } + inline + void VectorBase::extract_subvector_to (const std::vector &indices, + std::vector &values) const + { + extract_subvector_to(&(indices[0]), &(indices[0]) + indices.size(), &(values[0])); + } + + template + inline + void VectorBase::extract_subvector_to (const ForwardIterator indices_begin, + const ForwardIterator indices_end, + OutputIterator values_begin) const + { + const PetscInt n_idx = static_cast(indices_end - indices_begin); + if (n_idx == 0) + return; + + // if we are dealing + // with a parallel vector + if (ghosted ) + { + + int ierr; + + // there is the possibility + // that the vector has + // ghost elements. in that + // case, we first need to + // figure out which + // elements we own locally, + // then get a pointer to + // the elements that are + // stored here (both the + // ones we own as well as + // the ghost elements). in + // this array, the locally + // owned elements come + // first followed by the + // ghost elements whose + // position we can get from + // an index set + PetscInt begin, end, i; + ierr = VecGetOwnershipRange (vector, &begin, &end); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + Vec locally_stored_elements = PETSC_NULL; + ierr = VecGhostGetLocalForm(vector, &locally_stored_elements); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + PetscInt lsize; + ierr = VecGetSize(locally_stored_elements, &lsize); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + PetscScalar *ptr; + ierr = VecGetArray(locally_stored_elements, &ptr); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + for (i = 0; i < n_idx; i++) { + const unsigned int index = *(indices_begin+i); + if ( index>=static_cast(begin) + && index(end) ) + { + //local entry + *(values_begin+i) = *(ptr+index-begin); + } + else + { + //ghost entry + const unsigned int ghostidx + = ghost_indices.index_within_set(index); + + Assert(ghostidx+end-begin<(unsigned int)lsize, ExcInternalError()); + *(values_begin+i) = *(ptr+ghostidx+end-begin); + } + } + + ierr = VecRestoreArray(locally_stored_elements, &ptr); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + ierr = VecGhostRestoreLocalForm(vector, &locally_stored_elements); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + } + // if the vector is local or the + // caller, then simply access the + // element we are interested in + else + { + int ierr; + + PetscInt begin, end; + ierr = VecGetOwnershipRange (vector, &begin, &end); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + PetscScalar *ptr; + ierr = VecGetArray(vector, &ptr); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + for (PetscInt i = 0; i < n_idx; i++) { + const unsigned int index = *(indices_begin+i); + + Assert(index>=static_cast(begin) + && index(end), ExcInternalError()); + + *(values_begin+i) = *(ptr+index-begin); + } + + ierr = VecRestoreArray(vector, &ptr); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + } + } + #endif // DOXYGEN } diff --git a/deal.II/include/deal.II/lac/trilinos_vector_base.h b/deal.II/include/deal.II/lac/trilinos_vector_base.h index 674b6ca0d4..a721cdce8c 100644 --- a/deal.II/include/deal.II/lac/trilinos_vector_base.h +++ b/deal.II/include/deal.II/lac/trilinos_vector_base.h @@ -674,6 +674,28 @@ namespace TrilinosWrappers TrilinosScalar operator [] (const size_type index) const; + /** + * A collective get operation: instead + * of getting individual elements of a + * vector, this function allows to get + * a whole set of elements at once. The + * indices of the elements to be read + * are stated in the first argument, + * the corresponding values are returned in the + * second. + */ + void extract_subvector_to (const std::vector &indices, + std::vector &values) const; + + /** + * Just as the above, but with pointers. + * Useful in minimizing copying of data around. + */ + template + void extract_subvector_to (ForwardIterator indices_begin, + const ForwardIterator indices_end, + OutputIterator values_begin) const; + /** * Return the value of the vector * entry i. Note that this @@ -1316,6 +1338,30 @@ namespace TrilinosWrappers + inline + void VectorBase::extract_subvector_to (const std::vector &indices, + std::vector &values) const + { + for (size_type i = 0; i < indices.size(); ++i) + values[i] = operator()(indices[i]); + } + + + + template + inline + void VectorBase::extract_subvector_to (ForwardIterator indices_begin, + const ForwardIterator indices_end, + OutputIterator values_begin) const + { + while (indices_begin != indices_end) { + *values_begin = operator()(*indices_begin); + indices_begin++; values_begin++; + } + } + + + inline VectorBase::iterator VectorBase::begin() diff --git a/deal.II/include/deal.II/lac/vector.h b/deal.II/include/deal.II/lac/vector.h index 6866873d44..fb67e97c1c 100644 --- a/deal.II/include/deal.II/lac/vector.h +++ b/deal.II/include/deal.II/lac/vector.h @@ -717,6 +717,29 @@ public: * Exactly the same as operator(). */ Number &operator[] (const size_type i); + + /** + * A collective get operation: instead + * of getting individual elements of a + * vector, this function allows to get + * a whole set of elements at once. The + * indices of the elements to be read + * are stated in the first argument, + * the corresponding values are returned in the + * second. + */ + template + void extract_subvector_to (const std::vector &indices, + std::vector &values) const; + + /** + * Just as the above, but with pointers. + * Useful in minimizing copying of data around. + */ + template + void extract_subvector_to (ForwardIterator indices_begin, + const ForwardIterator indices_end, + OutputIterator values_begin) const; //@} @@ -1324,6 +1347,33 @@ Number &Vector::operator[] (const size_type i) +template +template +inline +void Vector::extract_subvector_to (const std::vector &indices, + std::vector &values) const +{ + for (size_type i = 0; i < indices.size(); ++i) + values[i] = operator()(indices[i]); +} + + + +template +template +inline +void Vector::extract_subvector_to (ForwardIterator indices_begin, + const ForwardIterator indices_end, + OutputIterator values_begin) const +{ + while (indices_begin != indices_end) { + *values_begin = operator()(*indices_begin); + indices_begin++; values_begin++; + } +} + + + template inline Vector &