// then do the gather operation
ierr = VecConvertMPIToSeqAll (static_cast<const Vec &>(v),
- &vector);
+ &vector);
AssertThrow (ierr == 0, ExcPETScError(ierr));
return *this;
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
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
#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
*@{
#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
/**
* present vector if necessary.
*/
Vector<Number> & operator= (const Vector<Number> &c);
-
+
/**
* Copy the given vector. Resize the
* present vector if necessary.
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
#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>
#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()),
// 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;
}
}
+#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
#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)
projects / dealii-svn.git / commitdiff