From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Fri, 21 Aug 2009 15:39:07 +0000 (+0000)
Subject: Provide parallelized versions of some more vector-vector operations like assignments... 
X-Git-Tag: v8.0.0~7241
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9954a53b8e5c55b8ec92c9f1a557812afbad9337;p=dealii.git

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

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 f3ec986367..7c9abe9647 100644
--- 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 5916637e57..6aafe6bc39 100644
--- 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 f37e3c8b64..7a41a6c426 100644
--- 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 39caf6945a..77dea0609b 100644
--- 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>&);
   }