From 40d69743a19d3e08fbbb59f9b30708094638e5aa Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Fri, 5 May 2017 17:32:54 +0200
Subject: [PATCH] [PETsc] Disabled VecGetValues by using a global scatter
 vector

---
 MathLib/LinAlg/PETSc/PETScVector.cpp | 46 +++++++++++++++++++++++++++-
 MathLib/LinAlg/PETSc/PETScVector.h   | 41 ++++++++++++++++---------
 2 files changed, 72 insertions(+), 15 deletions(-)

diff --git a/MathLib/LinAlg/PETSc/PETScVector.cpp b/MathLib/LinAlg/PETSc/PETScVector.cpp
index e5ea7875ce0..6bbf6d41c93 100644
--- a/MathLib/LinAlg/PETSc/PETScVector.cpp
+++ b/MathLib/LinAlg/PETSc/PETScVector.cpp
@@ -38,6 +38,8 @@ PETScVector::PETScVector(const PetscInt vec_size, const bool is_global_size)
     }
 
     config();
+
+    _global_v = std::unique_ptr<PetscScalar[]>{ new PetscScalar[_size] };
 }
 
 PETScVector::PETScVector(const PetscInt vec_size,
@@ -61,6 +63,8 @@ PETScVector::PETScVector(const PetscInt vec_size,
     }
 
     config();
+
+    _global_v = std::unique_ptr<PetscScalar[]>{ new PetscScalar[_size] };
 }
 
 PETScVector::PETScVector(const PETScVector &existing_vec, const bool deep_copy)
@@ -72,6 +76,13 @@ PETScVector::PETScVector(const PETScVector &existing_vec, const bool deep_copy)
     {
         VecCopy(existing_vec._v, _v);
     }
+
+    if (!_global_v)
+    {
+        _global_v = std::unique_ptr<PetscScalar[]>{ new PetscScalar[_size] };
+    }
+
+    getGlobalVector(_global_v.get());
 }
 
 PETScVector::PETScVector(PETScVector &&other)
@@ -83,7 +94,14 @@ PETScVector::PETScVector(PETScVector &&other)
     , _size_loc{other._size_loc}
     , _size_ghosts{other._size_ghosts}
     , _has_ghost_id{other._has_ghost_id}
-{}
+{
+    if (!_global_v)
+    {
+        _global_v = std::unique_ptr<PetscScalar[]>{ new PetscScalar[_size] };
+    }
+
+    getGlobalVector(_global_v.get());
+}
 
 void PETScVector::config()
 {
@@ -101,6 +119,13 @@ void PETScVector::finalizeAssembly()
 {
     VecAssemblyBegin(_v);
     VecAssemblyEnd(_v);
+
+    if (!_global_v)
+    {
+        _global_v = std::unique_ptr<PetscScalar[]>{ new PetscScalar[_size] };
+    }
+
+    getGlobalVector(_global_v.get());
 }
 
 void PETScVector::gatherLocalVectors( PetscScalar local_array[],
@@ -165,6 +190,24 @@ void PETScVector::copyValues(std::vector<double>& u) const
     restoreArray(loc_x);
 }
 
+std::vector<double> PETScVector::get(std::vector<IndexType> const& indices) const
+{
+    std::vector<double> local_x(indices.size());
+    // If VecGetValues can get values from different processors,
+    // use VecGetValues(_v, indices.size(), indices.data(),
+    //                    local_x.data());
+
+    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] = _global_v[id_p];
+    }
+
+    return local_x;
+}
+
 PetscScalar* PETScVector::getLocalVector() const
 {
     PetscScalar *loc_array;
@@ -228,6 +271,7 @@ void PETScVector::shallowCopy(const PETScVector &v)
 void finalizeVectorAssembly(PETScVector &vec)
 {
     vec.finalizeAssembly();
+    vec.setGlobalVector();
 }
 
 } //end of namespace
diff --git a/MathLib/LinAlg/PETSc/PETScVector.h b/MathLib/LinAlg/PETSc/PETScVector.h
index cba4692313e..0ea4cb51e8f 100644
--- a/MathLib/LinAlg/PETSc/PETScVector.h
+++ b/MathLib/LinAlg/PETSc/PETScVector.h
@@ -16,6 +16,7 @@
 
 #pragma once
 
+#include <memory>
 #include <string>
 #include <vector>
 #include <petscvec.h>
@@ -155,21 +156,14 @@ class PETScVector
         }
 
         //! Get several entries
-        std::vector<double> get(std::vector<IndexType> const& indices) const
-        {
-            std::vector<double> local_x(indices.size());
-
-            VecGetValues(_v, indices.size(), indices.data(), local_x.data());
-
-            return local_x;
-        }
+        std::vector<double> get(std::vector<IndexType> const& indices) const;
 
         // TODO preliminary
         double operator[] (PetscInt idx) const
         {
-            double value;
-            VecGetValues(_v, 1, &idx, &value);
-            return value;
+            const PetscInt id_p = (idx == -_size) ?  0 : std::abs(idx);
+            //VecGetValues(_v, 1, &id_p, &value);
+            return _global_v[id_p];
         }
 
         /*!
@@ -188,9 +182,11 @@ class PETScVector
         /// 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);
-            return x;
+            //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];
         }
 
         // TODO eliminate in favour of getRawVector()
@@ -271,6 +267,13 @@ class PETScVector
         /// Flag to indicate whether the vector is created with ghost entry indices
         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.
+        */
+        std::unique_ptr<PetscScalar[]> _global_v = nullptr;
+
         /*!
               \brief  Collect local vectors
               \param  local_array Local array
@@ -293,6 +296,16 @@ class PETScVector
         /// A funtion called by constructors to configure members
         void config();
 
+        /// Set _global_v
+        void setGlobalVector()
+        {
+            if (!_global_v)
+            {
+                _global_v = std::unique_ptr<PetscScalar[]>{ new PetscScalar[_size] };;
+            }
+            getGlobalVector(_global_v.get());
+        }
+
         friend void finalizeVectorAssembly(PETScVector &vec);
 };
 
-- 
GitLab