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,
57 if (format !=
nullptr)
58 std::printf (format, t);
60 std::printf (
" %5.2f",
double(t));
66 void print (
const std::complex<T> &t,
69 if (format !=
nullptr)
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::VectorImplementation::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::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::VectorImplementation::minimum_parallel_grain_size &&
173 Assert(partitioner.get() !=
nullptr,
175 "not set the TBB partitioner to a usable state."));
176 std::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)
210 void operator() (
const size_type begin,
const size_type end)
const 214 if (value == Number())
216 #ifdef DEAL_II_WITH_CXX17 217 if constexpr (std::is_trivial<Number>::value)
219 if (std::is_trivial<Number>::value)
222 std::memset(dst+begin, 0,
sizeof(Number)*(end-begin));
226 std::fill (dst+begin, dst+end, value);
233 template <
typename Number,
typename OtherNumber>
236 Vector_copy(
const OtherNumber *src, Number *dst)
245 void operator() (
const size_type begin,
const size_type end)
const 249 #if __GNUG__ && __GNUC__ < 5 250 if ( __has_trivial_copy(Number) && std::is_same<Number,OtherNumber>::value)
252 # ifdef DEAL_II_WITH_CXX17 253 if constexpr (std::is_trivially_copyable<Number>() && std::is_same<Number,OtherNumber>::value)
255 if (std::is_trivially_copyable<Number>() && std::is_same<Number,OtherNumber>::value)
258 std::memcpy(dst+begin, src+begin, (end-begin)*
sizeof(Number));
261 DEAL_II_OPENMP_SIMD_PRAGMA
262 for (size_type i=begin; i<end; ++i)
267 const OtherNumber *
const src;
271 template <
typename Number>
272 struct Vectorization_multiply_factor
274 Vectorization_multiply_factor(Number *val, Number factor)
280 void operator() (
const size_type begin,
const size_type end)
const 284 DEAL_II_OPENMP_SIMD_PRAGMA
285 for (size_type i=begin; i<end; ++i)
290 for (size_type i=begin; i<end; ++i)
299 template <
typename Number>
300 struct Vectorization_add_av
302 Vectorization_add_av(Number *val, Number *v_val, Number factor)
309 void operator() (
const size_type begin,
const size_type end)
const 313 DEAL_II_OPENMP_SIMD_PRAGMA
314 for (size_type i=begin; i<end; ++i)
315 val[i] += factor*v_val[i];
319 for (size_type i=begin; i<end; ++i)
320 val[i] += factor*v_val[i];
329 template <
typename Number>
330 struct Vectorization_sadd_xav
332 Vectorization_sadd_xav(Number *val, Number *v_val, Number a, Number x)
340 void operator() (
const size_type begin,
const size_type end)
const 344 DEAL_II_OPENMP_SIMD_PRAGMA
345 for (size_type i=begin; i<end; ++i)
346 val[i] = x*val[i] + a*v_val[i];
350 for (size_type i=begin; i<end; ++i)
351 val[i] = x*val[i] + a*v_val[i];
361 template <
typename Number>
362 struct Vectorization_subtract_v
364 Vectorization_subtract_v(Number *val, Number *v_val)
370 void operator() (
const size_type begin,
const size_type end)
const 374 DEAL_II_OPENMP_SIMD_PRAGMA
375 for (size_type i=begin; i<end; ++i)
380 for (size_type i=begin; i<end; ++i)
389 template <
typename Number>
390 struct Vectorization_add_factor
392 Vectorization_add_factor(Number *val, Number factor)
398 void operator() (
const size_type begin,
const size_type end)
const 402 DEAL_II_OPENMP_SIMD_PRAGMA
403 for (size_type i=begin; i<end; ++i)
408 for (size_type i=begin; i<end; ++i)
417 template <
typename Number>
418 struct Vectorization_add_v
420 Vectorization_add_v(Number *val, Number *v_val)
426 void operator() (
const size_type begin,
const size_type end)
const 430 DEAL_II_OPENMP_SIMD_PRAGMA
431 for (size_type i=begin; i<end; ++i)
436 for (size_type i=begin; i<end; ++i)
445 template <
typename Number>
446 struct Vectorization_add_avpbw
448 Vectorization_add_avpbw(Number *val, Number *v_val, Number *w_val, Number a, Number b)
457 void operator() (
const size_type begin,
const size_type end)
const 461 DEAL_II_OPENMP_SIMD_PRAGMA
462 for (size_type i=begin; i<end; ++i)
463 val[i] = val[i] + a*v_val[i] + b*w_val[i];
467 for (size_type i=begin; i<end; ++i)
468 val[i] = val[i] + a*v_val[i] + b*w_val[i];
479 template <
typename Number>
480 struct Vectorization_sadd_xv
482 Vectorization_sadd_xv(Number *val, Number *v_val, Number x)
489 void operator() (
const size_type begin,
const size_type end)
const 493 DEAL_II_OPENMP_SIMD_PRAGMA
494 for (size_type i=begin; i<end; ++i)
495 val[i] = x*val[i] + v_val[i];
499 for (size_type i=begin; i<end; ++i)
500 val[i] = x*val[i] + v_val[i];
509 template <
typename Number>
510 struct Vectorization_sadd_xavbw
512 Vectorization_sadd_xavbw(Number *val, Number *v_val, Number *w_val,
513 Number x, Number a, Number b)
523 void operator() (
const size_type begin,
const size_type end)
const 527 DEAL_II_OPENMP_SIMD_PRAGMA
528 for (size_type i=begin; i<end; ++i)
529 val[i] = x*val[i] + a*v_val[i] + b*w_val[i];
533 for (size_type i=begin; i<end; ++i)
534 val[i] = x*val[i] + a*v_val[i] + b*w_val[i];
546 template <
typename Number>
547 struct Vectorization_scale
549 Vectorization_scale(Number *val, Number *v_val)
555 void operator() (
const size_type begin,
const size_type end)
const 559 DEAL_II_OPENMP_SIMD_PRAGMA
560 for (size_type i=begin; i<end; ++i)
565 for (size_type i=begin; i<end; ++i)
574 template <
typename Number>
575 struct Vectorization_equ_au
577 Vectorization_equ_au(Number *val, Number *u_val, Number a)
584 void operator() (
const size_type begin,
const size_type end)
const 588 DEAL_II_OPENMP_SIMD_PRAGMA
589 for (size_type i=begin; i<end; ++i)
594 for (size_type i=begin; i<end; ++i)
604 template <
typename Number>
605 struct Vectorization_equ_aubv
607 Vectorization_equ_aubv(Number *val, Number *u_val, Number *v_val,
617 void operator() (
const size_type begin,
const size_type end)
const 621 DEAL_II_OPENMP_SIMD_PRAGMA
622 for (size_type i=begin; i<end; ++i)
623 val[i] = a*u_val[i] + b*v_val[i];
627 for (size_type i=begin; i<end; ++i)
628 val[i] = a*u_val[i] + b*v_val[i];
639 template <
typename Number>
640 struct Vectorization_equ_aubvcw
642 Vectorization_equ_aubvcw(Number *val, Number *u_val, Number *v_val,
643 Number *w_val, Number a, Number b, Number c)
654 void operator() (
const size_type begin,
const size_type end)
const 658 DEAL_II_OPENMP_SIMD_PRAGMA
659 for (size_type i=begin; i<end; ++i)
660 val[i] = a*u_val[i] + b*v_val[i] + c*w_val[i];
664 for (size_type i=begin; i<end; ++i)
665 val[i] = a*u_val[i] + b*v_val[i] + c*w_val[i];
678 template <
typename Number>
679 struct Vectorization_ratio
681 Vectorization_ratio(Number *val, Number *a_val, Number *b_val)
688 void operator() (
const size_type begin,
const size_type end)
const 692 DEAL_II_OPENMP_SIMD_PRAGMA
693 for (size_type i=begin; i<end; ++i)
694 val[i] = a_val[i]/b_val[i];
698 for (size_type i=begin; i<end; ++i)
699 val[i] = a_val[i]/b_val[i];
715 template <
typename Number,
typename Number2>
718 static const bool vectorizes = std::is_same<Number,Number2>::value &&
721 Dot(
const Number *X,
const Number2 *Y)
728 operator() (
const size_type i)
const 734 do_vectorized(
const size_type i)
const 748 "This operation is not correctly implemented for " 749 "complex-valued objects.");
757 template <
typename Number,
typename RealType>
762 Norm2(
const Number *X)
768 operator() (
const size_type i)
const 774 do_vectorized(
const size_type i)
const 784 template <
typename Number,
typename RealType>
789 Norm1(
const Number *X)
795 operator() (
const size_type i)
const 801 do_vectorized(
const size_type i)
const 811 template <
typename Number,
typename RealType>
816 NormP(
const Number *X, RealType p)
823 operator() (
const size_type i)
const 829 do_vectorized(
const size_type i)
const 833 return std::pow(std::abs(x),p);
840 template <
typename Number>
845 MeanValue(
const Number *X)
851 operator() (
const size_type i)
const 857 do_vectorized(
const size_type i)
const 867 template <
typename Number>
872 AddAndDot(Number *X,
const Number *V,
const Number *W, Number a)
881 operator() (
const size_type i)
const 888 do_vectorized(
const size_type i)
const 907 "This operation is not correctly implemented for " 908 "complex-valued objects.");
958 const unsigned int vector_accumulation_recursion_threshold = 128;
960 template <
typename Operation,
typename ResultType>
961 void accumulate_recursive (
const Operation &op,
962 const size_type first,
963 const size_type last,
967 if (vec_size <= vector_accumulation_recursion_threshold * 32)
973 ResultType outer_results [vector_accumulation_recursion_threshold];
977 outer_results[0] = ResultType();
984 const size_type remainder = vec_size % 32;
985 Assert (remainder == 0 || n_chunks < vector_accumulation_recursion_threshold,
992 accumulate_regular(op, n_chunks, index, outer_results,
993 std::integral_constant<bool, Operation::vectorizes>());
1005 const size_type inner_chunks = remainder / 8;
1007 const size_type remainder_inner = remainder % 8;
1008 ResultType r0 = ResultType(), r1 = ResultType(),
1010 switch (inner_chunks)
1014 for (size_type j=1; j<8; ++j)
1016 DEAL_II_FALLTHROUGH;
1019 for (size_type j=1; j<8; ++j)
1022 DEAL_II_FALLTHROUGH;
1025 for (size_type j=1; j<8; ++j)
1027 DEAL_II_FALLTHROUGH;
1029 for (size_type j=0; j<remainder_inner; ++j)
1033 if (n_chunks == vector_accumulation_recursion_threshold)
1034 outer_results[vector_accumulation_recursion_threshold-1] += r0;
1037 outer_results[n_chunks] = r0;
1048 while (n_chunks > 1)
1050 if (n_chunks % 2 == 1)
1051 outer_results[n_chunks++] = ResultType();
1052 for (size_type i=0; i<n_chunks; i+=2)
1053 outer_results[i/2] = outer_results[i] + outer_results[i+1];
1056 result = outer_results[0];
1064 (vec_size / (vector_accumulation_recursion_threshold * 32)) *
1065 vector_accumulation_recursion_threshold * 8;
1066 Assert (first+3*new_size < last,
1068 ResultType r0, r1, r2, r3;
1069 accumulate_recursive (op, first, first+new_size, r0);
1070 accumulate_recursive (op, first+new_size, first+2*new_size, r1);
1071 accumulate_recursive (op, first+2*new_size, first+3*new_size, r2);
1072 accumulate_recursive (op, first+3*new_size, last, r3);
1084 template <
typename Operation,
typename ResultType>
1086 accumulate_regular(
const Operation &op,
1087 size_type &n_chunks,
1089 ResultType (&outer_results)[vector_accumulation_recursion_threshold],
1090 std::integral_constant<bool, false>)
1094 for (size_type i=0; i<n_chunks; ++i)
1096 ResultType r0 = op(index);
1097 ResultType r1 = op(index+1);
1098 ResultType r2 = op(index+2);
1099 ResultType r3 = op(index+3);
1101 for (size_type j=1; j<8; ++j, index += 4)
1110 outer_results[i] = r0 + r2;
1121 template <
typename Operation,
typename Number>
1123 accumulate_regular(
const Operation &op,
1124 size_type &n_chunks,
1126 Number (&outer_results)[vector_accumulation_recursion_threshold],
1127 std::integral_constant<bool, true>)
1136 const size_type regular_chunks = n_chunks/nvecs;
1137 for (size_type i=0; i<regular_chunks; ++i)
1144 for (size_type j=1; j<8; ++j, index += nvecs*4)
1146 r0 += op.do_vectorized(index);
1147 r1 += op.do_vectorized(index+nvecs);
1148 r2 += op.do_vectorized(index+2*nvecs);
1149 r3 += op.do_vectorized(index+3*nvecs);
1170 const size_type start_irreg = regular_chunks * nvecs;
1171 for (size_type c=start_irreg; c<n_chunks; ++c)
1172 for (size_type j=0; j<32; j+=2*nvecs, index+=2*nvecs)
1174 r0 += op.do_vectorized(index);
1175 r1 += op.do_vectorized(index+nvecs);
1178 r0.
store(&outer_results[start_irreg]);
1187 #ifdef DEAL_II_WITH_THREADS 1216 template <
typename Operation,
typename ResultType>
1219 static const unsigned int threshold_array_allocate = 512;
1222 const size_type start,
1223 const size_type end)
1229 const size_type vec_size = end-start;
1231 const unsigned int gs = internal::VectorImplementation::minimum_parallel_grain_size;
1234 chunk_size = vec_size / n_chunks;
1240 if (chunk_size > 512)
1241 chunk_size = ((chunk_size + 511)/512)*512;
1242 n_chunks = (vec_size + chunk_size - 1) / chunk_size;
1246 if (n_chunks > threshold_array_allocate)
1250 large_array.resize(2*((n_chunks+1)/2));
1251 array_ptr = large_array.data();
1254 array_ptr = &small_array[0];
1261 void operator() (
const tbb::blocked_range<size_type> &range)
const 1263 for (size_type i = range.begin(); i < range.end(); ++i)
1264 accumulate_recursive(op, start+i*chunk_size, std::min(start+(i+1)*chunk_size, end),
1268 ResultType do_sum()
const 1270 while (n_chunks > 1)
1272 if (n_chunks % 2 == 1)
1273 array_ptr[n_chunks++] = ResultType();
1274 for (size_type i=0; i<n_chunks; i+=2)
1275 array_ptr[i/2] = array_ptr[i] + array_ptr[i+1];
1278 return array_ptr[0];
1281 const Operation &op;
1282 const size_type start;
1283 const size_type end;
1285 mutable unsigned int n_chunks;
1286 unsigned int chunk_size;
1287 ResultType small_array [threshold_array_allocate];
1288 std::vector<ResultType> large_array;
1291 mutable ResultType *array_ptr;
1301 template <
typename Operation,
typename ResultType>
1302 void parallel_reduce (
const Operation &op,
1303 const size_type start,
1304 const size_type end,
1306 std::shared_ptr<parallel::internal::TBBPartitioner> &partitioner)
1308 #ifdef DEAL_II_WITH_THREADS 1312 if (vec_size >= 4*internal::VectorImplementation::minimum_parallel_grain_size &&
1315 Assert(partitioner.get() !=
nullptr,
1317 "not set the TBB partitioner to a usable state."));
1318 std::shared_ptr<tbb::affinity_partitioner> tbb_partitioner =
1319 partitioner->acquire_one_partitioner();
1321 TBBReduceFunctor<Operation,ResultType> generic_functor(op, start, end);
1322 tbb::parallel_for (tbb::blocked_range<size_type> (0,
1323 generic_functor.n_chunks,
1327 partitioner->release_one_partitioner(tbb_partitioner);
1328 result = generic_functor.do_sum();
1331 accumulate_recursive(op,start,end,result);
1333 accumulate_recursive(op,start,end,result);
1340 DEAL_II_NAMESPACE_CLOSE
types::global_dof_index size_type
void store(Number *ptr) const
#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)
static ::ExceptionBase & ExcMessage(std::string arg1)
unsigned int global_dof_index
static real_type abs_square(const number &x)
#define Assert(cond, exc)
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)
void load(const Number *ptr)
static unsigned int n_threads()
static ::ExceptionBase & ExcInternalError()