Provide parallelized versions of some more vector-vector operations like assignments...

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

diff --git a/deal.II/lac/include/lac/vector.h b/deal.II/lac/include/lac/vector.h
index f3ec98636754cdf30916d133f19130557e25580f..7c9abe9647783f8070746db688de2e4bcda6634e 100644 (file)
--- a/deal.II/lac/include/lac/vector.h
+++ b/deal.II/lac/include/lac/vector.h
@@ -16,7 +16,9 @@
 
 #include <base/config.h>
 #include <base/exceptions.h>
+#include <base/parallel.h>
 #include <base/subscriptor.h>
+#include <boost/lambda/lambda.hpp>
 
 #include <cstdio>
 
@@ -767,6 +769,17 @@ class Vector : public Subscriptor
                                      */
     Vector<Number> & operator /= (const Number factor);
     
+                                    /**
+                                     * Scale each element of this
+                                     * vector by the corresponding
+                                     * element in the argument. This
+                                     * function is mostly meant to
+                                     * simulate multiplication (and
+                                     * immediate re-assignment) by a
+                                     * diagonal scaling matrix.
+                                     */
+    void scale (const Vector<Number> &scaling_factors);
+    
                                     /**
                                      * Scale each element of this
                                      * vector by the corresponding
@@ -778,7 +791,12 @@ class Vector : public Subscriptor
                                      */
     template <typename Number2>
     void scale (const Vector<Number2> &scaling_factors);
-    
+
+                                    /**
+                                     * Assignment <tt>*this = a*u</tt>.
+                                     */
+    void equ (const Number a, const Vector<Number>& u);
+
                                     /**
                                      * Assignment <tt>*this = a*u</tt>.
                                      */
@@ -1033,7 +1051,10 @@ Vector<Number> & Vector<Number>::operator = (const Number s)
   if (s != Number())
     Assert (vec_size!=0, ExcEmptyObject());
   if (vec_size!=0)
-    std::fill (begin(), end(), s);
+    parallel::apply_to_subranges (begin(), end(), 
+                                 std_cxx1x::bind(&std::fill<Number*,Number>,
+                                                 _1, _2, s),
+                                 internal::Vector::minimum_parallel_grain_size);
   return *this;
 }
 
@@ -1147,10 +1168,11 @@ Vector<Number>::scale (const Number factor)
 
   Assert (vec_size!=0, ExcEmptyObject());
 
-  iterator             ptr  = begin();
-  const const_iterator eptr = end();
-  while (ptr!=eptr)
-    *ptr++ *= factor;
+  parallel::transform (val,
+                      val+vec_size,
+                      val,
+                      (factor*boost::lambda::_1),
+                      internal::Vector::minimum_parallel_grain_size);
 }
 
 
@@ -1167,11 +1189,12 @@ Vector<Number>::add (const Number a,
   Assert (vec_size!=0, ExcEmptyObject());
   Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
 
-  iterator i_ptr = begin(),
-          i_end = end();
-  const_iterator v_ptr = v.begin();
-  while (i_ptr!=i_end)
-    *i_ptr++ += a * *v_ptr++;
+  parallel::transform (val,
+                      val+vec_size,
+                      v.val,
+                      val,
+                      (boost::lambda::_1 + a*boost::lambda::_2),
+                      internal::Vector::minimum_parallel_grain_size);
 }
 
 
@@ -1190,11 +1213,13 @@ Vector<Number>::sadd (const Number x,
 
   Assert (vec_size!=0, ExcEmptyObject());
   Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
-  iterator i_ptr = begin(),
-          i_end = end();
-  const_iterator v_ptr = v.begin();
-  for (; i_ptr!=i_end; ++i_ptr)
-    *i_ptr = x * *i_ptr  +  a * *v_ptr++;
+
+  parallel::transform (val,
+                      val+vec_size,
+                      v.val,
+                      val,
+                      (x*boost::lambda::_1 + a*boost::lambda::_2),
+                      internal::Vector::minimum_parallel_grain_size);
 }
 
 

diff --git a/deal.II/lac/include/lac/vector.templates.h b/deal.II/lac/include/lac/vector.templates.h
index 5916637e5715fe6d5d66c73ad5a10eafe6f95246..6aafe6bc39425e2841a7afb8e9535e04762c29ab 100644 (file)
--- a/deal.II/lac/include/lac/vector.templates.h
+++ b/deal.II/lac/include/lac/vector.templates.h
@@ -16,7 +16,6 @@
 
 #include <base/template_constraints.h>
 #include <base/numbers.h>
-#include <base/parallel.h>
 #include <lac/vector.h>
 #include <lac/block_vector.h>
 
@@ -29,8 +28,6 @@
 #  include <lac/trilinos_vector.h>
 #endif
 
-#include <boost/lambda/lambda.hpp>
-
 #include <cmath>
 #include <cstring>
 #include <algorithm>
@@ -528,22 +525,24 @@ void Vector<Number>::sadd (const Number x, const Number a,
                           const Vector<Number>& w, const Number c,
                            const Vector<Number>& y)
 {
-  Assert (numbers::is_finite(x), 
-          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
-  Assert (numbers::is_finite(a), 
-          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
-  Assert (numbers::is_finite(b), 
-          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
-  Assert (numbers::is_finite(c), 
-          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+  sadd (x, a, v, b, w);
+  add (c, y);
+}
+
+
 
+template <typename Number>
+void Vector<Number>::scale (const Vector<Number> &s)
+{
   Assert (vec_size!=0, ExcEmptyObject());
-  Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
-  Assert (vec_size == w.vec_size, ExcDimensionMismatch(vec_size, w.vec_size));
-  Assert (vec_size == y.vec_size, ExcDimensionMismatch(vec_size, y.vec_size));
+  Assert (vec_size == s.vec_size, ExcDimensionMismatch(vec_size, s.vec_size));
 
-  for (unsigned int i=0; i<vec_size; ++i)
-    val[i] = x*val[i] + a*v.val[i] + b*w.val[i] +c*y.val[i];
+  parallel::transform (val,
+                      val+vec_size,
+                      s.val,
+                      val,
+                      (boost::lambda::_1*boost::lambda::_2),
+                      internal::Vector::minimum_parallel_grain_size);
 }
 
 
@@ -554,12 +553,32 @@ void Vector<Number>::scale (const Vector<Number2> &s)
 {
   Assert (vec_size!=0, ExcEmptyObject());
   Assert (vec_size == s.vec_size, ExcDimensionMismatch(vec_size, s.vec_size));
-  
+
   for (unsigned int i=0; i<vec_size; ++i)
     val[i] *= s.val[i];
 }
 
 
+
+template <typename Number>
+void Vector<Number>::equ (const Number a,
+                         const Vector<Number>& u)
+{
+  Assert (numbers::is_finite(a), 
+         ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
+  Assert (vec_size!=0, ExcEmptyObject());
+  Assert (vec_size == u.vec_size, ExcDimensionMismatch(vec_size, u.vec_size));
+
+  parallel::transform (u.val,
+                      u.val+u.vec_size,
+                      val,
+                      (a*boost::lambda::_1),
+                      internal::Vector::minimum_parallel_grain_size);
+}
+
+
+
 template <typename Number>
 template <typename Number2>
 void Vector<Number>::equ (const Number a,
@@ -582,6 +601,7 @@ void Vector<Number>::equ (const Number a,
 }
 
 
+
 template <typename Number>
 void Vector<Number>::equ (const Number a, const Vector<Number>& u,
                          const Number b, const Vector<Number>& v)
@@ -595,9 +615,13 @@ void Vector<Number>::equ (const Number a, const Vector<Number>& u,
   Assert (vec_size == u.vec_size, ExcDimensionMismatch(vec_size, u.vec_size));
   Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
 
-  for (unsigned int i=0; i<vec_size; ++i)
-    val[i] = a * u.val[i] + b * v.val[i];
-}
+  parallel::transform (u.val,
+                      u.val+u.vec_size,
+                      v.val,
+                      val,
+                      (a*boost::lambda::_1 + b*boost::lambda::_2),
+                      internal::Vector::minimum_parallel_grain_size);
+  }
 
 
 template <typename Number>
@@ -610,8 +634,13 @@ void Vector<Number>::equ (const Number a, const Vector<Number>& u,
   Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
   Assert (vec_size == w.vec_size, ExcDimensionMismatch(vec_size, w.vec_size));
 
-  for (unsigned int i=0; i<vec_size; ++i)
-    val[i] = a * u.val[i] + b * v.val[i] + c * w.val[i];
+  parallel::transform (u.val,
+                      u.val+u.vec_size,
+                      v.val,
+                      w.val,
+                      val,
+                      (a*boost::lambda::_1 + b*boost::lambda::_2 + c*boost::lambda::_3),
+                      internal::Vector::minimum_parallel_grain_size);
 }
 
 
@@ -627,8 +656,12 @@ void Vector<Number>::ratio (const Vector<Number> &a,
                                   // we overwrite them anyway
   reinit (a.size(), true);
 
-  for (unsigned int i=0; i<vec_size; ++i)
-    val[i] = a.val[i] / b.val[i];
+  parallel::transform (a.val,
+                      a.val+a.vec_size,
+                      b.val,
+                      val,
+                      (boost::lambda::_1 / boost::lambda::_2),
+                      internal::Vector::minimum_parallel_grain_size);
 }
 
 
@@ -640,7 +673,11 @@ Vector<Number>::operator = (const Vector<Number>& v)
   if (v.vec_size != vec_size)
     reinit (v.vec_size, true);
   if (vec_size!=0)
-    std::copy (v.begin(), v.end(), begin());
+    parallel::transform (v.val,
+                        v.val+v.vec_size,
+                        val,
+                        (boost::lambda::_1),
+                        internal::Vector::minimum_parallel_grain_size);
   
   return *this;
 }

diff --git a/deal.II/lac/source/vector.cc b/deal.II/lac/source/vector.cc
index f37e3c8b64a477680bc85048eab9f5ada525391a..7a41a6c4268fd366dbd4547768a9c54df8916ed9 100644 (file)
--- a/deal.II/lac/source/vector.cc
+++ b/deal.II/lac/source/vector.cc
@@ -76,7 +76,9 @@ TEMPL_COPY_CONSTRUCTOR(std::complex<float>,std::complex<long double>);
   template \
   Vector<S1>& \
   Vector<S1>::DEAL_II_MEMBER_OP_TEMPLATE_INST \
-  operator=<S2>(const Vector<S2>&)
+  operator=<S2>(const Vector<S2>&);                  \
+  template void Vector<S1>::scale (const Vector<S2>&); \
+  template void Vector<S1>::equ (const S1, const Vector<S2>&)    
 
 TEMPL_OP_EQ(double,float);
 TEMPL_OP_EQ(float,double);

diff --git a/deal.II/lac/source/vector.inst.in b/deal.II/lac/source/vector.inst.in
index 39caf6945a36dd2feb694e2ae9f544e53bc50a64..77dea0609bd5a1368740de6416e9eb396d4a781d 100644 (file)
--- a/deal.II/lac/source/vector.inst.in
+++ b/deal.II/lac/source/vector.inst.in
@@ -29,10 +29,6 @@ for (S1, S2 : REAL_SCALARS)
       operator*<S2>(const Vector<S2>&) const;
     template
       void Vector<S1>::reinit<S2>(const Vector<S2>&, const bool);
-    template
-      void Vector<S1>::equ<S2>(const S1, const Vector<S2>&);
-    template
-      void Vector<S1>::scale<S2>(const Vector<S2>&);
   }
 
 
@@ -54,8 +50,4 @@ for (S1, S2 : COMPLEX_SCALARS)
       operator*<S2>(const Vector<S2>&) const;
     template
       void Vector<S1>::reinit<S2>(const Vector<S2>&, const bool);
-    template
-      void Vector<S1>::equ<S2>(const S1, const Vector<S2>&);
-    template
-      void Vector<S1>::scale<S2>(const Vector<S2>&);
   }