* multiplication <tt>vmult(VectorBase &dst, const VectorBase &src)</tt>
* which is pure virtual and must be reimplemented in derived classes.
* Besides the usual interface, this class has a matrix-vector
- * multiplication <tt>vmult(Vec &dst, const Vec &src)</tt> taking PETSc
- * Vec objects, which will be called by <tt>matrix_free_mult(Mat A, Vec src,
- * Vec dst)</tt> 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 <tt>Vec*</tt> vectors into a
- * deal.II vector, calls the usual vmult function with the usual interface
- * and converts the result back to PETSc <tt>Vec*</tt>. This could be made
- * much more efficient in derived classes without allocating new memory.
+ * multiplication <tt>vmult(Vec &dst, const Vec &src)</tt> taking PETSc Vec
+ * objects, which will be called by <tt>matrix_free_mult(Mat A, Vec src, Vec
+ * dst)</tt> 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
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<const Vec &>(*y), dst);
- AssertThrow (ierr == 0, ExcPETScError(ierr));
-
- delete (x);
- delete (y);
+ // This is implemented by derived classes
+ vmult (y, x);
}
AssertThrow (ierr == 0, ExcPETScError(ierr));
return value;
}
- // else see if we are dealing
- // with a parallel vector
- else if (dynamic_cast<const PETScWrappers::VectorBase *>(&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<size_type>(begin)
- && index<static_cast<size_type>(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<size_type>(begin)
+ && index<static_cast<size_type>(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<size_type>(begin)) &&
- (index < static_cast<size_type>(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<size_type>(begin)) &&
+ (index < static_cast<size_type>(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;
}
}
projects / dealii.git / commitdiff