From 185dd0dc85a0afa402806a053256b1757786e813 Mon Sep 17 00:00:00 2001 From: wolf Date: Sun, 4 Apr 2004 23:31:29 +0000 Subject: [PATCH] Have a few more conversion possibilities. Fix a bug where we are only allowed to access local elements, but happily did not. git-svn-id: https://svn.dealii.org/trunk@8963 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/petsc_vector.h | 2 +- deal.II/lac/include/lac/petsc_vector_base.h | 42 +++--------- deal.II/lac/include/lac/vector.h | 75 ++++++++++++++++++--- deal.II/lac/include/lac/vector.templates.h | 74 +++++++++++++++++++- deal.II/lac/source/petsc_vector_base.cc | 70 +++++++++++++++++++ 5 files changed, 221 insertions(+), 42 deletions(-) diff --git a/deal.II/lac/include/lac/petsc_vector.h b/deal.II/lac/include/lac/petsc_vector.h index 6ccee550b5..6598b16012 100644 --- a/deal.II/lac/include/lac/petsc_vector.h +++ b/deal.II/lac/include/lac/petsc_vector.h @@ -223,7 +223,7 @@ namespace PETScWrappers // then do the gather operation ierr = VecConvertMPIToSeqAll (static_cast(v), - &vector); + &vector); AssertThrow (ierr == 0, ExcPETScError(ierr)); return *this; diff --git a/deal.II/lac/include/lac/petsc_vector_base.h b/deal.II/lac/include/lac/petsc_vector_base.h index 796d1d7171..7d26ec2715 100644 --- a/deal.II/lac/include/lac/petsc_vector_base.h +++ b/deal.II/lac/include/lac/petsc_vector_base.h @@ -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(vector), &ptr); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - const PetscScalar value = *(ptr+index); - - ierr = VecRestoreArray (static_cast(vector), &ptr); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - return value; - } } + inline diff --git a/deal.II/lac/include/lac/vector.h b/deal.II/lac/include/lac/vector.h index 553216e83e..0c6ca30127 100644 --- a/deal.II/lac/include/lac/vector.h +++ b/deal.II/lac/include/lac/vector.h @@ -17,11 +17,21 @@ #include #include +#include + + #ifdef DEAL_II_USE_PETSC -# include +namespace PETScWrappers +{ + class Vector; + namespace MPI + { + class Vector; + } +} #endif -#include + /*! @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 & operator= (const Vector &c); - + /** * Copy the given vector. Resize the * present vector if necessary. @@ -271,6 +299,37 @@ class Vector template Vector & operator= (const Vector &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 & + 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 & + operator = (const PETScWrappers::MPI::Vector &v); +#endif + /** * Test for equality. This function * assumes that the present vector and diff --git a/deal.II/lac/include/lac/vector.templates.h b/deal.II/lac/include/lac/vector.templates.h index f6d695d567..4912173dbf 100644 --- a/deal.II/lac/include/lac/vector.templates.h +++ b/deal.II/lac/include/lac/vector.templates.h @@ -15,6 +15,12 @@ #include + +#ifdef DEAL_II_USE_PETSC +# include +# include +#endif + #include #include #include @@ -75,7 +81,7 @@ Vector::Vector (const Vector& v) #ifdef DEAL_II_USE_PETSC template -Vector::Vector (const PETScWrappers::VectorBase &v) +Vector::Vector (const PETScWrappers::Vector &v) : dim(v.size()), maxdim(v.size()), @@ -97,7 +103,26 @@ Vector::Vector (const PETScWrappers::VectorBase &v) // restore the representation of the // vector ierr = VecRestoreArray (static_cast(v), &start_ptr); - AssertThrow (ierr == 0, PETScWrappers::VectorBase::ExcPETScError(ierr)); + AssertThrow (ierr == 0, PETScWrappers::VectorBase::ExcPETScError(ierr)); + } +} + + + +template +Vector::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::operator = (const Vector& v) +#ifdef DEAL_II_USE_PETSC + +template +Vector & +Vector::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(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(v), &start_ptr); + AssertThrow (ierr == 0, PETScWrappers::VectorBase::ExcPETScError(ierr)); + } + + 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 + + template template bool diff --git a/deal.II/lac/source/petsc_vector_base.cc b/deal.II/lac/source/petsc_vector_base.cc index f251402640..c35f77aec6 100644 --- a/deal.II/lac/source/petsc_vector_base.cc +++ b/deal.II/lac/source/petsc_vector_base.cc @@ -13,6 +13,8 @@ #include +#include +#include #include @@ -21,6 +23,74 @@ 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(&vector) != 0) + { + PetscScalar *ptr; + int ierr + = VecGetArray (static_cast(vector), &ptr); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + const PetscScalar value = *(ptr+index); + + ierr = VecRestoreArray (static_cast(vector), &ptr); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + return value; + } + else if (dynamic_cast(&vector) != 0) + { + // first verify that the requested + // element is actually locally + // available + int ierr; + int begin, end; + ierr = VecGetOwnershipRange (static_cast(vector), + &begin, &end); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + AssertThrow ((index >= static_cast(begin)) && + (index < static_cast(end)), + ExcAccessToNonlocalElement (index, begin, end-1)); + + // then access it + PetscScalar *ptr; + ierr = VecGetArray (static_cast(vector), &ptr); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + const PetscScalar value = *(ptr+index-begin); + + ierr = VecRestoreArray (static_cast(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) -- 2.39.5