From 88b34fd76beaa45fc6192018148c5c1d930e781f Mon Sep 17 00:00:00 2001 From: David Wells Date: Thu, 20 Apr 2017 07:11:08 -0400 Subject: [PATCH] Make the PETSc to deal.II vector copy generic. We may as well allow copying instances of the base class: the operations are the same. --- include/deal.II/lac/petsc_parallel_vector.h | 3 + include/deal.II/lac/vector.h | 51 +++------ include/deal.II/lac/vector.templates.h | 121 ++++++++------------ 3 files changed, 70 insertions(+), 105 deletions(-) diff --git a/include/deal.II/lac/petsc_parallel_vector.h b/include/deal.II/lac/petsc_parallel_vector.h index 804bb079ca..2f07585ff3 100644 --- a/include/deal.II/lac/petsc_parallel_vector.h +++ b/include/deal.II/lac/petsc_parallel_vector.h @@ -36,6 +36,9 @@ DEAL_II_NAMESPACE_OPEN */ namespace PETScWrappers { + // forward declaration for copy constructor from sequential vector + class Vector; + /** * Namespace for PETSc classes that work in parallel over MPI, such as * distributed vectors and matrices. diff --git a/include/deal.II/lac/vector.h b/include/deal.II/lac/vector.h index 99a89bc601..b34b9d5cef 100644 --- a/include/deal.II/lac/vector.h +++ b/include/deal.II/lac/vector.h @@ -45,11 +45,7 @@ DEAL_II_NAMESPACE_OPEN #ifdef DEAL_II_WITH_PETSC namespace PETScWrappers { - class Vector; - namespace MPI - { - class Vector; - } + class VectorBase; } #endif @@ -208,23 +204,17 @@ public: #ifdef DEAL_II_WITH_PETSC /** - * Another copy constructor: copy the values from a sequential PETSc wrapper - * vector class. This copy constructor is only available if PETSc was - * detected during configuration time. - */ - explicit Vector (const PETScWrappers::Vector &v); - - /** - * Another copy constructor: copy the values from a parallel PETSc wrapper - * vector class. This copy constructor is only available if PETSc was - * detected during configuration time. + * Another copy constructor: copy the values from a PETSc vector class. This + * copy constructor is only available if PETSc was detected during + * configuration time. * * Note that due to the communication model used in MPI, this operation can - * only succeed if all processes do it at the same time. I.e., it is not - * possible for only one process to obtain a copy of a parallel vector while - * the other jobs do something else. + * only succeed if all processes do it at the same time when v + * is a distributed vector: It is not possible for only one process to + * obtain a copy of a parallel vector while the other jobs do something + * else. */ - explicit Vector (const PETScWrappers::MPI::Vector &v); + explicit Vector (const PETScWrappers::VectorBase &v); #endif #ifdef DEAL_II_WITH_TRILINOS @@ -382,25 +372,18 @@ public: #ifdef DEAL_II_WITH_PETSC /** - * Another copy operator: copy the values from a sequential PETSc wrapper - * vector class. This operator is only available if PETSc was detected - * during configuration time. - */ - Vector & - operator= (const PETScWrappers::Vector &v); - - /** - * Another copy operator: copy the values from a parallel PETSc wrapper - * vector class. This operator is only available if PETSc was detected - * during configuration time. + * Another copy operator: copy the values from a PETSc wrapper vector + * class. This operator is only available if PETSc was detected during + * configuration time. * * Note that due to the communication model used in MPI, this operation can - * only succeed if all processes do it at the same time. I.e., it is not - * possible for only one process to obtain a copy of a parallel vector while - * the other jobs do something else. + * only succeed if all processes do it at the same time when v + * is a distributed vector: It is not possible for only one process to + * obtain a copy of a parallel vector while the other jobs do something + * else. */ Vector & - operator= (const PETScWrappers::MPI::Vector &v); + operator= (const PETScWrappers::VectorBase &v); #endif diff --git a/include/deal.II/lac/vector.templates.h b/include/deal.II/lac/vector.templates.h index 38554edc26..df73ec9ed9 100644 --- a/include/deal.II/lac/vector.templates.h +++ b/include/deal.II/lac/vector.templates.h @@ -25,8 +25,7 @@ #include #ifdef DEAL_II_WITH_PETSC -# include -# include +# include #endif #ifdef DEAL_II_WITH_TRILINOS @@ -34,11 +33,12 @@ #endif +#include #include #include -#include -#include #include +#include +#include DEAL_II_NAMESPACE_OPEN @@ -95,38 +95,56 @@ Vector::Vector (const Vector &v) #ifdef DEAL_II_WITH_PETSC - -template -Vector::Vector (const PETScWrappers::Vector &v) - : - Subscriptor(), - vec_size(v.size()), - max_vec_size(v.size()), - val(nullptr) +namespace internal { - if (vec_size != 0) + template + void + copy_petsc_vector(const PETScWrappers::VectorBase &v, + ::dealii::Vector &out) + { + // Create a sequential PETSc vector and then copy over the entries into + // the deal.II vector. + // + // Wrap the sequential vector with a custom deleter. + std::unique_ptr> sequential_vector + (new Vec(), [](Vec *ptr) { - allocate(); - - // get a representation of the vector - // and copy it - PetscScalar *start_ptr; - PetscErrorCode ierr = VecGetArray (static_cast(v), &start_ptr); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - internal::VectorOperations::copy (start_ptr, start_ptr+vec_size, begin()); - - // restore the representation of the - // vector - ierr = VecRestoreArray (static_cast(v), &start_ptr); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - } +#if DEAL_II_PETSC_VERSION_LT(3,2,0) + const PetscErrorCode ierr = VecDestroy (*ptr); +#else + const PetscErrorCode ierr = VecDestroy (ptr); +#endif + (void)ierr; + AssertNothrow (ierr == 0, ExcPETScError(ierr)); + }); + + PetscErrorCode ierr = VecCreateSeq (v.get_mpi_communicator(), + v.size(), + sequential_vector.get()); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + ierr = VecCopy (v, *sequential_vector); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + PetscScalar *start_ptr; + ierr = VecGetArray(*sequential_vector, &start_ptr); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + const auto v_size = v.size(); + if (out.size() != v_size) + out.reinit (v_size, true); + + internal::VectorOperations::copy (start_ptr, + start_ptr + out.size(), + out.begin()); + ierr = VecRestoreArray (*sequential_vector, &start_ptr); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + } } template -Vector::Vector (const PETScWrappers::MPI::Vector &v) +Vector::Vector (const PETScWrappers::VectorBase &v) : Subscriptor(), vec_size(0), @@ -135,14 +153,9 @@ Vector::Vector (const PETScWrappers::MPI::Vector &v) { if (v.size() != 0) { - // do this in a two-stage process: - // first convert to a sequential petsc - // vector, then copy that - PETScWrappers::Vector seq (v); - *this = seq; + internal::copy_petsc_vector(v, *this); } } - #endif @@ -862,47 +875,13 @@ Vector::operator= (const BlockVector &v) #ifdef DEAL_II_WITH_PETSC - template Vector & -Vector::operator= (const PETScWrappers::Vector &v) +Vector::operator= (const PETScWrappers::VectorBase &v) { - if (v.size() != vec_size) - reinit (v.size(), true); - if (vec_size != 0) - { - // get a representation of the vector - // and copy it - PetscScalar *start_ptr; - PetscErrorCode ierr = VecGetArray (static_cast(v), &start_ptr); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - internal::VectorOperations::copy (start_ptr, start_ptr+vec_size, begin()); - - // restore the representation of the - // vector - ierr = VecRestoreArray (static_cast(v), &start_ptr); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - } - + internal::copy_petsc_vector(v, *this); return *this; } - - - -template -Vector & -Vector::operator= (const PETScWrappers::MPI::Vector &v) -{ - // do this in a two-stage process: - // first convert to a sequential petsc - // vector, then copy that - PETScWrappers::Vector seq (v); - *this = seq; - - return *this; -} - #endif -- 2.39.5