]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Have a few more conversion possibilities. Fix a bug where we are only allowed to...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 4 Apr 2004 23:31:29 +0000 (23:31 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 4 Apr 2004 23:31:29 +0000 (23:31 +0000)
git-svn-id: https://svn.dealii.org/trunk@8963 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/petsc_vector.h
deal.II/lac/include/lac/petsc_vector_base.h
deal.II/lac/include/lac/vector.h
deal.II/lac/include/lac/vector.templates.h
deal.II/lac/source/petsc_vector_base.cc

index 6ccee550b5808a4704a7ee97be08746b4ee835da..6598b160129e2c794c1c50d064a2e61936e414a6 100644 (file)
@@ -223,7 +223,7 @@ namespace PETScWrappers
 
                                      // then do the gather operation
     ierr = VecConvertMPIToSeqAll (static_cast<const Vec &>(v),
-                                            &vector);
+                                  &vector);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
     
     return *this;
index 796d1d71714ccad9f77108aa8103b85b7c20b79b..7d26ec27155c73ee293ce9649aacd5790e605aa9 100644 (file)
@@ -119,6 +119,16 @@ namespace PETScWrappers
                         int,
                         << "An error with error number " << arg1
                         << " occured while calling a PETSc function");
+                                         /**
+                                          * Exception
+                                          */
+        DeclException3 (ExcAccessToNonlocalElement,
+                        int, int, int,
+                        << "You tried to access element " << arg1
+                        << " of a distributed vector, but only elements "
+                        << arg2 << " through " << arg3
+                        << " are stored locally and can be accessed.");
+        
       private:
                                          /**
                                           * Point to the vector we are
@@ -762,38 +772,8 @@ namespace PETScWrappers
       
       return *this;
     }
-
-
-
-    inline
-    VectorReference::operator PetscScalar () const
-    {
-                                       // this is clumsy: there is no simple
-                                       // way in PETSc read an element from a
-                                       // vector, i.e. there is no function
-                                       // VecGetValue or so. The only way is
-                                       // to obtain a pointer to a contiguous
-                                       // representation of the vector and
-                                       // read from it. Subsequently, the
-                                       // vector representation has to be
-                                       // restored. If the vector has some
-                                       // kind of non-standard format, such as
-                                       // for parallel vectors, then this is a
-                                       // costly operation, just for a single
-                                       // read access..
-      PetscScalar *ptr;
-      int ierr
-        = VecGetArray (static_cast<const Vec &>(vector), &ptr);
-      AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-      const PetscScalar value = *(ptr+index);
-
-      ierr = VecRestoreArray (static_cast<const Vec &>(vector), &ptr);
-      AssertThrow (ierr == 0, ExcPETScError(ierr));
-      
-      return value;
-    }
   }
+  
 
 
   inline
index 553216e83ea69252584bc6e052a02439676b5b60..0c6ca30127e871b528037d7b65e57383461835d1 100644 (file)
 #include <base/config.h>
 #include <base/exceptions.h>
 
+#include <cstdio>
+
+
 #ifdef DEAL_II_USE_PETSC
-#  include <lac/petsc_vector_base.h>
+namespace PETScWrappers
+{
+  class Vector;
+  namespace MPI
+  {
+    class Vector;
+  }
+}
 #endif
 
-#include <cstdio>
+
 
 /*! @addtogroup Vectors
  *@{
@@ -114,12 +124,30 @@ class Vector
 #ifdef DEAL_II_USE_PETSC
                                      /**
                                       * Another copy constructor: copy the
-                                      * values from a PETSc wrapper vector
-                                      * class. This copy constructor is only
-                                      * available if PETSc was detected during
-                                      * configuration time.
+                                      * values from a sequential PETSc wrapper
+                                      * vector class. This copy constructor is
+                                      * only available if PETSc was detected
+                                      * during configuration time.
+                                      */
+    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.
+                                      *
+                                      * 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.
                                       */
-    Vector (const PETScWrappers::VectorBase &v);
+    Vector (const PETScWrappers::MPI::Vector &v);
 #endif
     
                                     /**
@@ -263,7 +291,7 @@ class Vector
                                      * present vector if necessary.
                                      */
     Vector<Number> & operator= (const Vector<Number> &c);
-
+    
                                     /**
                                      * Copy the given vector. Resize the
                                      * present vector if necessary.
@@ -271,6 +299,37 @@ class Vector
     template <typename Number2>
     Vector<Number> & operator= (const Vector<Number2> &v);
 
+#ifdef DEAL_II_USE_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<Number> &
+    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.
+                                      *
+                                      * 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.
+                                      */
+    Vector<Number> &
+    operator = (const PETScWrappers::MPI::Vector &v);
+#endif
+
                                      /**
                                       * Test for equality. This function
                                       * assumes that the present vector and
index f6d695d567f4bc5e15289befe2f95d039dd39e51..4912173dbf2622b16ba885046aa40a2852d24a47 100644 (file)
 
 
 #include <lac/vector.h>
+
+#ifdef DEAL_II_USE_PETSC
+#  include <lac/petsc_vector.h>
+#  include <lac/petsc_parallel_vector.h>
+#endif
+
 #include <cmath>
 #include <algorithm>
 #include <iostream>
@@ -75,7 +81,7 @@ Vector<Number>::Vector (const Vector<Number>& v)
 #ifdef DEAL_II_USE_PETSC
 
 template <typename Number>
-Vector<Number>::Vector (const PETScWrappers::VectorBase &v)
+Vector<Number>::Vector (const PETScWrappers::Vector &v)
                 :
                dim(v.size()),
                maxdim(v.size()),
@@ -97,7 +103,26 @@ Vector<Number>::Vector (const PETScWrappers::VectorBase &v)
                                        // restore the representation of the
                                        // vector
       ierr = VecRestoreArray (static_cast<const Vec&>(v), &start_ptr);
-      AssertThrow (ierr == 0, PETScWrappers::VectorBase::ExcPETScError(ierr));      
+      AssertThrow (ierr == 0, PETScWrappers::VectorBase::ExcPETScError(ierr));
+    }
+}
+
+
+
+template <typename Number>
+Vector<Number>::Vector (const PETScWrappers::MPI::Vector &v)
+                :
+               dim(v.size()),
+               maxdim(v.size()),
+               val(0)
+{
+  if (dim)
+    {
+                                       // do this in a two-stage process:
+                                       // first convert to a sequential petsc
+                                       // vector, then copy that
+      PETScWrappers::Vector seq (v);
+      *this = seq;
     }
 }
 
@@ -605,6 +630,51 @@ Vector<Number>::operator = (const Vector<Number2>& v)
 
 
 
+#ifdef DEAL_II_USE_PETSC
+
+template <typename Number>
+Vector<Number> &
+Vector<Number>::operator = (const PETScWrappers::Vector &v)
+{
+  if (v.size() != dim)
+    reinit (v.size(), true);
+  if (dim)
+    {
+                                       // get a representation of the vector
+                                       // and copy it
+      PetscScalar *start_ptr;
+      int ierr = VecGetArray (static_cast<const Vec&>(v), &start_ptr);
+      AssertThrow (ierr == 0, PETScWrappers::VectorBase::ExcPETScError(ierr));
+      
+      std::copy (start_ptr, start_ptr+dim, begin());
+
+                                       // restore the representation of the
+                                       // vector
+      ierr = VecRestoreArray (static_cast<const Vec&>(v), &start_ptr);
+      AssertThrow (ierr == 0, PETScWrappers::VectorBase::ExcPETScError(ierr));
+    }
+
+  return *this;
+}
+
+
+
+template <typename Number>
+Vector<Number> &
+Vector<Number>::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
+
+
 template <typename Number>
 template <typename Number2>
 bool
index f251402640ce5514eecb1d6f92ff625c5226bcb6..c35f77aec69fe618278397a3f36c27b011b26ffc 100644 (file)
@@ -13,6 +13,8 @@
 
 
 #include <lac/petsc_vector_base.h>
+#include <lac/petsc_vector.h>
+#include <lac/petsc_parallel_vector.h>
 
 #include <cmath>
 
 
 namespace PETScWrappers
 {
+  namespace internal
+  {
+    VectorReference::operator PetscScalar () const
+    {
+      Assert (index < vector.size(),
+              ExcIndexRange (index, 0, vector.size()));
+              
+                                       // this is clumsy: there is no simple
+                                       // way in PETSc read an element from a
+                                       // vector, i.e. there is no function
+                                       // VecGetValue or so. The only way is
+                                       // to obtain a pointer to a contiguous
+                                       // representation of the vector and
+                                       // read from it. Subsequently, the
+                                       // vector representation has to be
+                                       // restored. In addition, we can only
+                                       // get access to the local part of the
+                                       // vector, so we have to guard against
+                                       // that
+      if (dynamic_cast<const PETScWrappers::Vector *>(&vector) != 0)
+        {
+          PetscScalar *ptr;
+          int ierr
+            = VecGetArray (static_cast<const Vec &>(vector), &ptr);
+          AssertThrow (ierr == 0, ExcPETScError(ierr));
+          
+          const PetscScalar value = *(ptr+index);
+          
+          ierr = VecRestoreArray (static_cast<const Vec &>(vector), &ptr);
+          AssertThrow (ierr == 0, ExcPETScError(ierr));
+          
+          return value;
+        }
+      else if (dynamic_cast<const PETScWrappers::MPI::Vector *>(&vector) != 0)
+        {
+                                           // first verify that the requested
+                                           // element is actually locally
+                                           // available
+          int ierr;
+          int begin, end;
+          ierr = VecGetOwnershipRange (static_cast<const Vec &>(vector),
+                                       &begin, &end);
+          AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+          AssertThrow ((index >= static_cast<unsigned int>(begin)) &&
+                       (index < static_cast<unsigned int>(end)),
+                       ExcAccessToNonlocalElement (index, begin, end-1));
+
+                                           // then access it
+          PetscScalar *ptr;
+          ierr = VecGetArray (static_cast<const Vec &>(vector), &ptr);
+          AssertThrow (ierr == 0, ExcPETScError(ierr));
+          
+          const PetscScalar value = *(ptr+index-begin);
+          
+          ierr = VecRestoreArray (static_cast<const Vec &>(vector), &ptr);
+          AssertThrow (ierr == 0, ExcPETScError(ierr));
+          
+          return value;
+        }
+      else
+                                         // what? what other kind of vector
+                                         // exists there?
+        Assert (false, ExcInternalError());
+      return -1e20;
+    }  
+  }
+  
   VectorBase::VectorBase ()
                   :
                   last_action (LastAction::none)

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.