[PETSc] Added a functionality for accessing entries via local array

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);
 }

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

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,

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);
 }