From 1e675f2c2ab5a4ff47105d9026a70a43a2f3e52f Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 18 Apr 2016 16:26:08 +0200 Subject: [PATCH] Increment index at the position where we do the summation --- include/deal.II/lac/vector.templates.h | 115 +++++++++++-------------- 1 file changed, 52 insertions(+), 63 deletions(-) diff --git a/include/deal.II/lac/vector.templates.h b/include/deal.II/lac/vector.templates.h index d1b4a7f3a4..283c9ab407 100644 --- a/include/deal.II/lac/vector.templates.h +++ b/include/deal.II/lac/vector.templates.h @@ -575,23 +575,20 @@ namespace internal struct Dot { static const bool vectorizes = types_are_equal::value && - (types_are_equal::value || - types_are_equal::value); + (VectorizedArray::n_array_elements > 1); Number - operator() (size_type &i) const + operator() (const size_type i) const { - Number xi = X[i]; - return xi * Number(numbers::NumberTraits::conjugate(Y[i++])); + return X[i] * Number(numbers::NumberTraits::conjugate(Y[i])); } VectorizedArray - do_vectorized(size_type &i) const + do_vectorized(const size_type i) const { VectorizedArray x, y; x.load(X+i); y.load(Y+i); - i += VectorizedArray::n_array_elements; return x * y; } @@ -602,21 +599,19 @@ namespace internal template struct Norm2 { - static const bool vectorizes = types_are_equal::value || - types_are_equal::value; + static const bool vectorizes = VectorizedArray::n_array_elements > 1; RealType - operator() (size_type &i) const + operator() (const size_type i) const { - return numbers::NumberTraits::abs_square(X[i++]); + return numbers::NumberTraits::abs_square(X[i]); } VectorizedArray - do_vectorized(size_type &i) const + do_vectorized(const size_type i) const { VectorizedArray x; x.load(X+i); - i += VectorizedArray::n_array_elements; return x * x; } @@ -626,21 +621,19 @@ namespace internal template struct Norm1 { - static const bool vectorizes = types_are_equal::value || - types_are_equal::value; + static const bool vectorizes = VectorizedArray::n_array_elements > 1; RealType - operator() (size_type &i) const + operator() (const size_type i) const { - return numbers::NumberTraits::abs(X[i++]); + return numbers::NumberTraits::abs(X[i]); } VectorizedArray - do_vectorized(size_type &i) const + do_vectorized(const size_type i) const { VectorizedArray x; x.load(X+i); - i += VectorizedArray::n_array_elements; return std::abs(x); } @@ -650,21 +643,19 @@ namespace internal template struct NormP { - static const bool vectorizes = types_are_equal::value || - types_are_equal::value; + static const bool vectorizes = VectorizedArray::n_array_elements > 1; RealType - operator() (size_type &i) const + operator() (const size_type i) const { - return std::pow(numbers::NumberTraits::abs(X[i++]), p); + return std::pow(numbers::NumberTraits::abs(X[i]), p); } VectorizedArray - do_vectorized(size_type &i) const + do_vectorized(const size_type i) const { VectorizedArray x; x.load(X+i); - i += VectorizedArray::n_array_elements; return std::pow(std::abs(x),p); } @@ -675,21 +666,19 @@ namespace internal template struct MeanValue { - static const bool vectorizes = types_are_equal::value || - types_are_equal::value; + static const bool vectorizes = VectorizedArray::n_array_elements > 1; Number - operator() (size_type &i) const + operator() (const size_type i) const { - return X[i++]; + return X[i]; } VectorizedArray - do_vectorized(size_type &i) const + do_vectorized(const size_type i) const { VectorizedArray x; x.load(X+i); - i += VectorizedArray::n_array_elements; return x; } @@ -699,19 +688,17 @@ namespace internal template struct AddAndDot { - static const bool vectorizes = types_are_equal::value || - types_are_equal::value; + static const bool vectorizes = VectorizedArray::n_array_elements > 1; Number - operator() (size_type &i) const + operator() (const size_type i) const { X[i] += a * V[i]; - Number xi = X[i]; - return xi * Number(numbers::NumberTraits::conjugate(W[i++])); + return X[i] * Number(numbers::NumberTraits::conjugate(W[i])); } VectorizedArray - do_vectorized(size_type &i) const + do_vectorized(const size_type i) const { VectorizedArray x, w, v; x.load(X+i); @@ -721,7 +708,6 @@ namespace internal // may only load from W after storing in X because the pointers might // point to the same memory w.load(W+i); - i += VectorizedArray::n_array_elements; return x * w; } @@ -749,15 +735,16 @@ namespace internal for (size_type i=start_chunk; i) { - const size_type regular_chunks = n_chunks/VectorizedArray::n_array_elements; + const unsigned int nvecs = VectorizedArray::n_array_elements; + const size_type regular_chunks = n_chunks/nvecs; for (size_type i=0; i r0 = op.do_vectorized(index); - VectorizedArray r1 = op.do_vectorized(index); - VectorizedArray r2 = op.do_vectorized(index); - VectorizedArray r3 = op.do_vectorized(index); - for (size_type j=1; j<8; ++j) + VectorizedArray r1 = op.do_vectorized(index+nvecs); + VectorizedArray r2 = op.do_vectorized(index+2*nvecs); + VectorizedArray r3 = op.do_vectorized(index+3*nvecs); + index += nvecs*4; + for (size_type j=1; j<8; ++j, index += nvecs*4) { r0 += op.do_vectorized(index); - r1 += op.do_vectorized(index); - r2 += op.do_vectorized(index); - r3 += op.do_vectorized(index); + r1 += op.do_vectorized(index+nvecs); + r2 += op.do_vectorized(index+2*nvecs); + r3 += op.do_vectorized(index+3*nvecs); } r0 += r1; r2 += r3; @@ -808,12 +797,12 @@ namespace internal { VectorizedArray r0 = VectorizedArray(), r1 = VectorizedArray(); - const size_type start_irreg = regular_chunks * VectorizedArray::n_array_elements; + const size_type start_irreg = regular_chunks * nvecs; for (size_type c=start_irreg; c::n_array_elements) + for (size_type j=0; j<32; j+=2*nvecs, index+=2*nvecs) { r0 += op.do_vectorized(index); - r1 += op.do_vectorized(index); + r1 += op.do_vectorized(index+nvecs); } r0 += r1; r0.store(&outer_results[start_irreg]); @@ -883,24 +872,24 @@ namespace internal switch (inner_chunks) { case 3: - r2 = op(index); + r2 = op(index++); for (size_type j=1; j<8; ++j) - r2 += op(index); + r2 += op(index++); // no break case 2: - r1 = op(index); + r1 = op(index++); for (size_type j=1; j<8; ++j) - r1 += op(index); + r1 += op(index++); r1 += r2; // no break case 1: - r2 = op(index); + r2 = op(index++); for (size_type j=1; j<8; ++j) - r2 += op(index); + r2 += op(index++); // no break default: for (size_type j=0; j