diff --git a/MathLib/LinAlg/PETSc/PETScVector.cpp b/MathLib/LinAlg/PETSc/PETScVector.cpp
index a9a37c676ae85989a7d87d7d9c7a320f7dd392b8..ebd841c61a258bfe7114f36160a72da0875037aa 100644
--- a/MathLib/LinAlg/PETSc/PETScVector.cpp
+++ b/MathLib/LinAlg/PETSc/PETScVector.cpp
@@ -28,10 +28,13 @@ namespace MathLib
{
PETScVector::PETScVector(const PetscInt vec_size, const bool is_global_size)
{
- if( is_global_size ) {
+ if( is_global_size )
+ {
VecCreate(PETSC_COMM_WORLD, &_v);
VecSetSizes(_v, PETSC_DECIDE, vec_size);
- } else {
+ }
+ else
+ {
// Fix size partitioning
// the size can be associated to specific memory allocation of a matrix
VecCreateMPI(PETSC_COMM_WORLD, vec_size, PETSC_DECIDE, &_v);
@@ -43,8 +46,8 @@ PETScVector::PETScVector(const PetscInt vec_size, const bool is_global_size)
PETScVector::PETScVector(const PetscInt vec_size,
const std::vector<PetscInt>& ghost_ids,
const bool is_global_size)
- : _size_ghosts{static_cast<PetscInt>(ghost_ids.size())}
- , _has_ghost_id{true}
+ : _size_ghosts {static_cast<PetscInt>(ghost_ids.size())},
+ _has_ghost_id {true}
{
PetscInt nghosts = static_cast<PetscInt>( ghost_ids.size() );
if ( is_global_size )
@@ -61,6 +64,10 @@ PETScVector::PETScVector(const PetscInt vec_size,
}
config();
+
+ for (PetscInt i=0; i<nghosts; i++)
+ _global_ids2local_ids_ghost.insert
+ (std::make_pair(ghost_ids[i], _size_loc + i));
}
PETScVector::PETScVector(const PETScVector &existing_vec, const bool deep_copy)
@@ -75,14 +82,15 @@ PETScVector::PETScVector(const PETScVector &existing_vec, const bool deep_copy)
}
PETScVector::PETScVector(PETScVector &&other)
- : _v{std::move(other._v)}
- , _v_loc{std::move(other._v_loc)}
- , _start_rank{other._start_rank}
- , _end_rank{other._end_rank}
- , _size{other._size}
- , _size_loc{other._size_loc}
- , _size_ghosts{other._size_ghosts}
- , _has_ghost_id{other._has_ghost_id}
+ : _v {std::move(other._v)},
+ _v_loc {std::move(other._v_loc)},
+ _start_rank {other._start_rank},
+ _end_rank {other._end_rank},
+ _size {other._size},
+ _size_loc {other._size_loc},
+ _size_ghosts {other._size_ghosts},
+ _has_ghost_id {other._has_ghost_id},
+ _global_ids2local_ids_ghost {other._global_ids2local_ids_ghost}
{}
void PETScVector::config()
@@ -130,7 +138,7 @@ void PETScVector::gatherLocalVectors( PetscScalar local_array[],
}
-void PETScVector::getGlobalVector(PetscScalar u[]) const
+void PETScVector::getGlobalVector(std::vector<PetscScalar>& u) const
{
#ifdef TEST_MEM_PETSC
@@ -141,7 +149,7 @@ void PETScVector::getGlobalVector(PetscScalar u[]) const
PetscScalar *xp = nullptr;
VecGetArray(_v, &xp);
- gatherLocalVectors(xp, u);
+ gatherLocalVectors(xp, &u[0]);
//This following line may be needed late on
// for a communication load balance:
@@ -158,13 +166,23 @@ void PETScVector::getGlobalVector(PetscScalar u[]) const
void PETScVector::setLocalAccessibleVector() const
{
- // TODO: use getLocalVector
- if (!_global_v)
+ if (_entry_array.size() == 0)
{
- _global_v = std::unique_ptr<PetscScalar[]>{ new PetscScalar[_size] };
+ const PetscInt array_size
+ = _global_ids2local_ids_ghost.size() > 0 ?
+ _size_loc + _size_ghosts: _size;
+ _entry_array.resize(array_size);
}
- getGlobalVector(_global_v.get());
+ if (_global_ids2local_ids_ghost.size() > 0)
+ {
+ double* loc_x = getLocalVector();
+ std::copy_n(loc_x, _size_loc + _size_ghosts,
+ _entry_array.begin());
+ restoreArray(loc_x);
+ }
+ else
+ getGlobalVector(_entry_array);
}
void PETScVector::copyValues(std::vector<double>& u) const
@@ -176,6 +194,22 @@ void PETScVector::copyValues(std::vector<double>& u) const
restoreArray(loc_x);
}
+PetscScalar PETScVector::get(const PetscInt idx) const
+{
+ if (_global_ids2local_ids_ghost.size() > 0)
+ {
+ return _entry_array[getLocalIndex(idx)];
+ }
+ else
+ {
+ // Ghost entries, and its original index is 0.
+ const PetscInt id_p = (idx == -_size) ? 0 : std::abs(idx);
+ return _entry_array[id_p];
+ }
+ return 0; // avoid warning.
+}
+
+
std::vector<double> PETScVector::get(std::vector<IndexType> const& indices) const
{
std::vector<double> local_x(indices.size());
@@ -183,12 +217,22 @@ std::vector<double> PETScVector::get(std::vector<IndexType> const& indices) cons
// use VecGetValues(_v, indices.size(), indices.data(),
// local_x.data());
- for (std::size_t i=0; i<indices.size(); i++)
+ if (_global_ids2local_ids_ghost.size() > 0)
{
- // Ghost entries, and its original index is 0.
- const IndexType id_p = (indices[i] == -_size)
- ? 0 : std::abs(indices[i]);
- local_x[i] = _global_v[id_p];
+ for (std::size_t i=0; i<indices.size(); i++)
+ {
+ local_x[i] = _entry_array[getLocalIndex(indices[i])];
+ }
+ }
+ else
+ {
+ for (std::size_t i=0; i<indices.size(); i++)
+ {
+ // Ghost entries, and its original index is 0.
+ const IndexType id_p = (indices[i] == -_size)
+ ? 0 : std::abs(indices[i]);
+ local_x[i] = _entry_array[id_p];
+ }
}
return local_x;
@@ -205,10 +249,23 @@ PetscScalar* PETScVector::getLocalVector() const
VecGetArray(_v_loc, &loc_array);
}
else
- VecGetArray(_v, &loc_array);
+ VecGetArray(_v, &loc_array);
return loc_array;
}
+PetscInt PETScVector::getLocalIndex(const PetscInt global_index) const
+{
+ if (global_index >= 0) // non-ghost entry.
+ return global_index - _start_rank;
+
+ // A special case for a ghost location with global index equal to
+ // the size of the local vector:
+ PetscInt real_global_index =
+ (-global_index == _size) ? 0 : -global_index;
+
+ return _global_ids2local_ids_ghost.at(real_global_index);
+}
+
void PETScVector::restoreArray(PetscScalar* array) const
{
if (_has_ghost_id)
@@ -250,6 +307,7 @@ void PETScVector::shallowCopy(const PETScVector &v)
_size_loc = v._size_loc;
_size_ghosts = v._size_ghosts;
_has_ghost_id = v._has_ghost_id;
+ _global_ids2local_ids_ghost = v._global_ids2local_ids_ghost;
VecSetOption(_v, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
}
diff --git a/MathLib/LinAlg/PETSc/PETScVector.h b/MathLib/LinAlg/PETSc/PETScVector.h
index fb982912770a996066cb12c057f9663d2100cefa..187f8758f5fb97ff4714e114988315c2dccffdfd 100644
--- a/MathLib/LinAlg/PETSc/PETScVector.h
+++ b/MathLib/LinAlg/PETSc/PETScVector.h
@@ -16,9 +16,11 @@
#pragma once
+#include <map>
#include <memory>
#include <string>
#include <vector>
+
#include <petscvec.h>
namespace MathLib
@@ -162,18 +164,17 @@ class PETScVector
/// Get several entries
std::vector<double> get(std::vector<IndexType> const& indices) const;
- // TODO preliminary
+ // Get the value of an entry by [] operator
double operator[] (PetscInt idx) const
{
- const PetscInt id_p = (idx == -_size) ? 0 : std::abs(idx);
- return _global_v[id_p];
+ return get(idx);
}
/*!
Get global vector
\param u Array to store the global vector. Memory allocation is needed in advance
*/
- void getGlobalVector(PetscScalar u[]) const;
+ void getGlobalVector(std::vector<PetscScalar>& u) const;
/*!
Copy local entries including ghost ones to an array
@@ -183,14 +184,7 @@ class PETScVector
/// Get an entry value. This is an expensive operation,
/// and it only get local value. Use it for only test purpose
- PetscScalar get(const PetscInt idx) const
- {
- //PetscScalar x;
- //VecGetValues(_v, 1, &idx, &x);
- const PetscInt id_p = (idx == -_size) ? 0 : std::abs(idx);
- //VecGetValues(_v, 1, &id_p, &value);
- return _global_v[id_p];
- }
+ PetscScalar get(const PetscInt idx) const;
// TODO eliminate in favour of getRawVector()
/// Get PETsc vector. Use it only for test purpose
@@ -216,7 +210,7 @@ class PETScVector
* This method is dangerous insofar as you can do arbitrary things also
* with a const PETSc vector!
*/
- PETSc_Vec getRawVector() const {return _v; }
+ PETSc_Vec getRawVector() const { return _v; }
/*! View the global vector for test purpose. Do not use it for output a big vector.
@@ -252,7 +246,7 @@ class PETScVector
PETSc_Vec _v = nullptr;
/// Local vector, which is only for the case that _v is created
- /// with ghost entries.
+ /// with ghost entries.
mutable PETSc_Vec _v_loc = nullptr;
/// Starting index in a rank
@@ -271,11 +265,16 @@ class PETScVector
bool _has_ghost_id = false;
/*!
- Scattered global vector to all processors. Note: this member
- and its associated computations can be dropped if
- VecGetValues can get values from different processors.
+ \brief Array containing the entries of the vector.
+ If the vector is created without given ghost IDs, the array
+ contains all entries of the global vector, _v. Otherwise it
+ only contains the entries of the global vector owned by the
+ current rank.
*/
- mutable std::unique_ptr<PetscScalar[]> _global_v;
+ mutable std::vector<PetscScalar> _entry_array;
+
+ /// Map global indices of ghost enrties to local indices
+ mutable std::map<PetscInt, PetscInt> _global_ids2local_ids_ghost;
/*!
\brief Collect local vectors
@@ -290,6 +289,9 @@ class PETScVector
*/
PetscScalar* getLocalVector() const;
+ /// Get local index by a global index
+ PetscInt getLocalIndex(const PetscInt global_index) const;
+
/*!
Restore array after finish access local array
\param array Pointer to the local array fetched by VecGetArray
diff --git a/Tests/MathLib/TestGlobalVectorInterface.cpp b/Tests/MathLib/TestGlobalVectorInterface.cpp
index b0baba08a8f105d2f0b46bcacb3e8e5d1544b4dc..51582e59f1ad9f20857c41479910b892d751978f 100644
--- a/Tests/MathLib/TestGlobalVectorInterface.cpp
+++ b/Tests/MathLib/TestGlobalVectorInterface.cpp
@@ -135,7 +135,7 @@ void checkGlobalVectorInterfacePETSc()
ASSERT_NEAR(0.0, normy-norm2(y), 1.e-10);
- double x0[16];
+ std::vector<double> x0(16);
double z[] =
{
2.0000000000000000e+01,
diff --git a/Tests/MathLib/TestLinearSolver.cpp b/Tests/MathLib/TestLinearSolver.cpp
index 358cea04d444e8a6f576a840da55e1f9fe8cb46c..f1f1e92afa2923d5a37f734e18f69c7173b7e82c 100644
--- a/Tests/MathLib/TestLinearSolver.cpp
+++ b/Tests/MathLib/TestLinearSolver.cpp
@@ -166,8 +166,7 @@ void checkLinearSolverInterface(T_MATRIX& A, T_VECTOR& b,
local_vec[1] = 2. * (mrank+1);
x.set(row_pos, local_vec);
- double x0[6];
- double x1[6];
+ std::vector<double> x0(6);
x.getGlobalVector(x0);
MathLib::LinAlg::matMult(A, x, b);
@@ -194,6 +193,7 @@ void checkLinearSolverInterface(T_MATRIX& A, T_VECTOR& b,
EXPECT_GT(ls.getNumberOfIterations(), 0u);
+ std::vector<double> x1(6);
x.getGlobalVector(x1);
ASSERT_ARRAY_NEAR(x0, x1, 6, 1e-5);
}