]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Overhaul the PETSc MatrixFree vmult.
authorDavid Wells <wellsd2@rpi.edu>
Sat, 29 Apr 2017 02:12:42 +0000 (22:12 -0400)
committerDavid Wells <wellsd2@rpi.edu>
Thu, 4 May 2017 01:45:40 +0000 (21:45 -0400)
include/deal.II/lac/petsc_matrix_free.h
source/lac/petsc_matrix_free.cc
source/lac/petsc_vector_base.cc

index 01442f4d9b0d4941b19c3bbdb73461d7c4a4d3b0..fee301373f134776c6050ffa2813a66173fe655d 100644 (file)
@@ -48,14 +48,12 @@ namespace PETScWrappers
    * 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
index fba3879f563c55b0283c586a0680b5a483772d89..da9f8c5c7c39b4fc2d9a1c8988ab6f588c146514 100644 (file)
@@ -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<const Vec &>(*y), dst);
-    AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-    delete (x);
-    delete (y);
+    // This is implemented by derived classes
+    vmult (y, x);
   }
 
 
index 83c5fd273230b8e4b4f5f3babd47f0af1a78cda8..ceddd323d8cb1b8a8cb0b343b3c808ec6bac07e4 100644 (file)
@@ -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<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;
     }
   }
 

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.