From ea0d83a2b2fa524e9551083a511ad8ce38f80c14 Mon Sep 17 00:00:00 2001 From: David Wells Date: Fri, 28 Apr 2017 22:12:42 -0400 Subject: [PATCH] Overhaul the PETSc MatrixFree vmult. --- include/deal.II/lac/petsc_matrix_free.h | 14 +-- source/lac/petsc_matrix_free.cc | 59 +-------- source/lac/petsc_vector_base.cc | 156 +++++++++++------------- 3 files changed, 84 insertions(+), 145 deletions(-) diff --git a/include/deal.II/lac/petsc_matrix_free.h b/include/deal.II/lac/petsc_matrix_free.h index 01442f4d9b..fee301373f 100644 --- a/include/deal.II/lac/petsc_matrix_free.h +++ b/include/deal.II/lac/petsc_matrix_free.h @@ -48,14 +48,12 @@ namespace PETScWrappers * multiplication vmult(VectorBase &dst, const VectorBase &src) * which is pure virtual and must be reimplemented in derived classes. * Besides the usual interface, this class has a matrix-vector - * multiplication vmult(Vec &dst, const Vec &src) taking PETSc - * Vec objects, which will be called by matrix_free_mult(Mat A, Vec src, - * Vec dst) registered as matrix-vector multiplication of this PETSc - * matrix object. The default implementation of the vmult function in the - * base class translates the given PETSc Vec* vectors into a - * deal.II vector, calls the usual vmult function with the usual interface - * and converts the result back to PETSc Vec*. This could be made - * much more efficient in derived classes without allocating new memory. + * multiplication vmult(Vec &dst, const Vec &src) taking PETSc Vec + * objects, which will be called by matrix_free_mult(Mat A, Vec src, Vec + * dst) registered as matrix-vector multiplication of this PETSc matrix + * object. The default implementation of the vmult function in the base + * class wraps the given PETSc vectors with the PETScWrappers::VectorBase + * class and then calls the usual vmult function with the usual interface. * * @ingroup PETScWrappers * @ingroup Matrix1 diff --git a/source/lac/petsc_matrix_free.cc b/source/lac/petsc_matrix_free.cc index fba3879f56..da9f8c5c7c 100644 --- a/source/lac/petsc_matrix_free.cc +++ b/source/lac/petsc_matrix_free.cc @@ -173,61 +173,12 @@ namespace PETScWrappers void MatrixFree::vmult (Vec &dst, const Vec &src) const { + // VectorBase permits us to manipulate, but not own, a Vec + PETScWrappers::VectorBase x(src); + PETScWrappers::VectorBase y(dst); -//TODO: Translate the given PETSc Vec* vector into a deal.II -// vector so we can call the vmult function with the usual -// interface; then convert back. This could be much more -// efficient, if the PETScWrappers::*::Vector classes -// had a way to simply generate such a vector object from -// a given PETSc Vec* object without allocating new memory -// and without taking ownership of the Vec* - - VectorBase *x = nullptr; - VectorBase *y = nullptr; - // because we do not know, - // if dst and src are sequential - // or distributed vectors, - // we ask for the vector-type - // and reinit x and y with - // dealii::PETScWrappers::*::Vector: - const char *vec_type; - PetscErrorCode ierr = VecGetType (src, &vec_type); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - PetscInt local_size; - ierr = VecGetLocalSize (src, &local_size); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - if (strcmp(vec_type,"mpi") == 0) - { - PetscInt size; - ierr = VecGetSize (src, &size); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - x = new PETScWrappers::MPI::Vector (this->get_mpi_communicator (), size, local_size); - y = new PETScWrappers::MPI::Vector (this->get_mpi_communicator (), size, local_size); - } - else if (strcmp(vec_type,"seq") == 0) - { - x = new PETScWrappers::Vector (local_size); - y = new PETScWrappers::Vector (local_size); - } - else - AssertThrow (false, ExcMessage("PETScWrappers::MPI::MatrixFree::do_matrix_vector_action: " - "This only works for Petsc Vec Type = VECMPI | VECSEQ")); - - // copy src to x - x->equ(1., PETScWrappers::VectorBase(src)); - // and call vmult(x,y) which must - // be reimplemented in derived classes - vmult (*y, *x); - - // copy the result back to dst - ierr = VecCopy (static_cast(*y), dst); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - delete (x); - delete (y); + // This is implemented by derived classes + vmult (y, x); } diff --git a/source/lac/petsc_vector_base.cc b/source/lac/petsc_vector_base.cc index 83c5fd2732..ceddd323d8 100644 --- a/source/lac/petsc_vector_base.cc +++ b/source/lac/petsc_vector_base.cc @@ -47,103 +47,93 @@ namespace PETScWrappers AssertThrow (ierr == 0, ExcPETScError(ierr)); return value; } - // else see if we are dealing - // with a parallel vector - else if (dynamic_cast(&vector) != nullptr) + + // 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 + if (vector.ghosted) { - // 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 - if (vector.ghosted) - { - PetscInt begin, end; - PetscErrorCode ierr = VecGetOwnershipRange (vector.vector, &begin, - &end); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - Vec locally_stored_elements = PETSC_NULL; - ierr = VecGhostGetLocalForm(vector.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)); - - PetscScalar value; - - if ( index>=static_cast(begin) - && index(end) ) - { - //local entry - value = *(ptr+index-begin); - } - else - { - //ghost entry - const size_type ghostidx - = vector.ghost_indices.index_within_set(index); - - Assert(ghostidx+end-begin<(size_type)lsize, ExcInternalError()); - value = *(ptr+ghostidx+end-begin); - } - - ierr = VecRestoreArray(locally_stored_elements, &ptr); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - ierr = VecGhostRestoreLocalForm(vector.vector, &locally_stored_elements); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - return value; - } + PetscInt begin, end; + PetscErrorCode ierr = VecGetOwnershipRange (vector.vector, &begin, + &end); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + Vec locally_stored_elements = PETSC_NULL; + ierr = VecGhostGetLocalForm(vector.vector, &locally_stored_elements); + AssertThrow (ierr == 0, ExcPETScError(ierr)); - // first verify that the requested - // element is actually locally - // available - PetscInt begin, end; + PetscInt lsize; + ierr = VecGetSize(locally_stored_elements, &lsize); + AssertThrow (ierr == 0, ExcPETScError(ierr)); - PetscErrorCode ierr = VecGetOwnershipRange (vector.vector, &begin, &end); + PetscScalar *ptr; + ierr = VecGetArray(locally_stored_elements, &ptr); AssertThrow (ierr == 0, ExcPETScError(ierr)); + PetscScalar value; + + if ( index>=static_cast(begin) + && index(end) ) + { + //local entry + value = *(ptr+index-begin); + } + else + { + //ghost entry + const size_type ghostidx + = vector.ghost_indices.index_within_set(index); + Assert(ghostidx+end-begin<(size_type)lsize, ExcInternalError()); + value = *(ptr+ghostidx+end-begin); + } - AssertThrow ((index >= static_cast(begin)) && - (index < static_cast(end)), - ExcAccessToNonlocalElement (index, begin, end-1)); + ierr = VecRestoreArray(locally_stored_elements, &ptr); + AssertThrow (ierr == 0, ExcPETScError(ierr)); - // old version which only work with - // VecGetArray()... - PetscInt idx = index; - PetscScalar value; - ierr = VecGetValues(vector.vector, 1, &idx, &value); + ierr = VecGhostRestoreLocalForm(vector.vector, &locally_stored_elements); AssertThrow (ierr == 0, ExcPETScError(ierr)); return value; } - else - // what? what other kind of vector - // exists there? - Assert (false, ExcInternalError()); - return -1e20; + + // first verify that the requested + // element is actually locally + // available + PetscInt begin, end; + + PetscErrorCode ierr = VecGetOwnershipRange (vector.vector, &begin, &end); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + + + AssertThrow ((index >= static_cast(begin)) && + (index < static_cast(end)), + ExcAccessToNonlocalElement (index, begin, end-1)); + + // old version which only work with + // VecGetArray()... + PetscInt idx = index; + PetscScalar value; + ierr = VecGetValues(vector.vector, 1, &idx, &value); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + return value; } } -- 2.39.5