17 #ifndef dealii__vector_operations_internal_h 18 #define dealii__vector_operations_internal_h 20 #include <deal.II/base/multithread_info.h> 21 #include <deal.II/base/parallel.h> 22 #include <deal.II/base/thread_management.h> 23 #include <deal.II/base/vectorization.h> 28 DEAL_II_NAMESPACE_OPEN
32 namespace VectorOperations
37 bool is_non_negative (
const T &t)
44 bool is_non_negative (
const std::complex<T> &)
47 ExcMessage (
"Complex numbers do not have an ordering."));
54 void print (
const T &t,
58 std::printf (format, t);
60 std::printf (
" %5.2f",
double(t));
66 void print (
const std::complex<T> &t,
70 std::printf (format, t.real(), t.imag());
72 std::printf (
" %5.2f+%5.2fi",
73 double(t.real()),
double(t.imag()));
80 template <
typename T,
typename U>
81 void copy (
const T *begin,
85 std::copy (begin, end, dest);
88 template <
typename T,
typename U>
89 void copy (
const std::complex<T> *begin,
90 const std::complex<T> *end,
91 std::complex<U> *dest)
93 std::copy (begin, end, dest);
96 template <
typename T,
typename U>
97 void copy (
const std::complex<T> *,
98 const std::complex<T> *,
102 "into a vector of reals/doubles"));
107 #ifdef DEAL_II_WITH_THREADS 116 template <
typename Functor>
120 const size_type start,
127 const size_type vec_size = end-start;
129 const unsigned int gs = internal::Vector::minimum_parallel_grain_size;
132 chunk_size = vec_size / n_chunks;
138 if (chunk_size > 512)
139 chunk_size = ((chunk_size + 511)/512)*512;
140 n_chunks = (vec_size + chunk_size - 1) / chunk_size;
145 void operator() (
const tbb::blocked_range<size_type> &range)
const 147 const size_type r_begin = start + range.begin()*chunk_size;
148 const size_type r_end = std::min(start + range.end()*chunk_size, end);
149 functor(r_begin, r_end);
153 const size_type start;
155 unsigned int n_chunks;
156 size_type chunk_size;
160 template <
typename Functor>
161 void parallel_for(Functor &functor,
164 std_cxx11::shared_ptr<parallel::internal::TBBPartitioner> &partitioner)
166 #ifdef DEAL_II_WITH_THREADS 167 size_type vec_size = end-start;
170 if (vec_size >= 4*internal::Vector::minimum_parallel_grain_size &&
173 Assert(partitioner.get() != NULL,
175 "not set the TBB partitioner to a usable state."));
176 std_cxx11::shared_ptr<tbb::affinity_partitioner> tbb_partitioner =
177 partitioner->acquire_one_partitioner();
180 tbb::parallel_for (tbb::blocked_range<size_type> (0,
181 generic_functor.n_chunks,
185 partitioner->release_one_partitioner(tbb_partitioner);
187 else if (vec_size > 0)
199 template <
typename Number>
202 Vector_set(Number value, Number *dst)
208 void operator() (
const size_type begin,
const size_type end)
const 210 if (value == Number())
211 std::memset (dst+begin,0,(end-begin)*
sizeof(Number));
213 std::fill (dst+begin, dst+end, value);
220 template <
typename Number,
typename OtherNumber>
223 Vector_copy(
const OtherNumber *src, Number *dst)
229 void operator() (
const size_type begin,
const size_type end)
const 232 std::memcpy(dst+begin, src+begin, (end-begin)*
sizeof(Number));
235 DEAL_II_OPENMP_SIMD_PRAGMA
236 for (size_type i=begin; i<end; ++i)
241 const OtherNumber *src;
245 template <
typename Number>
246 struct Vectorization_multiply_factor
248 Vectorization_multiply_factor(Number *val, Number factor)
254 void operator() (
const size_type begin,
const size_type end)
const 258 DEAL_II_OPENMP_SIMD_PRAGMA
259 for (size_type i=begin; i<end; ++i)
264 for (size_type i=begin; i<end; ++i)
273 template <
typename Number>
274 struct Vectorization_add_av
276 Vectorization_add_av(Number *val, Number *v_val, Number factor)
283 void operator() (
const size_type begin,
const size_type end)
const 287 DEAL_II_OPENMP_SIMD_PRAGMA
288 for (size_type i=begin; i<end; ++i)
289 val[i] += factor*v_val[i];
293 for (size_type i=begin; i<end; ++i)
294 val[i] += factor*v_val[i];
303 template <
typename Number>
304 struct Vectorization_sadd_xav
306 Vectorization_sadd_xav(Number *val, Number *v_val, Number a, Number x)
314 void operator() (
const size_type begin,
const size_type end)
const 318 DEAL_II_OPENMP_SIMD_PRAGMA
319 for (size_type i=begin; i<end; ++i)
320 val[i] = x*val[i] + a*v_val[i];
324 for (size_type i=begin; i<end; ++i)
325 val[i] = x*val[i] + a*v_val[i];
335 template <
typename Number>
336 struct Vectorization_subtract_v
338 Vectorization_subtract_v(Number *val, Number *v_val)
344 void operator() (
const size_type begin,
const size_type end)
const 348 DEAL_II_OPENMP_SIMD_PRAGMA
349 for (size_type i=begin; i<end; ++i)
354 for (size_type i=begin; i<end; ++i)
363 template <
typename Number>
364 struct Vectorization_add_factor
366 Vectorization_add_factor(Number *val, Number factor)
372 void operator() (
const size_type begin,
const size_type end)
const 376 DEAL_II_OPENMP_SIMD_PRAGMA
377 for (size_type i=begin; i<end; ++i)
382 for (size_type i=begin; i<end; ++i)
391 template <
typename Number>
392 struct Vectorization_add_v
394 Vectorization_add_v(Number *val, Number *v_val)
400 void operator() (
const size_type begin,
const size_type end)
const 404 DEAL_II_OPENMP_SIMD_PRAGMA
405 for (size_type i=begin; i<end; ++i)
410 for (size_type i=begin; i<end; ++i)
419 template <
typename Number>
420 struct Vectorization_add_avpbw
422 Vectorization_add_avpbw(Number *val, Number *v_val, Number *w_val, Number a, Number b)
431 void operator() (
const size_type begin,
const size_type end)
const 435 DEAL_II_OPENMP_SIMD_PRAGMA
436 for (size_type i=begin; i<end; ++i)
437 val[i] = val[i] + a*v_val[i] + b*w_val[i];
441 for (size_type i=begin; i<end; ++i)
442 val[i] = val[i] + a*v_val[i] + b*w_val[i];
453 template <
typename Number>
454 struct Vectorization_sadd_xv
456 Vectorization_sadd_xv(Number *val, Number *v_val, Number x)
463 void operator() (
const size_type begin,
const size_type end)
const 467 DEAL_II_OPENMP_SIMD_PRAGMA
468 for (size_type i=begin; i<end; ++i)
469 val[i] = x*val[i] + v_val[i];
473 for (size_type i=begin; i<end; ++i)
474 val[i] = x*val[i] + v_val[i];
483 template <
typename Number>
484 struct Vectorization_sadd_xavbw
486 Vectorization_sadd_xavbw(Number *val, Number *v_val, Number *w_val,
487 Number x, Number a, Number b)
497 void operator() (
const size_type begin,
const size_type end)
const 501 DEAL_II_OPENMP_SIMD_PRAGMA
502 for (size_type i=begin; i<end; ++i)
503 val[i] = x*val[i] + a*v_val[i] + b*w_val[i];
507 for (size_type i=begin; i<end; ++i)
508 val[i] = x*val[i] + a*v_val[i] + b*w_val[i];
520 template <
typename Number>
521 struct Vectorization_scale
523 Vectorization_scale(Number *val, Number *v_val)
529 void operator() (
const size_type begin,
const size_type end)
const 533 DEAL_II_OPENMP_SIMD_PRAGMA
534 for (size_type i=begin; i<end; ++i)
539 for (size_type i=begin; i<end; ++i)
548 template <
typename Number>
549 struct Vectorization_equ_au
551 Vectorization_equ_au(Number *val, Number *u_val, Number a)
558 void operator() (
const size_type begin,
const size_type end)
const 562 DEAL_II_OPENMP_SIMD_PRAGMA
563 for (size_type i=begin; i<end; ++i)
568 for (size_type i=begin; i<end; ++i)
578 template <
typename Number>
579 struct Vectorization_equ_aubv
581 Vectorization_equ_aubv(Number *val, Number *u_val, Number *v_val,
591 void operator() (
const size_type begin,
const size_type end)
const 595 DEAL_II_OPENMP_SIMD_PRAGMA
596 for (size_type i=begin; i<end; ++i)
597 val[i] = a*u_val[i] + b*v_val[i];
601 for (size_type i=begin; i<end; ++i)
602 val[i] = a*u_val[i] + b*v_val[i];
613 template <
typename Number>
614 struct Vectorization_equ_aubvcw
616 Vectorization_equ_aubvcw(Number *val, Number *u_val, Number *v_val,
617 Number *w_val, Number a, Number b, Number c)
628 void operator() (
const size_type begin,
const size_type end)
const 632 DEAL_II_OPENMP_SIMD_PRAGMA
633 for (size_type i=begin; i<end; ++i)
634 val[i] = a*u_val[i] + b*v_val[i] + c*w_val[i];
638 for (size_type i=begin; i<end; ++i)
639 val[i] = a*u_val[i] + b*v_val[i] + c*w_val[i];
652 template <
typename Number>
653 struct Vectorization_ratio
655 Vectorization_ratio(Number *val, Number *a_val, Number *b_val)
662 void operator() (
const size_type begin,
const size_type end)
const 666 DEAL_II_OPENMP_SIMD_PRAGMA
667 for (size_type i=begin; i<end; ++i)
668 val[i] = a_val[i]/b_val[i];
672 for (size_type i=begin; i<end; ++i)
673 val[i] = a_val[i]/b_val[i];
689 template <
typename Number,
typename Number2>
695 Dot(
const Number *X,
const Number2 *Y)
702 operator() (
const size_type i)
const 708 do_vectorized(
const size_type i)
const 720 template <
typename Number,
typename RealType>
725 Norm2(
const Number *X)
731 operator() (
const size_type i)
const 737 do_vectorized(
const size_type i)
const 747 template <
typename Number,
typename RealType>
752 Norm1(
const Number *X)
758 operator() (
const size_type i)
const 764 do_vectorized(
const size_type i)
const 774 template <
typename Number,
typename RealType>
779 NormP(
const Number *X, RealType p)
786 operator() (
const size_type i)
const 792 do_vectorized(
const size_type i)
const 796 return std::pow(std::abs(x),p);
803 template <
typename Number>
808 MeanValue(
const Number *X)
814 operator() (
const size_type i)
const 820 do_vectorized(
const size_type i)
const 830 template <
typename Number>
835 AddAndDot(Number *X,
const Number *V,
const Number *W, Number a)
844 operator() (
const size_type i)
const 851 do_vectorized(
const size_type i)
const 910 const unsigned int vector_accumulation_recursion_threshold = 128;
912 template <
typename Operation,
typename ResultType>
913 void accumulate_recursive (
const Operation &op,
914 const size_type first,
915 const size_type last,
919 if (vec_size <= vector_accumulation_recursion_threshold * 32)
925 ResultType outer_results [vector_accumulation_recursion_threshold];
929 outer_results[0] = ResultType();
936 const size_type remainder = vec_size % 32;
937 Assert (remainder == 0 || n_chunks < vector_accumulation_recursion_threshold,
944 accumulate_regular(op, n_chunks, index, outer_results,
957 const size_type inner_chunks = remainder / 8;
959 const size_type remainder_inner = remainder % 8;
960 ResultType r0 = ResultType(), r1 = ResultType(),
962 switch (inner_chunks)
966 for (size_type j=1; j<8; ++j)
971 for (size_type j=1; j<8; ++j)
977 for (size_type j=1; j<8; ++j)
981 for (size_type j=0; j<remainder_inner; ++j)
985 if (n_chunks == vector_accumulation_recursion_threshold)
986 outer_results[vector_accumulation_recursion_threshold-1] += r0;
989 outer_results[n_chunks] = r0;
1000 while (n_chunks > 1)
1002 if (n_chunks % 2 == 1)
1003 outer_results[n_chunks++] = ResultType();
1004 for (size_type i=0; i<n_chunks; i+=2)
1005 outer_results[i/2] = outer_results[i] + outer_results[i+1];
1008 result = outer_results[0];
1016 (vec_size / (vector_accumulation_recursion_threshold * 32)) *
1017 vector_accumulation_recursion_threshold * 8;
1018 Assert (first+3*new_size < last,
1020 ResultType r0, r1, r2, r3;
1021 accumulate_recursive (op, first, first+new_size, r0);
1022 accumulate_recursive (op, first+new_size, first+2*new_size, r1);
1023 accumulate_recursive (op, first+2*new_size, first+3*new_size, r2);
1024 accumulate_recursive (op, first+3*new_size, last, r3);
1036 template <
typename Operation,
typename ResultType>
1038 accumulate_regular(
const Operation &op,
1039 size_type &n_chunks,
1041 ResultType (&outer_results)[vector_accumulation_recursion_threshold],
1046 for (size_type i=0; i<n_chunks; ++i)
1048 ResultType r0 = op(index);
1049 ResultType r1 = op(index+1);
1050 ResultType r2 = op(index+2);
1051 ResultType r3 = op(index+3);
1053 for (size_type j=1; j<8; ++j, index += 4)
1062 outer_results[i] = r0 + r2;
1073 template <
typename Operation,
typename Number>
1075 accumulate_regular(
const Operation &op,
1076 size_type &n_chunks,
1078 Number (&outer_results)[vector_accumulation_recursion_threshold],
1088 const size_type regular_chunks = n_chunks/nvecs;
1089 for (size_type i=0; i<regular_chunks; ++i)
1096 for (size_type j=1; j<8; ++j, index += nvecs*4)
1098 r0 += op.do_vectorized(index);
1099 r1 += op.do_vectorized(index+nvecs);
1100 r2 += op.do_vectorized(index+2*nvecs);
1101 r3 += op.do_vectorized(index+3*nvecs);
1122 const size_type start_irreg = regular_chunks * nvecs;
1123 for (size_type c=start_irreg; c<n_chunks; ++c)
1124 for (size_type j=0; j<32; j+=2*nvecs, index+=2*nvecs)
1126 r0 += op.do_vectorized(index);
1127 r1 += op.do_vectorized(index+nvecs);
1130 r0.
store(&outer_results[start_irreg]);
1139 #ifdef DEAL_II_WITH_THREADS 1168 template <
typename Operation,
typename ResultType>
1171 static const unsigned int threshold_array_allocate = 512;
1174 const size_type start,
1175 const size_type end)
1181 const size_type vec_size = end-start;
1183 const unsigned int gs = internal::Vector::minimum_parallel_grain_size;
1186 chunk_size = vec_size / n_chunks;
1192 if (chunk_size > 512)
1193 chunk_size = ((chunk_size + 511)/512)*512;
1194 n_chunks = (vec_size + chunk_size - 1) / chunk_size;
1198 if (n_chunks > threshold_array_allocate)
1202 large_array.resize(2*((n_chunks+1)/2));
1203 array_ptr = &large_array[0];
1206 array_ptr = &small_array[0];
1213 void operator() (
const tbb::blocked_range<size_type> &range)
const 1215 for (size_type i = range.begin(); i < range.end(); ++i)
1216 accumulate_recursive(op, start+i*chunk_size, std::min(start+(i+1)*chunk_size, end),
1220 ResultType do_sum()
const 1222 while (n_chunks > 1)
1224 if (n_chunks % 2 == 1)
1225 array_ptr[n_chunks++] = ResultType();
1226 for (size_type i=0; i<n_chunks; i+=2)
1227 array_ptr[i/2] = array_ptr[i] + array_ptr[i+1];
1230 return array_ptr[0];
1233 const Operation &op;
1234 const size_type start;
1235 const size_type end;
1237 mutable unsigned int n_chunks;
1238 unsigned int chunk_size;
1239 ResultType small_array [threshold_array_allocate];
1240 std::vector<ResultType> large_array;
1243 mutable ResultType *array_ptr;
1253 template <
typename Operation,
typename ResultType>
1254 void parallel_reduce (
const Operation &op,
1255 const size_type start,
1256 const size_type end,
1258 std_cxx11::shared_ptr<parallel::internal::TBBPartitioner> &partitioner)
1260 #ifdef DEAL_II_WITH_THREADS 1264 if (vec_size >= 4*internal::Vector::minimum_parallel_grain_size &&
1267 Assert(partitioner.get() != NULL,
1269 "not set the TBB partitioner to a usable state."));
1270 std_cxx11::shared_ptr<tbb::affinity_partitioner> tbb_partitioner =
1271 partitioner->acquire_one_partitioner();
1273 TBBReduceFunctor<Operation,ResultType> generic_functor(op, start, end);
1274 tbb::parallel_for (tbb::blocked_range<size_type> (0,
1275 generic_functor.n_chunks,
1279 partitioner->release_one_partitioner(tbb_partitioner);
1280 result = generic_functor.do_sum();
1283 accumulate_recursive(op,start,end,result);
1285 accumulate_recursive(op,start,end,result);
1292 DEAL_II_NAMESPACE_CLOSE
types::global_dof_index size_type
#define AssertDimension(dim1, dim2)
void operator()(const tbb::blocked_range< size_type > &range) const
#define AssertIndexRange(index, range)
static real_type abs(const number &x)
DEAL_II_ALWAYS_INLINE void load(const Number *ptr)
static ::ExceptionBase & ExcMessage(std::string arg1)
unsigned int global_dof_index
static real_type abs_square(const number &x)
#define Assert(cond, exc)
DEAL_II_ALWAYS_INLINE void store(Number *ptr) const
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static unsigned int n_threads()
static ::ExceptionBase & ExcInternalError()