Use a smarter implementation of a parallel = operator. Work around a problem with...

git-svn-id: https://svn.dealii.org/trunk@19413 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/lac/include/lac/vector.h b/deal.II/lac/include/lac/vector.h
index ad1ad3e760ef31c67d5b83f339b087f067eb4319..38ed5b720cf732534d809b1d0f4271669383d783 100644 (file)
--- a/deal.II/lac/include/lac/vector.h
+++ b/deal.II/lac/include/lac/vector.h
@@ -1015,6 +1015,49 @@ void Vector<Number>::reinit (const unsigned int n, const bool fast)
 
 
 
+namespace internal
+{
+  namespace Vector
+  {
+    template<typename T>
+    void set_subrange (const T            s,
+                      const unsigned int begin,
+                      const unsigned int end,
+                      dealii::Vector<T> &dst)
+    {
+      if (s == T())
+       memset ((dst.begin()+begin),0,(end-begin)*sizeof(T));
+      else
+       std::fill (&*(dst.begin()+begin), &*(dst.begin()+end), s);
+    }
+
+    template<typename T>
+    void copy_subrange (const dealii::Vector<T>&src,
+                       const unsigned int      begin,
+                       const unsigned int      end,
+                       dealii::Vector<T>      &dst)
+    {
+      memcpy(&*(dst.begin()+begin), &*(src.begin()+begin),
+            (end-begin)*sizeof(T));
+    }
+
+    template<typename T, typename U>
+    void copy_subrange_ext (const dealii::Vector<T>&src,
+                           const unsigned int      begin,
+                           const unsigned int      end,
+                           dealii::Vector<U>      &dst)
+    {
+      const T* q = src.begin()+begin;
+      const T* const end_q = src.begin()+end;
+      U* p = dst.begin()+begin;
+      for (; q!=end_q; ++q, ++p)
+       *p = *q;
+    }
+  }
+}
+
+
+
 template <typename Number>
 inline
 Vector<Number> & Vector<Number>::operator = (const Number s)
@@ -1025,10 +1068,70 @@ Vector<Number> & Vector<Number>::operator = (const Number s)
   if (s != Number())
     Assert (vec_size!=0, ExcEmptyObject());
   if (vec_size!=0)
-    parallel::apply_to_subranges (begin(), end(),
-                                 std_cxx1x::bind(&std::fill<Number*,Number>,
-                                                 _1, _2, s),
+    parallel::apply_to_subranges (0U, size(),
+                                 std_cxx1x::bind(&internal::Vector::template
+                                                 set_subrange<Number>,
+                                                 s, _1, _2, std_cxx1x::ref(*this)),
+                                 internal::Vector::minimum_parallel_grain_size);
+
+  return *this;
+}
+
+
+
+template <>
+inline
+Vector<std::complex<float> > & Vector<std::complex<float> >::operator = (const std::complex<float> s)
+{
+  Assert (numbers::is_finite(s),
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
+  if (s != std::complex<float>())
+    Assert (vec_size!=0, ExcEmptyObject());
+  if (vec_size!=0)
+    std::fill (begin(), end(), s);
+
+  return *this;
+}
+
+
+
+template <typename Number>
+inline
+Vector<Number> &
+Vector<Number>::operator = (const Vector<Number>& v)
+{
+  if (v.vec_size != vec_size)
+    reinit (v.vec_size, true);
+  if (vec_size!=0)
+    parallel::apply_to_subranges (0U, size(),
+                                 std_cxx1x::bind(&internal::Vector::template
+                                                 copy_subrange<Number>,
+                                                 std_cxx1x::cref(v), _1, _2, 
+                                                 std_cxx1x::ref(*this)),
+                                 internal::Vector::minimum_parallel_grain_size);
+
+  return *this;
+}
+
+
+
+template <typename Number>
+template <typename Number2>
+inline
+Vector<Number> &
+Vector<Number>::operator = (const Vector<Number2>& v)
+{
+  if (v.size() != vec_size)
+    reinit (v.size(), true);
+  if (vec_size!=0)
+    parallel::apply_to_subranges (0U, size(),
+                                 std_cxx1x::bind(&internal::Vector::template
+                                                 copy_subrange_ext<Number2,Number>,
+                                                 std_cxx1x::cref(v), _1, _2, 
+                                                 std_cxx1x::ref(*this)),
                                  internal::Vector::minimum_parallel_grain_size);
+  
   return *this;
 }
 

diff --git a/deal.II/lac/include/lac/vector.templates.h b/deal.II/lac/include/lac/vector.templates.h
index 6aafe6bc39425e2841a7afb8e9535e04762c29ab..c1d9f089ee5c0c6a2e11c0544d29f7e560ced8d9 100644 (file)
--- a/deal.II/lac/include/lac/vector.templates.h
+++ b/deal.II/lac/include/lac/vector.templates.h
@@ -666,39 +666,6 @@ void Vector<Number>::ratio (const Vector<Number> &a,
 
 
 
-template <typename Number>
-Vector<Number> &
-Vector<Number>::operator = (const Vector<Number>& v)
-{
-  if (v.vec_size != vec_size)
-    reinit (v.vec_size, true);
-  if (vec_size!=0)
-    parallel::transform (v.val,
-                        v.val+v.vec_size,
-                        val,
-                        (boost::lambda::_1),
-                        internal::Vector::minimum_parallel_grain_size);
-  
-  return *this;
-}
-
-
-
-template <typename Number>
-template <typename Number2>
-Vector<Number> &
-Vector<Number>::operator = (const Vector<Number2>& v)
-{
-  if (v.size() != vec_size)
-    reinit (v.size(), true);
-  if (vec_size!=0)
-    std::copy (v.begin(), v.end(), begin());
-  
-  return *this;
-}
-
-
-
 template <typename Number>
 Vector<Number> &
 Vector<Number>::operator = (const BlockVector<Number>& v)

diff --git a/deal.II/lac/source/vector.cc b/deal.II/lac/source/vector.cc
index 39607ff73ab920dd587876524b6f22ce2297fa73..0ce029e10f86628fc118769dd669d6c2113385a0 100644 (file)
--- a/deal.II/lac/source/vector.cc
+++ b/deal.II/lac/source/vector.cc
@@ -52,10 +52,6 @@ TEMPL_COPY_CONSTRUCTOR(std::complex<float>,std::complex<long double>);
 
 
 #define TEMPL_OP_EQ(S1,S2) \
-  template \
-  Vector<S1>& \
-  Vector<S1>::DEAL_II_MEMBER_OP_TEMPLATE_INST \
-  operator=<S2>(const Vector<S2>&);                  \
   template void Vector<S1>::scale (const Vector<S2>&); \
   template void Vector<S1>::equ (const S1, const Vector<S2>&)