From: Peter Munch <peterrmuench@gmail.com>
Date: Sun, 22 Mar 2020 14:02:01 +0000 (+0100)
Subject: Fill many functions
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=090b8b09dc56dc5dbc37d9f30f5ac51ca708398f;p=dealii.git

Fill many functions
---

diff --git a/include/deal.II/lac/la_sm_vector.h b/include/deal.II/lac/la_sm_vector.h
index d911f3a532..87ba0660fa 100644
--- a/include/deal.II/lac/la_sm_vector.h
+++ b/include/deal.II/lac/la_sm_vector.h
@@ -831,7 +831,6 @@ namespace LinearAlgebra
     inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
     Vector<Number, MemorySpace>::get_partitioner() const
     {
-      Assert(false, ExcNotImplemented());
       return partitioner;
     }
 
diff --git a/include/deal.II/lac/la_sm_vector.templates.h b/include/deal.II/lac/la_sm_vector.templates.h
index 54a8697c28..44dfbea92f 100644
--- a/include/deal.II/lac/la_sm_vector.templates.h
+++ b/include/deal.II/lac/la_sm_vector.templates.h
@@ -318,9 +318,11 @@ namespace LinearAlgebra
     Vector<Number, MemorySpaceType>::
     operator=(const Vector<Number, MemorySpaceType> &c)
     {
-      Assert(false, ExcNotImplemented());
-      (void)c;
-      return *this;
+#ifdef _MSC_VER
+      return this->operator=<Number>(c);
+#else
+      return this->template operator=<Number>(c);
+#endif
     }
 
 
@@ -331,8 +333,58 @@ namespace LinearAlgebra
     Vector<Number, MemorySpaceType>::
     operator=(const Vector<Number2, MemorySpaceType> &c)
     {
-      Assert(false, ExcNotImplemented());
-      (void)c;
+      Assert(c.partitioner.get() != nullptr, ExcNotInitialized());
+
+      // we update ghost values whenever one of the input or output vector
+      // already held ghost values or when we import data from a vector with
+      // the same local range but different ghost layout
+      bool must_update_ghost_values = c.vector_is_ghosted;
+
+      // check whether the two vectors use the same parallel partitioner. if
+      // not, check if all local ranges are the same (that way, we can
+      // exchange data between different parallel layouts). One variant which
+      // is included here and necessary for compatibility with the other
+      // distributed vector classes (Trilinos, PETSc) is the case when vector
+      // c does not have any ghosts (constructed without ghost elements given)
+      // but the current vector does: In that case, we need to exchange data
+      // also when none of the two vector had updated its ghost values before.
+      if (partitioner.get() == nullptr)
+        reinit(c, true);
+      else if (partitioner.get() != c.partitioner.get())
+        {
+          // local ranges are also the same if both partitioners are empty
+          // (even if they happen to define the empty range as [0,0) or [c,c)
+          // for some c!=0 in a different way).
+          int local_ranges_are_identical =
+            (partitioner->local_range() == c.partitioner->local_range() ||
+             (partitioner->local_range().second ==
+                partitioner->local_range().first &&
+              c.partitioner->local_range().second ==
+                c.partitioner->local_range().first));
+          if ((c.partitioner->n_mpi_processes() > 1 &&
+               Utilities::MPI::min(local_ranges_are_identical,
+                                   c.partitioner->get_mpi_communicator()) ==
+                 0) ||
+              !local_ranges_are_identical)
+            reinit(c, true);
+          else
+            must_update_ghost_values |= vector_is_ghosted;
+
+          must_update_ghost_values |=
+            (c.partitioner->ghost_indices_initialized() == false &&
+             partitioner->ghost_indices_initialized() == true);
+        }
+      else
+        must_update_ghost_values |= vector_is_ghosted;
+
+      const size_type this_size = partitioner->local_size();
+      if (this_size > 0)
+        std::memcpy(this->begin(), c.begin(), this_size * sizeof(Number));
+
+      if (must_update_ghost_values)
+        update_ghost_values();
+      else
+        zero_out_ghosts();
       return *this;
     }
 
@@ -369,7 +421,6 @@ namespace LinearAlgebra
     Vector<Number, MemorySpaceType>::compress(
       ::dealii::VectorOperation::values operation)
     {
-      Assert(false, ExcNotImplemented());
       compress_start(0, operation);
       compress_finish(operation);
     }
@@ -380,7 +431,6 @@ namespace LinearAlgebra
     void
     Vector<Number, MemorySpaceType>::update_ghost_values() const
     {
-      Assert(false, ExcNotImplemented());
       update_ghost_values_start();
       update_ghost_values_finish();
     }
@@ -721,8 +771,22 @@ namespace LinearAlgebra
     void
     Vector<Number, MemorySpaceType>::scale(const VectorSpaceVector<Number> &vv)
     {
-      Assert(false, ExcNotImplemented());
-      (void)vv;
+      // Downcast. Throws an exception if invalid.
+      using VectorType = Vector<Number, MemorySpaceType>;
+      Assert(dynamic_cast<const VectorType *>(&vv) != nullptr,
+             ExcVectorTypeNotCompatible());
+      const VectorType &v = dynamic_cast<const VectorType &>(vv);
+
+      AssertDimension(local_size(), v.local_size());
+
+      auto       values       = this->begin();
+      const auto values_other = v.begin();
+
+      for (unsigned int i = 0; i < partitioner->local_size(); i++)
+        values[i] *= values_other[i];
+
+      if (vector_is_ghosted)
+        update_ghost_values();
     }
 
 
@@ -732,9 +796,23 @@ namespace LinearAlgebra
     Vector<Number, MemorySpaceType>::equ(const Number                     a,
                                          const VectorSpaceVector<Number> &vv)
     {
-      Assert(false, ExcNotImplemented());
-      (void)a;
-      (void)vv;
+      // Downcast. Throws an exception if invalid.
+      using VectorType = Vector<Number, MemorySpaceType>;
+      Assert(dynamic_cast<const VectorType *>(&vv) != nullptr,
+             ExcVectorTypeNotCompatible());
+      const VectorType &v = dynamic_cast<const VectorType &>(vv);
+
+      AssertIsFinite(a);
+      AssertDimension(local_size(), v.local_size());
+
+      auto       values       = this->begin();
+      const auto values_other = v.begin();
+
+      for (unsigned int i = 0; i < partitioner->local_size(); i++)
+        values[i] = a * values_other[i];
+
+      if (vector_is_ghosted)
+        update_ghost_values();
     }
 
 
@@ -760,8 +838,7 @@ namespace LinearAlgebra
     bool
     Vector<Number, MemorySpaceType>::all_zero() const
     {
-      Assert(false, ExcNotImplemented());
-      return false;
+      return (linfty_norm() == 0) ? true : false;
     }
 
 
@@ -772,8 +849,20 @@ namespace LinearAlgebra
     Vector<Number, MemorySpaceType>::inner_product_local(
       const Vector<Number2, MemorySpaceType> &v) const
     {
-      Assert(false, ExcNotImplemented());
-      (void)v;
+      if (PointerComparison::equal(this, &v))
+        return norm_sqr_local();
+
+      AssertDimension(partitioner->local_size(), v.partitioner->local_size());
+
+      real_type sum = Number();
+
+      auto values       = this->begin();
+      auto values_other = v.begin();
+
+      for (unsigned int i = 0; i < partitioner->local_size(); i++)
+        sum += values[i] * values_other[i];
+
+      return sum;
     }
 
 
@@ -782,9 +871,18 @@ namespace LinearAlgebra
     Number Vector<Number, MemorySpaceType>::
            operator*(const VectorSpaceVector<Number> &vv) const
     {
-      Assert(false, ExcNotImplemented());
-      (void)vv;
-      return 0;
+      // Downcast. Throws an exception if invalid.
+      using VectorType = Vector<Number, MemorySpaceType>;
+      Assert((dynamic_cast<const VectorType *>(&vv) != nullptr),
+             ExcVectorTypeNotCompatible());
+      const VectorType &v = dynamic_cast<const VectorType &>(vv);
+
+      Number local_result = inner_product_local(v);
+      if (partitioner->n_mpi_processes() > 1)
+        return Utilities::MPI::sum(local_result,
+                                   partitioner->get_mpi_communicator());
+      else
+        return local_result;
     }
 
 
@@ -793,8 +891,16 @@ namespace LinearAlgebra
     typename Vector<Number, MemorySpaceType>::real_type
     Vector<Number, MemorySpaceType>::norm_sqr_local() const
     {
-      Assert(false, ExcNotImplemented());
-      return 0;
+      real_type sum = Number();
+
+      auto values = this->begin();
+
+      for (unsigned int i = 0; i < partitioner->local_size(); i++)
+        sum += values[i] * values[i];
+
+      AssertIsFinite(sum);
+
+      return sum;
     }
 
 
@@ -843,8 +949,12 @@ namespace LinearAlgebra
     typename Vector<Number, MemorySpaceType>::real_type
     Vector<Number, MemorySpaceType>::norm_sqr() const
     {
-      Assert(false, ExcNotImplemented());
-      return 0;
+      real_type local_result = norm_sqr_local();
+      if (partitioner->n_mpi_processes() > 1)
+        return Utilities::MPI::sum(local_result,
+                                   partitioner->get_mpi_communicator());
+      else
+        return local_result;
     }
 
 
@@ -853,8 +963,7 @@ namespace LinearAlgebra
     typename Vector<Number, MemorySpaceType>::real_type
     Vector<Number, MemorySpaceType>::l2_norm() const
     {
-      Assert(false, ExcNotImplemented());
-      return 0;
+      return std::sqrt(norm_sqr());
     }
 
 
@@ -918,12 +1027,26 @@ namespace LinearAlgebra
       const Vector<Number, MemorySpaceType> &v,
       const Vector<Number, MemorySpaceType> &w)
     {
-      Assert(false, ExcNotImplemented());
-      (void)a;
-      (void)v;
-      (void)w;
+      const size_type vec_size = partitioner->local_size();
+      AssertDimension(vec_size, v.local_size());
+      AssertDimension(vec_size, w.local_size());
 
-      return 0;
+      real_type sum = Number();
+
+      auto values   = this->begin();
+      auto values_v = v.begin();
+      auto values_w = w.begin();
+
+      for (unsigned int i = 0; i < partitioner->local_size(); i++)
+        {
+          const auto temp = values[i] + a * values_v[i];
+          values[i]       = temp;
+          sum += temp * values_w[i];
+        }
+
+      AssertIsFinite(sum);
+
+      return sum;
     }
 
 
@@ -935,11 +1058,21 @@ namespace LinearAlgebra
       const VectorSpaceVector<Number> &vv,
       const VectorSpaceVector<Number> &ww)
     {
-      Assert(false, ExcNotImplemented());
-      (void)a;
-      (void)vv;
-      (void)ww;
-      return 0;
+      // Downcast. Throws an exception if invalid.
+      using VectorType = Vector<Number, MemorySpaceType>;
+      Assert((dynamic_cast<const VectorType *>(&vv) != nullptr),
+             ExcVectorTypeNotCompatible());
+      const VectorType &v = dynamic_cast<const VectorType &>(vv);
+      Assert((dynamic_cast<const VectorType *>(&ww) != nullptr),
+             ExcVectorTypeNotCompatible());
+      const VectorType &w = dynamic_cast<const VectorType &>(ww);
+
+      Number local_result = add_and_dot_local(a, v, w);
+      if (partitioner->n_mpi_processes() > 1)
+        return Utilities::MPI::sum(local_result,
+                                   partitioner->get_mpi_communicator());
+      else
+        return local_result;
     }
 
 
diff --git a/source/lac/la_sm_vector.cc b/source/lac/la_sm_vector.cc
index 1334c456f2..acfab48a20 100644
--- a/source/lac/la_sm_vector.cc
+++ b/source/lac/la_sm_vector.cc
@@ -34,8 +34,9 @@ namespace LinearAlgebra
     &Vector<S1, ::dealii::MemorySpace::Host>::operator= \
       <S2>(const Vector<S2, ::dealii::MemorySpace::Host> &)
 
-    TEMPL_COPY_CONSTRUCTOR(double, float);
-    TEMPL_COPY_CONSTRUCTOR(float, double);
+    // TODO
+    // TEMPL_COPY_CONSTRUCTOR(double, float);
+    // TEMPL_COPY_CONSTRUCTOR(float, double);
 
 #undef TEMPL_COPY_CONSTRUCTOR
   } // namespace SharedMPI
diff --git a/source/lac/vector_memory.cc b/source/lac/vector_memory.cc
index 461c1cf50c..0cc491dce8 100644
--- a/source/lac/vector_memory.cc
+++ b/source/lac/vector_memory.cc
@@ -16,6 +16,7 @@
 #include <deal.II/lac/block_vector.h>
 #include <deal.II/lac/la_parallel_block_vector.h>
 #include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/la_sm_vector.h>
 #include <deal.II/lac/la_vector.h>
 #include <deal.II/lac/petsc_block_vector.h>
 #include <deal.II/lac/petsc_vector.h>
@@ -31,6 +32,8 @@ DEAL_II_NAMESPACE_OPEN
 
 #include "vector_memory.inst"
 
+template class GrowingVectorMemory<LinearAlgebra::SharedMPI::Vector<double>>;
+
 namespace internal
 {
   namespace GrowingVectorMemoryImplementation