]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Patch by Fahad Alrasched: Add extract_subvector_to() function to all vectors. Use...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 2 Sep 2013 22:03:06 +0000 (22:03 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 2 Sep 2013 22:03:06 +0000 (22:03 +0000)
git-svn-id: https://svn.dealii.org/trunk@30555 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/dofs/dof_accessor.templates.h
deal.II/include/deal.II/lac/block_vector_base.h
deal.II/include/deal.II/lac/parallel_vector.h
deal.II/include/deal.II/lac/petsc_vector_base.h
deal.II/include/deal.II/lac/trilinos_vector_base.h
deal.II/include/deal.II/lac/vector.h

index c99ce95f9d627e1755679c9143dfe9903b65e3e8..c918de3ea1df515a52811fc4422e45b4d521d313 100644 (file)
@@ -68,6 +68,13 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li>
+  New: All vector classes now have functions <code>extract_subvector_to()</code>
+  that allow extracting not just a single value but a whole set.
+  <br>
+  (Fahad Alrasched, 2013/09/02)
+  </li>
+
   <li>
   Fixed: <code>common/Make.global_options</code> now exports enable-threads
   correctly, furthermore, <code>lib-suffix</code>, <code>shared-lib-suffix</code>
index 3b17b0fa522dba789041d8ce72011d1d6567ba38..4cef5607b27935566c4a72a342b7a357882af14c 100644 (file)
@@ -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<types::global_dof_index> local_dof_indices (dofs_per_cell);
         get_dof_indices (accessor, local_dof_indices);
 
-        for (unsigned int i=0; i<dofs_per_cell; ++i, ++local_values_begin)
-          *local_values_begin = values(local_dof_indices[i]);
+        types::global_dof_index* local_indices_begin = &(local_dof_indices[0]);
+        values.extract_subvector_to (local_indices_begin,
+                                    local_indices_begin + dofs_per_cell,
+                                    local_values_begin);
       }
 
 
index c66dfd47306e58380ec7f5db523bc4bb20f15f3c..899112340815188fc595d0007dc8c4210b792de8 100644 (file)
@@ -927,6 +927,29 @@ public:
    */
   reference 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 <typename OtherNumber>
+  void extract_subvector_to (const std::vector<size_type> &indices,
+                             std::vector<OtherNumber> &values) const;
+
+  /**
+   * Just as the above, but with pointers.
+   * Useful in minimizing copying of data around.
+   */
+  template <typename ForwardIterator, typename OutputIterator>
+  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<VectorType>::operator[] (const size_type i)
   return operator()(i);
 }
 
+
+
+template <typename VectorType>
+template <typename OtherNumber>
+inline
+void BlockVectorBase<VectorType>::extract_subvector_to (const std::vector<size_type> &indices,
+                                           std::vector<OtherNumber> &values) const
+{
+  for (size_type i = 0; i < indices.size(); ++i)
+    values[i] = operator()(indices[i]);
+}
+
+
+
+template <typename VectorType>
+template <typename ForwardIterator, typename OutputIterator>
+inline
+void BlockVectorBase<VectorType>::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
index 613fb61ab4c2048aa3f19c9723792ccfe1bc1db9..6d966154a07b9740f801f1f8d6b7e4efcbd3388d 100644 (file)
@@ -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 <typename OtherNumber>
+      void extract_subvector_to (const std::vector<size_type> &indices,
+                                 std::vector<OtherNumber> &values) const;
+
+      /**
+       * Just as the above, but with pointers.
+       * Useful in minimizing copying of data around.
+       */
+      template <typename ForwardIterator, typename OutputIterator>
+      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 <typename Number>
+    template <typename OtherNumber>
+    inline
+    void Vector<Number>::extract_subvector_to (const std::vector<size_type> &indices,
+                                               std::vector<OtherNumber> &values) const
+    {
+      for (size_type i = 0; i < indices.size(); ++i)
+        values[i] = operator()(indices[i]);
+    }
+
+
+
+    template <typename Number>
+    template <typename ForwardIterator, typename OutputIterator>
+    inline
+    void Vector<Number>::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 <typename Number>
     inline
     Number
index 8d0f7e056f166af8f17470d952eb69b948b32e2a..709e1e61859ab34fd4ccce6ded37945cbc72cf53 100644 (file)
@@ -473,6 +473,28 @@ namespace PETScWrappers
     void set (const std::vector<size_type>   &indices,
               const std::vector<PetscScalar>  &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<size_type> &indices,
+                               std::vector<PetscScalar> &values) const;
+
+    /**
+     * Just as the above, but with pointers.
+     * Useful in minimizing copying of data around.
+     */
+    template <typename ForwardIterator, typename OutputIterator>
+    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<size_type> &indices,
+                                         std::vector<PetscScalar> &values) const
+  {
+    extract_subvector_to(&(indices[0]), &(indices[0]) + indices.size(),  &(values[0]));
+  }
+
+  template <typename ForwardIterator, typename OutputIterator>
+  inline
+  void VectorBase::extract_subvector_to (const ForwardIterator    indices_begin,
+                                         const ForwardIterator    indices_end,
+                                         OutputIterator           values_begin) const
+  {
+    const PetscInt n_idx = static_cast<PetscInt>(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<unsigned int>(begin)
+            && index<static_cast<unsigned int>(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<unsigned int>(begin)
+            && index<static_cast<unsigned int>(end), ExcInternalError());
+
+        *(values_begin+i) = *(ptr+index-begin);
+      }
+
+      ierr = VecRestoreArray(vector, &ptr);
+      AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    }
+  }
+
 #endif // DOXYGEN
 }
 
index 674b6ca0d4f208d816bab78d8c7e3f343082ed25..a721cdce8c95ebd6135697994e81a8b4d8ff4da6 100644 (file)
@@ -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<size_type> &indices,
+                               std::vector<TrilinosScalar> &values) const;
+
+    /**
+     * Just as the above, but with pointers.
+     * Useful in minimizing copying of data around.
+     */
+    template <typename ForwardIterator, typename OutputIterator>
+    void extract_subvector_to (ForwardIterator          indices_begin,
+                               const ForwardIterator    indices_end,
+                               OutputIterator           values_begin) const;
+
     /**
      * Return the value of the vector
      * entry <i>i</i>. Note that this
@@ -1316,6 +1338,30 @@ namespace TrilinosWrappers
 
 
 
+  inline
+  void VectorBase::extract_subvector_to (const std::vector<size_type> &indices,
+                                         std::vector<TrilinosScalar>  &values) const
+  {
+    for (size_type i = 0; i < indices.size(); ++i)
+      values[i] = operator()(indices[i]);
+  }
+
+
+
+  template <typename ForwardIterator, typename OutputIterator>
+  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()
index 6866873d44641e46d87147218df312a7a3a71814..fb67e97c1cd0f7c65d10d23a9b2cc5ceed26b926 100644 (file)
@@ -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 <typename OtherNumber>
+  void extract_subvector_to (const std::vector<size_type> &indices,
+                             std::vector<OtherNumber> &values) const;
+
+  /**
+   * Just as the above, but with pointers.
+   * Useful in minimizing copying of data around.
+   */
+  template <typename ForwardIterator, typename OutputIterator>
+  void extract_subvector_to (ForwardIterator       indices_begin,
+                             const ForwardIterator indices_end,
+                             OutputIterator        values_begin) const;
   //@}
 
 
@@ -1324,6 +1347,33 @@ Number &Vector<Number>::operator[] (const size_type i)
 
 
 
+template <typename Number>
+template <typename OtherNumber>
+inline
+void Vector<Number>::extract_subvector_to (const std::vector<size_type> &indices,
+                                           std::vector<OtherNumber> &values) const
+{
+  for (size_type i = 0; i < indices.size(); ++i)
+    values[i] = operator()(indices[i]);
+}
+
+
+
+template <typename Number>
+template <typename ForwardIterator, typename OutputIterator>
+inline
+void Vector<Number>::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 <typename Number>
 inline
 Vector<Number> &

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.