Align the underlying array of Vector and add runtime switched between add and vectori...

diff --git a/include/deal.II/lac/vector.h b/include/deal.II/lac/vector.h
index 59205f0dc252c87d55d094aea6cf16973af66c9a..36076effb1ede0b00676d82c334e773af7704dc3 100644 (file)
--- a/include/deal.II/lac/vector.h
+++ b/include/deal.II/lac/vector.h
@@ -31,6 +31,10 @@
 #include <cstring>
 #include <vector>
 
+#if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0
+#include <emmintrin.h>
+#endif
+
 DEAL_II_NAMESPACE_OPEN
 
 
@@ -626,14 +630,6 @@ public:
    */
   void add (const Vector<Number> &V);
 
-  /**
-   * Simple vectorized  (SIMD) vector addition, equal to the <tt>operator +=</tt>. 
-   * This function should not be used to add the vector to itself, i.e to
-   * multiply the current vector by two.
-   *
-   * @dealiiOperationIsMultithreaded
-   */
-  void vectorized_add (const Vector<Number> &V);
 
   /**
    * Multiple addition of scaled vectors, i.e. <tt>*this += a*V+b*W</tt>.
@@ -962,6 +958,18 @@ protected:
    * VectorView will access the pointer.
    */
   friend class VectorView<Number>;
+
+private:
+
+  /**
+   * Allocate @p v and, if possible, align along 64-byte boundaries.
+   */
+  void initialize_val(const size_type n);
+
+  /**
+   * Free @p val.
+   */
+  void clear_val();
 };
 
 /*@}*/
@@ -1017,7 +1025,7 @@ Vector<Number>::~Vector ()
 {
   if (val)
     {
-      delete[] val;
+      clear_val();
       val=0;
     }
 }
@@ -1030,7 +1038,7 @@ void Vector<Number>::reinit (const size_type n, const bool fast)
 {
   if (n==0)
     {
-      if (val) delete[] val;
+      if (val) clear_val();
       val = 0;
       max_vec_size = vec_size = 0;
       return;
@@ -1038,8 +1046,8 @@ void Vector<Number>::reinit (const size_type n, const bool fast)
 
   if (n>max_vec_size)
     {
-      if (val) delete[] val;
-      val = new value_type[n];
+      if (val) clear_val();
+      initialize_val(n);
       Assert (val != 0, ExcOutOfMemory());
       max_vec_size = n;
     };
@@ -1341,11 +1349,38 @@ Vector<Number>::load (Archive &ar, const unsigned int)
 
   ar &vec_size &max_vec_size ;
 
-  val = new Number[max_vec_size];
+  initialize_val(max_vec_size);
   ar &boost::serialization::make_array(val, max_vec_size);
 }
 
 
+
+template <typename Number>
+inline 
+void 
+Vector<Number>::initialize_val(const size_type size)
+{
+#if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0
+  val = static_cast<Number*>(_mm_malloc (size, 64));
+#else
+  val = static_cast<Number*>(malloc (size));
+#endif
+}
+
+
+
+template <typename Number>
+inline
+void
+Vector<Number>::clear_val()
+{
+#if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0
+  _mm_free(val);
+#else
+  free(val);
+#endif
+}
+
 #endif
 
 

diff --git a/include/deal.II/lac/vector.templates.h b/include/deal.II/lac/vector.templates.h
index 0289f655b73e0082bcfa6783d4bff0030790e209..ec034a507e8aa673b835ba73f68d644952b3339e 100644 (file)
--- a/include/deal.II/lac/vector.templates.h
+++ b/include/deal.II/lac/vector.templates.h
@@ -148,7 +148,7 @@ Vector<Number>::Vector (const Vector<Number> &v)
 {
   if (vec_size != 0)
     {
-      val = new Number[max_vec_size];
+      initialize_val(max_vec_size);
       Assert (val != 0, ExcOutOfMemory());
       *this = v;
     }
@@ -168,7 +168,7 @@ Vector<Number>::Vector (const Vector<OtherNumber> &v)
 {
   if (vec_size != 0)
     {
-      val = new Number[max_vec_size];
+      initialize_val(max_vec_size);
       Assert (val != 0, ExcOutOfMemory());
       std::copy (v.begin(), v.end(), begin());
     }
@@ -189,7 +189,7 @@ Vector<Number>::Vector (const PETScWrappers::Vector &v)
 {
   if (vec_size != 0)
     {
-      val = new Number[max_vec_size];
+      iniatilize_val(max_vec_size);
       Assert (val != 0, ExcOutOfMemory());
 
       // get a representation of the vector
@@ -241,7 +241,7 @@ Vector<Number>::Vector (const TrilinosWrappers::MPI::Vector &v)
 {
   if (vec_size != 0)
     {
-      val = new Number[max_vec_size];
+      initialize_val(max_vec_size);
       Assert (val != 0, ExcOutOfMemory());
 
       // Copy the distributed vector to
@@ -275,7 +275,7 @@ Vector<Number>::Vector (const TrilinosWrappers::Vector &v)
 {
   if (vec_size != 0)
     {
-      val = new Number[max_vec_size];
+      initialize_val(max_vec_size);
       Assert (val != 0, ExcOutOfMemory());
 
       // get a representation of the vector
@@ -972,16 +972,45 @@ void Vector<Number>::add (const Vector<Number> &v)
   Assert (vec_size!=0, ExcEmptyObject());
   Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
 
-  if (vec_size>internal::Vector::minimum_parallel_grain_size)
-    parallel::transform (val,
-                         val+vec_size,
-                         v.val,
-                         val,
-                         (boost::lambda::_1 + boost::lambda::_2),
-                         internal::Vector::minimum_parallel_grain_size);
-  else if (vec_size > 0)
+  // Vectorization can only be if v is not the current vector otherwise there
+  // will be collision
+  if (&(v.val[0])!=&val[0])
+  {
+#ifndef DEAL_II_WITH_THREADS
+#pragma omp simd
     for (size_type i=0; i<vec_size; ++i)
       val[i] += v.val[i];
+#else
+    if (vec_size>internal::Vector::minimum_parallel_grain_size)
+    {
+      internal::Vectorization_add_1<Number> vector_add;
+      vector_add.val = val;
+      vector_add.v_val = v.val;
+      tbb::parallel_for (tbb::blocked_range<size_type> (0, 
+            vec_size, 
+            internal::Vector::minimum_parallel_grain_size),
+          vector_add,
+          tbb::auto_partitioner());
+    }
+    else if (vec_size > 0)
+#pragma omp simd
+      for (size_type i=0; i<vec_size; ++i)
+        val[i] += v.val[i];
+#endif
+  }
+  else
+  {
+    if (vec_size>internal::Vector::minimum_parallel_grain_size)
+      parallel::transform (val,
+          val+vec_size,
+          v.val,
+          val,
+          (boost::lambda::_1 + boost::lambda::_2),
+          internal::Vector::minimum_parallel_grain_size);
+    else if (vec_size > 0)
+      for (size_type i=0; i<vec_size; ++i)
+        val[i] += v.val[i];
+  }
 }
 
 
@@ -1010,35 +1039,6 @@ void Vector<Number>::add (const Number a, const Vector<Number> &v,
 }
 
 
-template <typename Number>
-void Vector<Number>::vectorized_add (const Vector<Number> &v)
-{
-  Assert (vec_size!=0, ExcEmptyObject());
-  Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
-
-#ifndef DEAL_II_WITH_THREADS
-#pragma omp simd
-  for (size_type i=0; i<vec_size; ++i)
-    val[i] += v.val[i];
-#else
-  if (vec_size>internal::Vector::minimum_parallel_grain_size)
-  {
-    internal::Vectorization_add_1<Number> vector_add;
-    vector_add.val = val;
-    vector_add.v_val = v.val;
-    tbb::parallel_for (tbb::blocked_range<size_type> (0, 
-                                                      vec_size, 
-                                                      internal::Vector::minimum_parallel_grain_size),
-                       vector_add,
-                       tbb::auto_partitioner());
-  }
-  else if (vec_size > 0)
-#pragma omp simd
-    for (size_type i=0; i<vec_size; ++i)
-      val[i] += v.val[i];
-#endif
-}
-
 
 template <typename Number>
 void Vector<Number>::sadd (const Number x,