Dear Roy,
On Wed, 25 Feb 2009, Tim Kroeger wrote:
On Wed, 25 Feb 2009, Tim Kroeger wrote:
Besides, this can't be the right fix, can it? Why should any parallel
communication be required to copy one consistent ghosted vector to
another? It seems as if VecCopy is simply failing to copy the whole
vec, and I don't like not understanding why.
Good question. Perhaps one should apply VecCopy to the local form
(obtained via VecGhostGetLocalForm())? Also, I wonder whether other
methods, such as VecScale(), might possibly suffer from similar
problems?
I'll ask in the petsc-users list and let you know what they say.
The PETSc people say that the local form has to be used (or
VecGhostUpdate*(), as we used in PetscVector::close()). See also
Jed's reply.
I'll go through the PetscVector code once more and change this in all
methods that are expected to work without parallel communication.
This will most probably also happen next week.
I've done this now, see attachment. I explicitly used the local form
for the ghosted case in the following methods:
PetscVector::zero(void)
PetscVector::add(const T)
PetscVector::add(const T, const NumericVector&)
PetscVector::abs(void)
PetscVector::operator=(const T)
PetscVector::operator=(const PetscVector&)
PetscVector::pointwiseMult(const NumericVector&, const NumericVector&)
I wonder whether we should do this in PetscVector::operator=(const
std::vector<T>&). At least in its "case 1", it should be possible to
do without parallel communication using the local form, but in "case
2", it's of course impossible. I left this method unchanged, since
the version using this->close() should certainly work as well.
All the other functions seem not to require the local form to me,
because they either (a) require closing the vector afterwards, or (b)
perform parallel communication, or (c) do not modify the vector, or
(d) do not directly use PETSc methods.
It would be nice if you would double-check my changes before
committing. One never knows.
Best Regards,
Tim
--
Dr. Tim Kroeger
tim.kroe...@mevis.fraunhofer.de Phone +49-421-218-7710
tim.kroe...@cevis.uni-bremen.de Fax +49-421-218-4236
Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany
Index: include/numerics/petsc_vector.h
===================================================================
--- include/numerics/petsc_vector.h (Revision 3298)
+++ include/numerics/petsc_vector.h (Arbeitskopie)
@@ -821,19 +821,37 @@
PetscScalar z=0.;
+ if(this->type() != GHOSTED)
+ {
#if PETSC_VERSION_LESS_THAN(2,3,0)
-
- // 2.2.x & earlier style
- ierr = VecSet (&z, _vec);
- CHKERRABORT(libMesh::COMM_WORLD,ierr);
-
+ // 2.2.x & earlier style
+ ierr = VecSet (&z, _vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
#else
-
- // 2.3.x & newer
- ierr = VecSet (_vec, z);
- CHKERRABORT(libMesh::COMM_WORLD,ierr);
-
+ // 2.3.x & newer
+ ierr = VecSet (_vec, z);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
#endif
+ }
+ else
+ {
+ /* Vectors that include ghost values require a special
+ handling. */
+ Vec loc_vec;
+ ierr = VecGhostGetLocalForm (_vec,&loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+#if PETSC_VERSION_LESS_THAN(2,3,0)
+ // 2.2.x & earlier style
+ ierr = VecSet (&z, loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+#else
+ // 2.3.x & newer
+ ierr = VecSet (loc_vec, z);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+#endif
+ ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
}
Index: src/numerics/petsc_vector.C
===================================================================
--- src/numerics/petsc_vector.C (Revision 3298)
+++ src/numerics/petsc_vector.C (Arbeitskopie)
@@ -234,23 +234,56 @@
int ierr=0;
PetscScalar* values;
const PetscScalar v = static_cast<PetscScalar>(v_in);
- const int n = static_cast<int>(this->local_size());
- const int fli = static_cast<int>(this->first_local_index());
-
- for (int i=0; i<n; i++)
+
+ if(this->type() != GHOSTED)
{
- ierr = VecGetArray (_vec, &values);
- CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ const int n = static_cast<int>(this->local_size());
+ const int fli = static_cast<int>(this->first_local_index());
- int ig = fli + i;
+ for (int i=0; i<n; i++)
+ {
+ ierr = VecGetArray (_vec, &values);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ int ig = fli + i;
+
+ PetscScalar value = (values[i] + v);
+
+ ierr = VecRestoreArray (_vec, &values);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ ierr = VecSetValues (_vec, 1, &ig, &value, INSERT_VALUES);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
+ }
+ else
+ {
+ /* Vectors that include ghost values require a special
+ handling. */
+ Vec loc_vec;
+ ierr = VecGhostGetLocalForm (_vec,&loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ int n=0;
+ ierr = VecGetSize(loc_vec, &n);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
- PetscScalar value = (values[i] + v);
-
- ierr = VecRestoreArray (_vec, &values);
- CHKERRABORT(libMesh::COMM_WORLD,ierr);
-
- ierr = VecSetValues (_vec, 1, &ig, &value, INSERT_VALUES);
- CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ for (int i=0; i<n; i++)
+ {
+ ierr = VecGetArray (loc_vec, &values);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ PetscScalar value = (values[i] + v);
+
+ ierr = VecRestoreArray (loc_vec, &values);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ ierr = VecSetValues (loc_vec, 1, &i, &value, INSERT_VALUES);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
+
+ ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
}
this->_is_closed = false;
@@ -277,20 +310,42 @@
libmesh_assert(this->size() == v->size());
+ if(this->type() != GHOSTED)
+ {
+#if PETSC_VERSION_LESS_THAN(2,3,0)
+ // 2.2.x & earlier style
+ ierr = VecAXPY(&a, v->_vec, _vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+#else
+ // 2.3.x & later style
+ ierr = VecAXPY(_vec, a, v->_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+#endif
+ }
+ else
+ {
+ Vec loc_vec;
+ Vec v_loc_vec;
+ ierr = VecGhostGetLocalForm (_vec,&loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = VecGhostGetLocalForm (v->_vec,&v_loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
#if PETSC_VERSION_LESS_THAN(2,3,0)
-
- // 2.2.x & earlier style
- ierr = VecAXPY(&a, v->_vec, _vec);
- CHKERRABORT(libMesh::COMM_WORLD,ierr);
-
+ // 2.2.x & earlier style
+ ierr = VecAXPY(&a, v_loc_vec, loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
#else
-
- // 2.3.x & later style
- ierr = VecAXPY(_vec, a, v->_vec);
- CHKERRABORT(libMesh::COMM_WORLD,ierr);
-
+ // 2.3.x & later style
+ ierr = VecAXPY(loc_vec, a, v_loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
#endif
+
+ ierr = VecGhostRestoreLocalForm (v->_vec,&v_loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
}
@@ -349,27 +404,61 @@
int ierr = 0;
PetscScalar factor = static_cast<PetscScalar>(factor_in);
+ if(this->type() != GHOSTED)
+ {
#if PETSC_VERSION_LESS_THAN(2,3,0)
+ // 2.2.x & earlier style
+ ierr = VecScale(&factor, _vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+#else
+ // 2.3.x & later style
+ ierr = VecScale(_vec, factor);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+#endif
+ }
+ else
+ {
+ Vec loc_vec;
+ ierr = VecGhostGetLocalForm (_vec,&loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
- // 2.2.x & earlier style
- ierr = VecScale(&factor, _vec);
- CHKERRABORT(libMesh::COMM_WORLD,ierr);
-
+#if PETSC_VERSION_LESS_THAN(2,3,0)
+ // 2.2.x & earlier style
+ ierr = VecScale(&factor, loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
#else
-
- // 2.3.x & later style
- ierr = VecScale(_vec, factor);
- CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ // 2.3.x & later style
+ ierr = VecScale(loc_vec, factor);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+#endif
-#endif
+ ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
}
template <typename T>
void PetscVector<T>::abs()
{
int ierr = 0;
- ierr = VecAbs(_vec);
- CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ if(this->type() != GHOSTED)
+ {
+ ierr = VecAbs(_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
+ else
+ {
+ Vec loc_vec;
+ ierr = VecGhostGetLocalForm (_vec,&loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ ierr = VecAbs(loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
}
template <typename T>
@@ -405,19 +494,37 @@
if (this->size() != 0)
{
+ if(this->type() != GHOSTED)
+ {
#if PETSC_VERSION_LESS_THAN(2,3,0)
-
- // 2.2.x & earlier style
- ierr = VecSet(&s, _vec);
- CHKERRABORT(libMesh::COMM_WORLD,ierr);
-
+ // 2.2.x & earlier style
+ ierr = VecSet(&s, _vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
#else
+ // 2.3.x & later style
+ ierr = VecSet(_vec, s);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+#endif
+ }
+ else
+ {
+ Vec loc_vec;
+ ierr = VecGhostGetLocalForm (_vec,&loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
- // 2.3.x & later style
- ierr = VecSet(_vec, s);
- CHKERRABORT(libMesh::COMM_WORLD,ierr);
-
+#if PETSC_VERSION_LESS_THAN(2,3,0)
+ // 2.2.x & earlier style
+ ierr = VecSet(&s, loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+#else
+ // 2.3.x & later style
+ ierr = VecSet(loc_vec, s);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
#endif
+
+ ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
}
return *this;
@@ -453,14 +560,31 @@
if (v.size() != 0)
{
- int ierr = VecCopy (v._vec, this->_vec);
- CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ int ierr = 0;
+ if(this->type() != GHOSTED)
+ {
+ ierr = VecCopy (v._vec, this->_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
+ else
+ {
+ Vec loc_vec;
+ Vec v_loc_vec;
+ ierr = VecGhostGetLocalForm (_vec,&loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = VecGhostGetLocalForm (v._vec,&v_loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ ierr = VecCopy (v_loc_vec, loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ ierr = VecGhostRestoreLocalForm (v._vec,&v_loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
}
- // Make sure ghost dofs are copied
- if (this->type() == GHOSTED)
- this->close();
-
return *this;
}
@@ -934,12 +1058,37 @@
libmesh_error();
#else
+
+ if(this->type() != GHOSTED)
+ {
+ ierr = VecPointwiseMult(this->vec(),
+ const_cast<PetscVector<T>*>(vec1_petsc)->vec(),
+ const_cast<PetscVector<T>*>(vec2_petsc)->vec());
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
+ else
+ {
+ Vec loc_vec;
+ Vec v1_loc_vec;
+ Vec v2_loc_vec;
+ ierr = VecGhostGetLocalForm (_vec,&loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = VecGhostGetLocalForm
(const_cast<PetscVector<T>*>(vec1_petsc)->vec(),&v1_loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = VecGhostGetLocalForm
(const_cast<PetscVector<T>*>(vec2_petsc)->vec(),&v2_loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
- ierr = VecPointwiseMult(this->vec(),
- const_cast<PetscVector<T>*>(vec1_petsc)->vec(),
- const_cast<PetscVector<T>*>(vec2_petsc)->vec());
- CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = VecPointwiseMult(loc_vec,v1_loc_vec,v2_loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = VecGhostRestoreLocalForm
(const_cast<PetscVector<T>*>(vec1_petsc)->vec(),&v1_loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = VecGhostRestoreLocalForm
(const_cast<PetscVector<T>*>(vec2_petsc)->vec(),&v2_loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
+
#endif
}
------------------------------------------------------------------------------
Open Source Business Conference (OSBC), March 24-25, 2009, San Francisco, CA
-OSBC tackles the biggest issue in open source: Open Sourcing the Enterprise
-Strategies to boost innovation and cut costs with open source participation
-Receive a $600 discount off the registration fee with the source code: SFAD
http://p.sf.net/sfu/XcvMzF8H
_______________________________________________
Libmesh-devel mailing list
Libmesh-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-devel
