Reference documentation for deal.II version 9.3.3
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
vector_operations_internal.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2016 - 2020 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
17#ifndef dealii_vector_operations_internal_h
18#define dealii_vector_operations_internal_h
19
20#include <deal.II/base/config.h>
21
26#include <deal.II/base/types.h>
28
30#include <deal.II/lac/cuda_kernels.templates.h>
32
33#include <cstdio>
34#include <cstring>
35
37
38namespace internal
39{
40 namespace VectorOperations
41 {
43
44 template <typename T>
45 bool
47 {
48 return t >= 0;
49 }
50
51
52 template <typename T>
53 bool
54 is_non_negative(const std::complex<T> &)
55 {
56 Assert(false, ExcMessage("Complex numbers do not have an ordering."));
57
58 return false;
59 }
60
61
62 // call std::copy, except for in
63 // the case where we want to copy
64 // from std::complex to a
65 // non-complex type
66 template <typename T, typename U>
67 void
68 copy(const T *begin, const T *end, U *dest)
69 {
70 std::copy(begin, end, dest);
71 }
72
73 template <typename T, typename U>
74 void
75 copy(const std::complex<T> *begin,
76 const std::complex<T> *end,
77 std::complex<U> * dest)
78 {
79 std::copy(begin, end, dest);
80 }
81
82 template <typename T, typename U>
83 void
84 copy(const std::complex<T> *, const std::complex<T> *, U *)
85 {
86 Assert(false,
87 ExcMessage("Can't convert a vector of complex numbers "
88 "into a vector of reals/doubles"));
89 }
90
91
92
93#ifdef DEAL_II_WITH_TBB
102 template <typename Functor>
104 {
106 const size_type start,
107 const size_type end)
109 , start(start)
110 , end(end)
111 {
112 const size_type vec_size = end - start;
113 // set chunk size for sub-tasks
114 const unsigned int gs =
116 n_chunks =
118 vec_size / gs);
119 chunk_size = vec_size / n_chunks;
120
121 // round to next multiple of 512 (or minimum grain size if that happens
122 // to be smaller). this is advantageous because our accumulation
123 // algorithms favor lengths of a power of 2 due to pairwise summation ->
124 // at most one 'oddly' sized chunk
125 if (chunk_size > 512)
126 chunk_size = ((chunk_size + 511) / 512) * 512;
127 n_chunks = (vec_size + chunk_size - 1) / chunk_size;
128 AssertIndexRange((n_chunks - 1) * chunk_size, vec_size);
129 AssertIndexRange(vec_size, n_chunks * chunk_size + 1);
130 }
131
132 void
133 operator()(const tbb::blocked_range<size_type> &range) const
134 {
135 const size_type r_begin = start + range.begin() * chunk_size;
136 const size_type r_end = std::min(start + range.end() * chunk_size, end);
137 functor(r_begin, r_end);
138 }
139
140 Functor & functor;
143 unsigned int n_chunks;
145 };
146#endif
147
148 template <typename Functor>
149 void
151 Functor & functor,
152 const size_type start,
153 const size_type end,
154 const std::shared_ptr<::parallel::internal::TBBPartitioner>
155 &partitioner)
156 {
157#ifdef DEAL_II_WITH_TBB
158 const size_type vec_size = end - start;
159 // only go to the parallel function in case there are at least 4 parallel
160 // items, otherwise the overhead is too large
161 if (vec_size >=
164 {
165 Assert(partitioner.get() != nullptr,
167 "Unexpected initialization of Vector that does "
168 "not set the TBB partitioner to a usable state."));
169 std::shared_ptr<tbb::affinity_partitioner> tbb_partitioner =
170 partitioner->acquire_one_partitioner();
171
172 TBBForFunctor<Functor> generic_functor(functor, start, end);
173 // We use a minimum grain size of 1 here since the grains at this
174 // stage of dividing the work refer to the number of vector chunks
175 // that are processed by (possibly different) threads in the
176 // parallelized for loop (i.e., they do not refer to individual
177 // vector entries). The number of chunks here is calculated inside
178 // TBBForFunctor. See also GitHub issue #2496 for further discussion
179 // of this strategy.
181 static_cast<size_type>(0),
182 static_cast<size_type>(generic_functor.n_chunks),
183 generic_functor,
184 1,
185 tbb_partitioner);
186 partitioner->release_one_partitioner(tbb_partitioner);
187 }
188 else if (vec_size > 0)
189 functor(start, end);
190#else
191 functor(start, end);
192 (void)partitioner;
193#endif
194 }
195
196
197 // Define the functors necessary to use SIMD with TBB. we also include the
198 // simple copy and set operations
199
200 template <typename Number>
202 {
203 Vector_set(const Number value, Number *const dst)
204 : value(value)
205 , dst(dst)
206 {
207 Assert(dst != nullptr, ExcInternalError());
208 }
209
210 void
212 {
214
215 if (value == Number())
216 {
217#ifdef DEAL_II_HAVE_CXX17
218 if constexpr (std::is_trivial<Number>::value)
219#else
220 if (std::is_trivial<Number>::value)
221#endif
222 {
223 std::memset(dst + begin, 0, sizeof(Number) * (end - begin));
224 return;
225 }
226 }
227 std::fill(dst + begin, dst + end, value);
228 }
229
230 const Number value;
231 Number *const dst;
232 };
233
234 template <typename Number, typename OtherNumber>
236 {
237 Vector_copy(const OtherNumber *const src, Number *const dst)
238 : src(src)
239 , dst(dst)
240 {
241 Assert(src != nullptr, ExcInternalError());
242 Assert(dst != nullptr, ExcInternalError());
243 }
244
245 void
247 {
249
250#if __GNUG__ && __GNUC__ < 5
251 if (__has_trivial_copy(Number) &&
252 std::is_same<Number, OtherNumber>::value)
253#else
254# ifdef DEAL_II_HAVE_CXX17
255 if constexpr (std::is_trivially_copyable<Number>() &&
256 std::is_same<Number, OtherNumber>::value)
257# else
258 if (std::is_trivially_copyable<Number>() &&
259 std::is_same<Number, OtherNumber>::value)
260# endif
261#endif
262 std::memcpy(dst + begin, src + begin, (end - begin) * sizeof(Number));
263 else
264 {
266 for (size_type i = begin; i < end; ++i)
267 dst[i] = src[i];
268 }
269 }
270
271 const OtherNumber *const src;
272 Number *const dst;
273 };
274
275 template <typename Number>
277 {
278 Vectorization_multiply_factor(Number *const val, const Number factor)
279 : val(val)
280 , factor(factor)
281 {}
282
283 void
285 {
287 {
289 for (size_type i = begin; i < end; ++i)
290 val[i] *= factor;
291 }
292 else
293 {
294 for (size_type i = begin; i < end; ++i)
295 val[i] *= factor;
296 }
297 }
298
299 Number *const val;
300 const Number factor;
301 };
302
303 template <typename Number>
305 {
307 const Number *const v_val,
308 const Number factor)
309 : val(val)
310 , v_val(v_val)
311 , factor(factor)
312 {}
313
314 void
316 {
318 {
320 for (size_type i = begin; i < end; ++i)
321 val[i] += factor * v_val[i];
322 }
323 else
324 {
325 for (size_type i = begin; i < end; ++i)
326 val[i] += factor * v_val[i];
327 }
328 }
329
330 Number *const val;
331 const Number *const v_val;
332 const Number factor;
333 };
334
335 template <typename Number>
337 {
339 const Number *const v_val,
340 const Number a,
341 const Number x)
342 : val(val)
343 , v_val(v_val)
344 , a(a)
345 , x(x)
346 {}
347
348 void
350 {
352 {
354 for (size_type i = begin; i < end; ++i)
355 val[i] = x * val[i] + a * v_val[i];
356 }
357 else
358 {
359 for (size_type i = begin; i < end; ++i)
360 val[i] = x * val[i] + a * v_val[i];
361 }
362 }
363
364 Number *const val;
365 const Number *const v_val;
366 const Number a;
367 const Number x;
368 };
369
370 template <typename Number>
372 {
373 Vectorization_subtract_v(Number *val, const Number *const v_val)
374 : val(val)
375 , v_val(v_val)
376 {}
377
378 void
380 {
382 {
384 for (size_type i = begin; i < end; ++i)
385 val[i] -= v_val[i];
386 }
387 else
388 {
389 for (size_type i = begin; i < end; ++i)
390 val[i] -= v_val[i];
391 }
392 }
393
394 Number *const val;
395 const Number *const v_val;
396 };
397
398 template <typename Number>
400 {
401 Vectorization_add_factor(Number *const val, const Number factor)
402 : val(val)
403 , factor(factor)
404 {}
405
406 void
408 {
410 {
412 for (size_type i = begin; i < end; ++i)
413 val[i] += factor;
414 }
415 else
416 {
417 for (size_type i = begin; i < end; ++i)
418 val[i] += factor;
419 }
420 }
421
422 Number *const val;
423 const Number factor;
424 };
425
426 template <typename Number>
428 {
429 Vectorization_add_v(Number *const val, const Number *const v_val)
430 : val(val)
431 , v_val(v_val)
432 {}
433
434 void
436 {
438 {
440 for (size_type i = begin; i < end; ++i)
441 val[i] += v_val[i];
442 }
443 else
444 {
445 for (size_type i = begin; i < end; ++i)
446 val[i] += v_val[i];
447 }
448 }
449
450 Number *const val;
451 const Number *const v_val;
452 };
453
454 template <typename Number>
456 {
458 const Number *const v_val,
459 const Number *const w_val,
460 const Number a,
461 const Number b)
462 : val(val)
463 , v_val(v_val)
464 , w_val(w_val)
465 , a(a)
466 , b(b)
467 {}
468
469 void
471 {
473 {
475 for (size_type i = begin; i < end; ++i)
476 val[i] = val[i] + a * v_val[i] + b * w_val[i];
477 }
478 else
479 {
480 for (size_type i = begin; i < end; ++i)
481 val[i] = val[i] + a * v_val[i] + b * w_val[i];
482 }
483 }
484
485 Number *const val;
486 const Number *const v_val;
487 const Number *const w_val;
488 const Number a;
489 const Number b;
490 };
491
492 template <typename Number>
494 {
496 const Number *const v_val,
497 const Number x)
498 : val(val)
499 , v_val(v_val)
500 , x(x)
501 {}
502
503 void
505 {
507 {
509 for (size_type i = begin; i < end; ++i)
510 val[i] = x * val[i] + v_val[i];
511 }
512 else
513 {
514 for (size_type i = begin; i < end; ++i)
515 val[i] = x * val[i] + v_val[i];
516 }
517 }
518
519 Number *const val;
520 const Number *const v_val;
521 const Number x;
522 };
523
524 template <typename Number>
526 {
528 const Number *v_val,
529 const Number *w_val,
530 Number x,
531 Number a,
532 Number b)
533 : val(val)
534 , v_val(v_val)
535 , w_val(w_val)
536 , x(x)
537 , a(a)
538 , b(b)
539 {}
540
541 void
543 {
545 {
547 for (size_type i = begin; i < end; ++i)
548 val[i] = x * val[i] + a * v_val[i] + b * w_val[i];
549 }
550 else
551 {
552 for (size_type i = begin; i < end; ++i)
553 val[i] = x * val[i] + a * v_val[i] + b * w_val[i];
554 }
555 }
556
557 Number *const val;
558 const Number *const v_val;
559 const Number *const w_val;
560 const Number x;
561 const Number a;
562 const Number b;
563 };
564
565 template <typename Number>
567 {
568 Vectorization_scale(Number *const val, const Number *const v_val)
569 : val(val)
570 , v_val(v_val)
571 {}
572
573 void
575 {
577 {
579 for (size_type i = begin; i < end; ++i)
580 val[i] *= v_val[i];
581 }
582 else
583 {
584 for (size_type i = begin; i < end; ++i)
585 val[i] *= v_val[i];
586 }
587 }
588
589 Number *const val;
590 const Number *const v_val;
591 };
592
593 template <typename Number>
595 {
597 const Number *const u_val,
598 const Number a)
599 : val(val)
600 , u_val(u_val)
601 , a(a)
602 {}
603
604 void
606 {
608 {
610 for (size_type i = begin; i < end; ++i)
611 val[i] = a * u_val[i];
612 }
613 else
614 {
615 for (size_type i = begin; i < end; ++i)
616 val[i] = a * u_val[i];
617 }
618 }
619
620 Number *const val;
621 const Number *const u_val;
622 const Number a;
623 };
624
625 template <typename Number>
627 {
629 const Number *const u_val,
630 const Number *const v_val,
631 const Number a,
632 const Number b)
633 : val(val)
634 , u_val(u_val)
635 , v_val(v_val)
636 , a(a)
637 , b(b)
638 {}
639
640 void
642 {
644 {
646 for (size_type i = begin; i < end; ++i)
647 val[i] = a * u_val[i] + b * v_val[i];
648 }
649 else
650 {
651 for (size_type i = begin; i < end; ++i)
652 val[i] = a * u_val[i] + b * v_val[i];
653 }
654 }
655
656 Number *const val;
657 const Number *const u_val;
658 const Number *const v_val;
659 const Number a;
660 const Number b;
661 };
662
663 template <typename Number>
665 {
667 const Number *u_val,
668 const Number *v_val,
669 const Number *w_val,
670 const Number a,
671 const Number b,
672 const Number c)
673 : val(val)
674 , u_val(u_val)
675 , v_val(v_val)
676 , w_val(w_val)
677 , a(a)
678 , b(b)
679 , c(c)
680 {}
681
682 void
684 {
686 {
688 for (size_type i = begin; i < end; ++i)
689 val[i] = a * u_val[i] + b * v_val[i] + c * w_val[i];
690 }
691 else
692 {
693 for (size_type i = begin; i < end; ++i)
694 val[i] = a * u_val[i] + b * v_val[i] + c * w_val[i];
695 }
696 }
697
698 Number *const val;
699 const Number *const u_val;
700 const Number *const v_val;
701 const Number *const w_val;
702 const Number a;
703 const Number b;
704 const Number c;
705 };
706
707 template <typename Number>
709 {
710 Vectorization_ratio(Number *val, const Number *a_val, const Number *b_val)
711 : val(val)
712 , a_val(a_val)
713 , b_val(b_val)
714 {}
715
716 void
718 {
720 {
722 for (size_type i = begin; i < end; ++i)
723 val[i] = a_val[i] / b_val[i];
724 }
725 else
726 {
727 for (size_type i = begin; i < end; ++i)
728 val[i] = a_val[i] / b_val[i];
729 }
730 }
731
732 Number *const val;
733 const Number *const a_val;
734 const Number *const b_val;
735 };
736
737
738
739 // All sums over all the vector entries (l2-norm, inner product, etc.) are
740 // performed with the same code, using a templated operation defined
741 // here. There are always two versions defined, a standard one that covers
742 // most cases and a vectorized one which is only for equal types and float
743 // and double.
744 template <typename Number, typename Number2>
745 struct Dot
746 {
747 static constexpr bool vectorizes = std::is_same<Number, Number2>::value &&
749
750 Dot(const Number *const X, const Number2 *const Y)
751 : X(X)
752 , Y(Y)
753 {}
754
755 Number
756 operator()(const size_type i) const
757 {
758 return X[i] * Number(numbers::NumberTraits<Number2>::conjugate(Y[i]));
759 }
760
763 {
765 x.load(X + i);
766 y.load(Y + i);
767
768 // the following operation in VectorizedArray does an element-wise
769 // scalar product without taking into account complex values and
770 // the need to take the complex-conjugate of one argument. this
771 // may be a bug, but because all VectorizedArray classes only
772 // work on real scalars, it doesn't really matter very much.
773 // in any case, assert that we really don't get here for
774 // complex-valued objects
775 static_assert(numbers::NumberTraits<Number>::is_complex == false,
776 "This operation is not correctly implemented for "
777 "complex-valued objects.");
778 return x * y;
779 }
780
781 const Number *const X;
782 const Number2 *const Y;
783 };
784
785 template <typename Number, typename RealType>
786 struct Norm2
787 {
788 static const bool vectorizes = VectorizedArray<Number>::size() > 1;
789
790 Norm2(const Number *const X)
791 : X(X)
792 {}
793
794 RealType
795 operator()(const size_type i) const
796 {
798 }
799
802 {
804 x.load(X + i);
805 return x * x;
806 }
807
808 const Number *const X;
809 };
810
811 template <typename Number, typename RealType>
812 struct Norm1
813 {
814 static const bool vectorizes = VectorizedArray<Number>::size() > 1;
815
816 Norm1(const Number *X)
817 : X(X)
818 {}
819
820 RealType
821 operator()(const size_type i) const
822 {
824 }
825
828 {
830 x.load(X + i);
831 return std::abs(x);
832 }
833
834 const Number *X;
835 };
836
837 template <typename Number, typename RealType>
838 struct NormP
839 {
840 static const bool vectorizes = VectorizedArray<Number>::size() > 1;
841
842 NormP(const Number *X, RealType p)
843 : X(X)
844 , p(p)
845 {}
846
847 RealType
848 operator()(const size_type i) const
849 {
851 }
852
855 {
857 x.load(X + i);
858 return std::pow(std::abs(x), p);
859 }
860
861 const Number * X;
862 const RealType p;
863 };
864
865 template <typename Number>
867 {
868 static const bool vectorizes = VectorizedArray<Number>::size() > 1;
869
870 MeanValue(const Number *X)
871 : X(X)
872 {}
873
874 Number
875 operator()(const size_type i) const
876 {
877 return X[i];
878 }
879
882 {
884 x.load(X + i);
885 return x;
886 }
887
888 const Number *X;
889 };
890
891 template <typename Number>
893 {
894 static const bool vectorizes = VectorizedArray<Number>::size() > 1;
895
896 AddAndDot(Number *const X,
897 const Number *const V,
898 const Number *const W,
899 const Number a)
900 : X(X)
901 , V(V)
902 , W(W)
903 , a(a)
904 {}
905
906 Number
907 operator()(const size_type i) const
908 {
909 X[i] += a * V[i];
910 return X[i] * Number(numbers::NumberTraits<Number>::conjugate(W[i]));
911 }
912
915 {
917 x.load(X + i);
918 v.load(V + i);
919 x += a * v;
920 x.store(X + i);
921 // may only load from W after storing in X because the pointers might
922 // point to the same memory
923 w.load(W + i);
924
925 // the following operation in VectorizedArray does an element-wise
926 // scalar product without taking into account complex values and
927 // the need to take the complex-conjugate of one argument. this
928 // may be a bug, but because all VectorizedArray classes only
929 // work on real scalars, it doesn't really matter very much.
930 // in any case, assert that we really don't get here for
931 // complex-valued objects
932 static_assert(numbers::NumberTraits<Number>::is_complex == false,
933 "This operation is not correctly implemented for "
934 "complex-valued objects.");
935 return x * w;
936 }
937
938 Number *const X;
939 const Number *const V;
940 const Number *const W;
941 const Number a;
942 };
943
944
945
946 // this is the main working loop for all vector sums using the templated
947 // operation above. it accumulates the sums using a block-wise summation
948 // algorithm with post-update. this blocked algorithm has been proposed in
949 // a similar form by Castaldo, Whaley and Chronopoulos (SIAM
950 // J. Sci. Comput. 31, 1156-1174, 2008) and we use the smallest possible
951 // block size, 2. Sometimes it is referred to as pairwise summation. The
952 // worst case error made by this algorithm is on the order O(eps *
953 // log2(vec_size)), whereas a naive summation is O(eps * vec_size). Even
954 // though the Kahan summation is even more accurate with an error O(eps)
955 // by carrying along remainders not captured by the main sum, that involves
956 // additional costs which are not worthwhile. See the Wikipedia article on
957 // the Kahan summation algorithm.
958
959 // The algorithm implemented here has the additional benefit that it is
960 // easily parallelized without changing the order of how the elements are
961 // added (floating point addition is not associative). For the same vector
962 // size and minimum_parallel_grainsize, the blocks are always the
963 // same and added pairwise.
964
965 // The depth of recursion is controlled by the 'magic' parameter
966 // vector_accumulation_recursion_threshold: If the length is below
967 // vector_accumulation_recursion_threshold * 32 (32 is the part of code we
968 // unroll), a straight loop instead of recursion will be used. At the
969 // innermost level, eight values are added consecutively in order to better
970 // balance multiplications and additions.
971
972 // Loops are unrolled as follows: the range [first,last) is broken into
973 // @p n_chunks each of size 32 plus the @p remainder.
974 // accumulate_regular() does the work on 32*n_chunks elements employing SIMD
975 // if possible and stores the result of the operation for each chunk in @p outer_results.
976
977 // The code returns the result as the last argument in order to make
978 // spawning tasks simpler and use automatic template deduction.
979
980
987
988 template <typename Operation, typename ResultType>
989 void
990 accumulate_recursive(const Operation &op,
991 const size_type first,
992 const size_type last,
993 ResultType & result)
994 {
995 const size_type vec_size = last - first;
996 if (vec_size <= vector_accumulation_recursion_threshold * 32)
997 {
998 // the vector is short enough so we perform the summation. first
999 // work on the regular part. The innermost 32 values are expanded in
1000 // order to obtain known loop bounds for most of the work.
1001 size_type index = first;
1002 ResultType outer_results[vector_accumulation_recursion_threshold];
1003
1004 // set the zeroth element to zero to correctly handle the case where
1005 // vec_size == 0
1006 outer_results[0] = ResultType();
1007
1008 // the variable serves two purposes: (i) number of chunks (each 32
1009 // indices) for the given size; all results are stored in
1010 // outer_results[0,n_chunks) (ii) in the SIMD case n_chunks is also a
1011 // next free index in outer_results[] to which we can write after
1012 // accumulate_regular() is executed.
1013 size_type n_chunks = vec_size / 32;
1014 const size_type remainder = vec_size % 32;
1015 Assert(remainder == 0 ||
1018
1019 // Select between the regular version and vectorized version based
1020 // on the number types we are given. To choose the vectorized
1021 // version often enough, we need to have all tasks but the last one
1022 // to be divisible by the vectorization length
1024 op,
1025 n_chunks,
1026 index,
1027 outer_results,
1028 std::integral_constant<bool, Operation::vectorizes>());
1029
1030 // now work on the remainder, i.e., the last up to 32 values. Use
1031 // switch statement with fall-through to work on these values.
1032 if (remainder > 0)
1033 {
1034 // if we got here, it means that (vec_size <=
1035 // vector_accumulation_recursion_threshold * 32), which is to say
1036 // that the domain can be split into n_chunks <=
1037 // vector_accumulation_recursion_threshold:
1038 AssertIndexRange(n_chunks,
1040 // split the remainder into chunks of 8, there could be up to 3
1041 // such chunks since remainder < 32.
1042 // Work on those chunks without any SIMD, that is we call
1043 // op(index).
1044 const size_type inner_chunks = remainder / 8;
1045 Assert(inner_chunks <= 3, ExcInternalError());
1046 const size_type remainder_inner = remainder % 8;
1047 ResultType r0 = ResultType(), r1 = ResultType(),
1048 r2 = ResultType();
1049 switch (inner_chunks)
1050 {
1051 case 3:
1052 r2 = op(index++);
1053 for (size_type j = 1; j < 8; ++j)
1054 r2 += op(index++);
1056 case 2:
1057 r1 = op(index++);
1058 for (size_type j = 1; j < 8; ++j)
1059 r1 += op(index++);
1060 r1 += r2;
1062 case 1:
1063 r2 = op(index++);
1064 for (size_type j = 1; j < 8; ++j)
1065 r2 += op(index++);
1067 default:
1068 for (size_type j = 0; j < remainder_inner; ++j)
1069 r0 += op(index++);
1070 r0 += r2;
1071 r0 += r1;
1074 1] += r0;
1075 else
1076 {
1077 outer_results[n_chunks] = r0;
1078 n_chunks++;
1079 }
1080 break;
1081 }
1082 }
1083 // make sure we worked through all indices
1084 AssertDimension(index, last);
1085
1086 // now sum the results from the chunks stored in
1087 // outer_results[0,n_chunks) recursively
1088 while (n_chunks > 1)
1089 {
1090 if (n_chunks % 2 == 1)
1091 outer_results[n_chunks++] = ResultType();
1092 for (size_type i = 0; i < n_chunks; i += 2)
1093 outer_results[i / 2] = outer_results[i] + outer_results[i + 1];
1094 n_chunks /= 2;
1095 }
1096 result = outer_results[0];
1097 }
1098 else
1099 {
1100 // split vector into four pieces and work on the pieces
1101 // recursively. Make pieces (except last) divisible by one fourth the
1102 // recursion threshold.
1103 const size_type new_size =
1104 (vec_size / (vector_accumulation_recursion_threshold * 32)) *
1106 Assert(first + 3 * new_size < last, ExcInternalError());
1107 ResultType r0, r1, r2, r3;
1108 accumulate_recursive(op, first, first + new_size, r0);
1109 accumulate_recursive(op, first + new_size, first + 2 * new_size, r1);
1111 first + 2 * new_size,
1112 first + 3 * new_size,
1113 r2);
1114 accumulate_recursive(op, first + 3 * new_size, last, r3);
1115 r0 += r1;
1116 r2 += r3;
1117 result = r0 + r2;
1118 }
1119 }
1120
1121
1122 // this is the inner working routine for the accumulation loops
1123 // below. This is the standard case where the loop bounds are known. We
1124 // pulled this function out of the regular accumulate routine because we
1125 // might do this thing vectorized (see specialized function below)
1126 template <typename Operation, typename ResultType>
1127 void
1129 const Operation &op,
1130 const size_type &n_chunks,
1131 size_type & index,
1132 ResultType (&outer_results)[vector_accumulation_recursion_threshold],
1133 std::integral_constant<bool, false>)
1134 {
1135 // note that each chunk is chosen to have a width of 32, thereby the index
1136 // is incremented by 4*8 for each @p i.
1137 for (size_type i = 0; i < n_chunks; ++i)
1138 {
1139 ResultType r0 = op(index);
1140 ResultType r1 = op(index + 1);
1141 ResultType r2 = op(index + 2);
1142 ResultType r3 = op(index + 3);
1143 index += 4;
1144 for (size_type j = 1; j < 8; ++j, index += 4)
1145 {
1146 r0 += op(index);
1147 r1 += op(index + 1);
1148 r2 += op(index + 2);
1149 r3 += op(index + 3);
1150 }
1151 r0 += r1;
1152 r2 += r3;
1153 outer_results[i] = r0 + r2;
1154 }
1155 }
1156
1157
1158
1159 // this is the inner working routine for the accumulation loops
1160 // below. This is the specialized case where the loop bounds are known and
1161 // where we can vectorize. In that case, we request the 'do_vectorized'
1162 // routine of the operation instead of the regular one which does several
1163 // operations at once.
1164 template <typename Operation, typename Number>
1165 void
1167 const Operation &op,
1168 size_type & n_chunks,
1169 size_type & index,
1170 Number (&outer_results)[vector_accumulation_recursion_threshold],
1171 std::integral_constant<bool, true>)
1172 {
1173 // we start from @p index and workout @p n_chunks each of size 32.
1174 // in order employ SIMD and work on @p nvecs at a time, we split this
1175 // loop yet again:
1176 // First we work on (n_chunks/nvecs) chunks, where each chunk processes
1177 // nvecs*(4*8) elements.
1178
1179 constexpr unsigned int nvecs = VectorizedArray<Number>::size();
1180 const size_type regular_chunks = n_chunks / nvecs;
1181 for (size_type i = 0; i < regular_chunks; ++i)
1182 {
1183 VectorizedArray<Number> r0 = op.do_vectorized(index);
1184 VectorizedArray<Number> r1 = op.do_vectorized(index + nvecs);
1185 VectorizedArray<Number> r2 = op.do_vectorized(index + 2 * nvecs);
1186 VectorizedArray<Number> r3 = op.do_vectorized(index + 3 * nvecs);
1187 index += nvecs * 4;
1188 for (size_type j = 1; j < 8; ++j, index += nvecs * 4)
1189 {
1190 r0 += op.do_vectorized(index);
1191 r1 += op.do_vectorized(index + nvecs);
1192 r2 += op.do_vectorized(index + 2 * nvecs);
1193 r3 += op.do_vectorized(index + 3 * nvecs);
1194 }
1195 r0 += r1;
1196 r2 += r3;
1197 r0 += r2;
1198 r0.store(&outer_results[i * nvecs]);
1199 }
1200
1201 // If we are treating a case where the vector length is not divisible by
1202 // the vectorization length, need a cleanup loop
1203 // The remaining chunks are processed one by one starting from
1204 // regular_chunks * nvecs; We do as much as possible with 2 SIMD
1205 // operations within each chunk. Here we assume that nvecs < 32/2 = 16 as
1206 // well as 16%nvecs==0.
1207 static_assert(
1210 "VectorizedArray::size() must be a power of 2 and not more than 16");
1211 Assert(16 % nvecs == 0, ExcInternalError());
1212 if (n_chunks % nvecs != 0)
1213 {
1216 const size_type start_irreg = regular_chunks * nvecs;
1217 for (size_type c = start_irreg; c < n_chunks; ++c)
1218 for (size_type j = 0; j < 32; j += 2 * nvecs, index += 2 * nvecs)
1219 {
1220 r0 += op.do_vectorized(index);
1221 r1 += op.do_vectorized(index + nvecs);
1222 }
1223 r0 += r1;
1224 r0.store(&outer_results[start_irreg]);
1225 // update n_chunks to denote unused element in outer_results[] from
1226 // which we can keep writing.
1227 n_chunks = start_irreg + VectorizedArray<Number>::size();
1228 }
1229 }
1230
1231
1232
1233#ifdef DEAL_II_WITH_TBB
1262 template <typename Operation, typename ResultType>
1264 {
1265 static const unsigned int threshold_array_allocate = 512;
1266
1267 TBBReduceFunctor(const Operation &op,
1268 const size_type start,
1269 const size_type end)
1270 : op(op)
1271 , start(start)
1272 , end(end)
1273 {
1274 const size_type vec_size = end - start;
1275 // set chunk size for sub-tasks
1276 const unsigned int gs =
1278 n_chunks =
1279 std::min(static_cast<size_type>(4 * MultithreadInfo::n_threads()),
1280 vec_size / gs);
1281 chunk_size = vec_size / n_chunks;
1282
1283 // round to next multiple of 512 (or leave it at the minimum grain size
1284 // if that happens to be smaller). this is advantageous because our
1285 // algorithm favors lengths of a power of 2 due to pairwise summation ->
1286 // at most one 'oddly' sized chunk
1287 if (chunk_size > 512)
1288 chunk_size = ((chunk_size + 511) / 512) * 512;
1289 n_chunks = (vec_size + chunk_size - 1) / chunk_size;
1290 AssertIndexRange((n_chunks - 1) * chunk_size, vec_size);
1291 AssertIndexRange(vec_size, n_chunks * chunk_size + 1);
1292
1294 {
1295 // make sure we allocate an even number of elements,
1296 // access to the new last element is needed in do_sum()
1297 large_array.resize(2 * ((n_chunks + 1) / 2));
1298 array_ptr = large_array.data();
1299 }
1300 else
1301 array_ptr = &small_array[0];
1302 }
1303
1308 void
1309 operator()(const tbb::blocked_range<size_type> &range) const
1310 {
1311 for (size_type i = range.begin(); i < range.end(); ++i)
1313 start + i * chunk_size,
1314 std::min(start + (i + 1) * chunk_size, end),
1315 array_ptr[i]);
1316 }
1317
1318 ResultType
1319 do_sum() const
1320 {
1321 while (n_chunks > 1)
1322 {
1323 if (n_chunks % 2 == 1)
1324 array_ptr[n_chunks++] = ResultType();
1325 for (size_type i = 0; i < n_chunks; i += 2)
1326 array_ptr[i / 2] = array_ptr[i] + array_ptr[i + 1];
1327 n_chunks /= 2;
1328 }
1329 return array_ptr[0];
1330 }
1331
1332 const Operation &op;
1335
1336 mutable unsigned int n_chunks;
1337 unsigned int chunk_size;
1339 std::vector<ResultType> large_array;
1340 // this variable either points to small_array or large_array depending on
1341 // the number of threads we want to feed
1342 mutable ResultType *array_ptr;
1343 };
1344#endif
1345
1346
1347
1352 template <typename Operation, typename ResultType>
1353 void
1355 const Operation &op,
1356 const size_type start,
1357 const size_type end,
1358 ResultType & result,
1359 const std::shared_ptr<::parallel::internal::TBBPartitioner>
1360 &partitioner)
1361 {
1362#ifdef DEAL_II_WITH_TBB
1363 const size_type vec_size = end - start;
1364 // only go to the parallel function in case there are at least 4 parallel
1365 // items, otherwise the overhead is too large
1366 if (vec_size >=
1369 {
1370 Assert(partitioner.get() != nullptr,
1372 "Unexpected initialization of Vector that does "
1373 "not set the TBB partitioner to a usable state."));
1374 std::shared_ptr<tbb::affinity_partitioner> tbb_partitioner =
1375 partitioner->acquire_one_partitioner();
1376
1377 TBBReduceFunctor<Operation, ResultType> generic_functor(op,
1378 start,
1379 end);
1380 // We use a minimum grain size of 1 here since the grains at this
1381 // stage of dividing the work refer to the number of vector chunks
1382 // that are processed by (possibly different) threads in the
1383 // parallelized for loop (i.e., they do not refer to individual
1384 // vector entries). The number of chunks here is calculated inside
1385 // TBBForFunctor. See also GitHub issue #2496 for further discussion
1386 // of this strategy.
1388 static_cast<size_type>(0),
1389 static_cast<size_type>(generic_functor.n_chunks),
1390 generic_functor,
1391 1,
1392 tbb_partitioner);
1393 partitioner->release_one_partitioner(tbb_partitioner);
1394 result = generic_functor.do_sum();
1395 }
1396 else
1397 accumulate_recursive(op, start, end, result);
1398#else
1399 accumulate_recursive(op, start, end, result);
1400 (void)partitioner;
1401#endif
1402 }
1403
1404
1405 template <typename Number, typename Number2, typename MemorySpace>
1407 {
1408 static void
1410 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1411 /*thread_loop_partitioner*/,
1412 const size_type /*size*/,
1413 const ::MemorySpace::MemorySpaceData<Number2, MemorySpace>
1414 & /*v_data*/,
1416 {
1417 static_assert(
1418 std::is_same<MemorySpace, ::MemorySpace::CUDA>::value &&
1419 std::is_same<Number, Number2>::value,
1420 "For the CUDA MemorySpace Number and Number2 should be the same type");
1421 }
1422
1423 static void
1425 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1426 /*thread_loop_partitioner*/,
1427 const size_type /*size*/,
1428 const Number /*s*/,
1430 {}
1431
1432 static void
1434 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1435 /*thread_loop_partitioner*/,
1436 const size_type /*size*/,
1437 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1438 & /*v_data*/,
1440 {}
1441
1442 static void
1444 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1445 /*thread_loop_partitioner*/,
1446 const size_type /*size*/,
1447 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1448 & /*v_data*/,
1450 {}
1451
1452 static void
1454 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1455 /*thread_loop_partitioner*/,
1456 const size_type /*size*/,
1457 Number /*a*/,
1459 {}
1460
1461 static void
1463 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1464 /*thread_loop_partitioner*/,
1465 const size_type /*size*/,
1466 const Number /*a*/,
1467 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1468 & /*v_data*/,
1470 {}
1471
1472 static void
1474 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1475 /*thread_loop_partitioner*/,
1476 const size_type /*size*/,
1477 const Number /*a*/,
1478 const Number /*b*/,
1479 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1480 & /*v_data*/,
1481 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1482 & /*w_data*/,
1484 {}
1485
1486 static void
1488 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1489 /*thread_loop_partitioner*/,
1490 const size_type /*size*/,
1491 const Number /*x*/,
1492 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1493 & /*v_data*/,
1495 {}
1496
1497 static void
1499 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1500 /*thread_loop_partitioner*/,
1501 const size_type /*size*/,
1502 const Number /*x*/,
1503 const Number /*a*/,
1504 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1505 & /*v_data*/,
1507 {}
1508
1509 static void
1511 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1512 /*thread_loop_partitioner*/,
1513 const size_type /*size*/,
1514 const Number /*x*/,
1515 const Number /*a*/,
1516 const Number /*b*/,
1517 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1518 & /*v_data*/,
1519 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1520 & /*w_data*/,
1522 {}
1523
1524 static void
1526 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1527 /*thread_loop_partitioner*/,
1528 const size_type /*size*/,
1529 const Number /*factor*/,
1531 {}
1532
1533 static void
1535 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1536 /*thread_loop_partitioner*/,
1537 const size_type /*size*/,
1538 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1539 & /*v_data*/,
1541 {}
1542
1543 static void
1545 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1546 /*thread_loop_partitioner*/,
1547 const size_type /*size*/,
1548 const Number /*a*/,
1549 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1550 & /*v_data*/,
1552 {}
1553
1554 static void
1556 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1557 /*thread_loop_partitioner*/,
1558 const size_type /*size*/,
1559 const Number /*a*/,
1560 const Number /*b*/,
1561 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1562 & /*v_data*/,
1563 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1564 & /*w_data*/,
1566 {}
1567
1568 static Number
1570 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1571 /*thread_loop_partitioner*/,
1572 const size_type /*size*/,
1573 const ::MemorySpace::MemorySpaceData<Number2, MemorySpace>
1574 & /*v_data*/,
1576 {
1577 return Number();
1578 }
1579
1580 template <typename real_type>
1581 static void
1583 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1584 /*thread_loop_partitioner*/,
1585 const size_type /*size*/,
1586 real_type & /*sum*/,
1587 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1588 & /*v_data*/,
1590 {}
1591
1592 static Number
1594 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1595 /*thread_loop_partitioner*/,
1596 const size_type /*size*/,
1597 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1598 & /*data*/)
1599 {
1600 return Number();
1601 }
1602
1603 template <typename real_type>
1604 static void
1606 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1607 /*thread_loop_partitioner*/,
1608 const size_type /*size*/,
1609 real_type & /*sum*/,
1610 Number * /*values*/,
1611 Number * /*values_dev*/)
1612 {}
1613
1614 template <typename real_type>
1615 static void
1617 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1618 /*thread_loop_partitioner*/,
1619 const size_type /*size*/,
1620 real_type & /*sum*/,
1621 real_type /*p*/,
1623 {}
1624
1625 static Number
1627 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1628 /*thread_loop_partitioner*/,
1629 const size_type /*size*/,
1630 const Number /*a*/,
1631 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1632 & /*v_data*/,
1633 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1634 & /*w_data*/,
1636 {
1637 return Number();
1638 }
1639
1640 template <typename MemorySpace2>
1641 static void
1643 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1644 /*thread_loop_partitioner*/,
1645 const size_type /*size*/,
1646 VectorOperation::values /*operation*/,
1647 const ::MemorySpace::MemorySpaceData<Number, MemorySpace2>
1648 & /*v_data*/,
1650 {}
1651 };
1652
1653
1654
1655 template <typename Number, typename Number2>
1656 struct functions<Number, Number2, ::MemorySpace::Host>
1657 {
1658 static void
1659 copy(const std::shared_ptr<::parallel::internal::TBBPartitioner>
1660 & thread_loop_partitioner,
1661 const size_type size,
1662 const ::MemorySpace::
1663 MemorySpaceData<Number2, ::MemorySpace::Host> &v_data,
1666 &data)
1667 {
1668 Vector_copy<Number, Number2> copier(v_data.values.get(),
1669 data.values.get());
1670 parallel_for(copier, 0, size, thread_loop_partitioner);
1671 }
1672
1673 static void
1674 set(const std::shared_ptr<::parallel::internal::TBBPartitioner>
1675 & thread_loop_partitioner,
1676 const size_type size,
1677 const Number s,
1680 &data)
1681 {
1682 Vector_set<Number> setter(s, data.values.get());
1683 parallel_for(setter, 0, size, thread_loop_partitioner);
1684 }
1685
1686 static void
1688 const std::shared_ptr<::parallel::internal::TBBPartitioner>
1689 & thread_loop_partitioner,
1690 const size_type size,
1691 const ::MemorySpace::
1692 MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1695 &data)
1696 {
1697 Vectorization_add_v<Number> vector_add(data.values.get(),
1698 v_data.values.get());
1699 parallel_for(vector_add, 0, size, thread_loop_partitioner);
1700 }
1701
1702 static void
1704 const std::shared_ptr<::parallel::internal::TBBPartitioner>
1705 & thread_loop_partitioner,
1706 const size_type size,
1707 const ::MemorySpace::
1708 MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1711 &data)
1712 {
1713 Vectorization_subtract_v<Number> vector_subtract(data.values.get(),
1714 v_data.values.get());
1715 parallel_for(vector_subtract, 0, size, thread_loop_partitioner);
1716 }
1717
1718 static void
1720 const std::shared_ptr<::parallel::internal::TBBPartitioner>
1721 & thread_loop_partitioner,
1722 const size_type size,
1723 Number a,
1726 &data)
1727 {
1728 Vectorization_add_factor<Number> vector_add(data.values.get(), a);
1729 parallel_for(vector_add, 0, size, thread_loop_partitioner);
1730 }
1731
1732 static void
1733 add_av(const std::shared_ptr<::parallel::internal::TBBPartitioner>
1734 & thread_loop_partitioner,
1735 const size_type size,
1736 const Number a,
1737 const ::MemorySpace::
1738 MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1741 &data)
1742 {
1743 Vectorization_add_av<Number> vector_add(data.values.get(),
1744 v_data.values.get(),
1745 a);
1746 parallel_for(vector_add, 0, size, thread_loop_partitioner);
1747 }
1748
1749 static void
1751 const std::shared_ptr<::parallel::internal::TBBPartitioner>
1752 & thread_loop_partitioner,
1753 const size_type size,
1754 const Number a,
1755 const Number b,
1756 const ::MemorySpace::
1757 MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1758 const ::MemorySpace::
1759 MemorySpaceData<Number, ::MemorySpace::Host> &w_data,
1762 &data)
1763 {
1765 data.values.get(), v_data.values.get(), w_data.values.get(), a, b);
1766 parallel_for(vector_add, 0, size, thread_loop_partitioner);
1767 }
1768
1769 static void
1771 const std::shared_ptr<::parallel::internal::TBBPartitioner>
1772 & thread_loop_partitioner,
1773 const size_type size,
1774 const Number x,
1775 const ::MemorySpace::
1776 MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1779 &data)
1780 {
1781 Vectorization_sadd_xv<Number> vector_sadd(data.values.get(),
1782 v_data.values.get(),
1783 x);
1784 parallel_for(vector_sadd, 0, size, thread_loop_partitioner);
1785 }
1786
1787 static void
1789 const std::shared_ptr<::parallel::internal::TBBPartitioner>
1790 & thread_loop_partitioner,
1791 const size_type size,
1792 const Number x,
1793 const Number a,
1794 const ::MemorySpace::
1795 MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1798 &data)
1799 {
1800 Vectorization_sadd_xav<Number> vector_sadd(data.values.get(),
1801 v_data.values.get(),
1802 a,
1803 x);
1804 parallel_for(vector_sadd, 0, size, thread_loop_partitioner);
1805 }
1806
1807 static void
1809 const std::shared_ptr<::parallel::internal::TBBPartitioner>
1810 & thread_loop_partitioner,
1811 const size_type size,
1812 const Number x,
1813 const Number a,
1814 const Number b,
1815 const ::MemorySpace::
1816 MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1817 const ::MemorySpace::
1818 MemorySpaceData<Number, ::MemorySpace::Host> &w_data,
1821 &data)
1822 {
1824 data.values.get(), v_data.values.get(), w_data.values.get(), x, a, b);
1825 parallel_for(vector_sadd, 0, size, thread_loop_partitioner);
1826 }
1827
1828 static void
1830 const std::shared_ptr<::parallel::internal::TBBPartitioner>
1831 & thread_loop_partitioner,
1832 const size_type size,
1833 const Number factor,
1836 &data)
1837 {
1838 Vectorization_multiply_factor<Number> vector_multiply(data.values.get(),
1839 factor);
1840 parallel_for(vector_multiply, 0, size, thread_loop_partitioner);
1841 }
1842
1843 static void
1844 scale(const std::shared_ptr<::parallel::internal::TBBPartitioner>
1845 & thread_loop_partitioner,
1846 const size_type size,
1847 const ::MemorySpace::
1848 MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1851 &data)
1852 {
1853 Vectorization_scale<Number> vector_scale(data.values.get(),
1854 v_data.values.get());
1855 parallel_for(vector_scale, 0, size, thread_loop_partitioner);
1856 }
1857
1858 static void
1859 equ_au(const std::shared_ptr<::parallel::internal::TBBPartitioner>
1860 & thread_loop_partitioner,
1861 const size_type size,
1862 const Number a,
1863 const ::MemorySpace::
1864 MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1867 &data)
1868 {
1869 Vectorization_equ_au<Number> vector_equ(data.values.get(),
1870 v_data.values.get(),
1871 a);
1872 parallel_for(vector_equ, 0, size, thread_loop_partitioner);
1873 }
1874
1875 static void
1877 const std::shared_ptr<::parallel::internal::TBBPartitioner>
1878 & thread_loop_partitioner,
1879 const size_type size,
1880 const Number a,
1881 const Number b,
1882 const ::MemorySpace::
1883 MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1884 const ::MemorySpace::
1885 MemorySpaceData<Number, ::MemorySpace::Host> &w_data,
1888 &data)
1889 {
1891 data.values.get(), v_data.values.get(), w_data.values.get(), a, b);
1892 parallel_for(vector_equ, 0, size, thread_loop_partitioner);
1893 }
1894
1895 static Number
1896 dot(const std::shared_ptr<::parallel::internal::TBBPartitioner>
1897 & thread_loop_partitioner,
1898 const size_type size,
1899 const ::MemorySpace::
1900 MemorySpaceData<Number2, ::MemorySpace::Host> &v_data,
1903 &data)
1904 {
1905 Number sum;
1907 data.values.get(), v_data.values.get());
1909 dot, 0, size, sum, thread_loop_partitioner);
1911
1912 return sum;
1913 }
1914
1915 template <typename real_type>
1916 static void
1917 norm_2(const std::shared_ptr<::parallel::internal::TBBPartitioner>
1918 & thread_loop_partitioner,
1919 const size_type size,
1920 real_type & sum,
1923 &data)
1924 {
1925 Norm2<Number, real_type> norm2(data.values.get());
1926 parallel_reduce(norm2, 0, size, sum, thread_loop_partitioner);
1927 }
1928
1929 static Number
1931 const std::shared_ptr<::parallel::internal::TBBPartitioner>
1932 & thread_loop_partitioner,
1933 const size_type size,
1934 const ::MemorySpace::
1935 MemorySpaceData<Number, ::MemorySpace::Host> &data)
1936 {
1937 Number sum;
1938 MeanValue<Number> mean(data.values.get());
1939 parallel_reduce(mean, 0, size, sum, thread_loop_partitioner);
1940
1941 return sum;
1942 }
1943
1944 template <typename real_type>
1945 static void
1946 norm_1(const std::shared_ptr<::parallel::internal::TBBPartitioner>
1947 & thread_loop_partitioner,
1948 const size_type size,
1949 real_type & sum,
1952 &data)
1953 {
1954 Norm1<Number, real_type> norm1(data.values.get());
1955 parallel_reduce(norm1, 0, size, sum, thread_loop_partitioner);
1956 }
1957
1958 template <typename real_type>
1959 static void
1960 norm_p(const std::shared_ptr<::parallel::internal::TBBPartitioner>
1961 & thread_loop_partitioner,
1962 const size_type size,
1963 real_type & sum,
1964 const real_type p,
1967 &data)
1968 {
1969 NormP<Number, real_type> normp(data.values.get(), p);
1970 parallel_reduce(normp, 0, size, sum, thread_loop_partitioner);
1971 }
1972
1973 static Number
1975 const std::shared_ptr<::parallel::internal::TBBPartitioner>
1976 & thread_loop_partitioner,
1977 const size_type size,
1978 const Number a,
1979 const ::MemorySpace::
1980 MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1981 const ::MemorySpace::
1982 MemorySpaceData<Number, ::MemorySpace::Host> &w_data,
1985 &data)
1986 {
1987 Number sum;
1988 AddAndDot<Number> adder(data.values.get(),
1989 v_data.values.get(),
1990 w_data.values.get(),
1991 a);
1992 parallel_reduce(adder, 0, size, sum, thread_loop_partitioner);
1993
1994 return sum;
1995 }
1996
1997 template <typename MemorySpace2>
1998 static void
2000 const std::shared_ptr<::parallel::internal::TBBPartitioner>
2001 & thread_loop_partitioner,
2002 const size_type size,
2003 VectorOperation::values operation,
2004 const ::MemorySpace::MemorySpaceData<Number, MemorySpace2>
2005 &v_data,
2008 &data,
2009 typename std::enable_if<
2010 std::is_same<MemorySpace2, ::MemorySpace::Host>::value,
2011 int>::type = 0)
2012 {
2013 if (operation == VectorOperation::insert)
2014 {
2015 copy(thread_loop_partitioner, size, v_data, data);
2016 }
2017 else if (operation == VectorOperation::add)
2018 {
2019 add_vector(thread_loop_partitioner, size, v_data, data);
2020 }
2021 else
2022 {
2024 }
2025 }
2026
2027#ifdef DEAL_II_COMPILER_CUDA_AWARE
2028 template <typename MemorySpace2>
2029 static void
2031 const std::shared_ptr<::parallel::internal::TBBPartitioner>
2032 & /*thread_loop_partitioner*/,
2033 const size_type size,
2034 VectorOperation::values operation,
2035 const ::MemorySpace::MemorySpaceData<Number, MemorySpace2>
2036 &v_data,
2039 &data,
2040 typename std::enable_if<
2041 std::is_same<MemorySpace2, ::MemorySpace::CUDA>::value,
2042 int>::type = 0)
2043 {
2044 if (operation == VectorOperation::insert)
2045 {
2046 cudaError_t cuda_error_code = cudaMemcpy(data.values.get(),
2047 v_data.values_dev.get(),
2048 size * sizeof(Number),
2049 cudaMemcpyDeviceToHost);
2050 AssertCuda(cuda_error_code);
2051 }
2052 else
2053 {
2055 }
2056 }
2057#endif
2058 };
2059
2060
2061
2062#ifdef DEAL_II_COMPILER_CUDA_AWARE
2063 template <typename Number>
2064 struct functions<Number, Number, ::MemorySpace::CUDA>
2065 {
2066 static const int block_size =
2068 static const int chunk_size =
2070
2071 static void
2073 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2074 const size_type size,
2075 const ::MemorySpace::
2076 MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2079 &data)
2080 {
2081 cudaError_t cuda_error_code = cudaMemcpy(data.values_dev.get(),
2082 v_data.values_dev.get(),
2083 size * sizeof(Number),
2084 cudaMemcpyDeviceToDevice);
2085 AssertCuda(cuda_error_code);
2086 }
2087
2088 static void
2089 set(const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2090 const size_type size,
2091 const Number s,
2094 &data)
2095 {
2096 const int n_blocks = 1 + size / (chunk_size * block_size);
2097 ::LinearAlgebra::CUDAWrappers::kernel::set<Number>
2098 <<<n_blocks, block_size>>>(data.values_dev.get(), s, size);
2100 }
2101
2102 static void
2104 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2105 const size_type size,
2106 const ::MemorySpace::
2107 MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2110 &data)
2111 {
2112 const int n_blocks = 1 + size / (chunk_size * block_size);
2113 ::LinearAlgebra::CUDAWrappers::kernel::add_aV<Number>
2114 <<<n_blocks, block_size>>>(data.values_dev.get(),
2115 1.,
2116 v_data.values_dev.get(),
2117 size);
2119 }
2120
2121 static void
2123 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2124 const size_type size,
2125 const ::MemorySpace::
2126 MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2129 &data)
2130 {
2131 const int n_blocks = 1 + size / (chunk_size * block_size);
2132 ::LinearAlgebra::CUDAWrappers::kernel::add_aV<Number>
2133 <<<n_blocks, block_size>>>(data.values_dev.get(),
2134 -1.,
2135 v_data.values_dev.get(),
2136 size);
2138 }
2139
2140 static void
2142 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2143 const size_type size,
2144 Number a,
2147 &data)
2148 {
2149 const int n_blocks = 1 + size / (chunk_size * block_size);
2150 ::LinearAlgebra::CUDAWrappers::kernel::vec_add<Number>
2151 <<<n_blocks, block_size>>>(data.values_dev.get(), a, size);
2153 }
2154
2155 static void
2157 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2158 const size_type size,
2159 const Number a,
2160 const ::MemorySpace::
2161 MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2164 &data)
2165 {
2166 const int n_blocks = 1 + size / (chunk_size * block_size);
2167 ::LinearAlgebra::CUDAWrappers::kernel::add_aV<Number>
2168 <<<n_blocks, block_size>>>(data.values_dev.get(),
2169 a,
2170 v_data.values_dev.get(),
2171 size);
2173 }
2174
2175 static void
2177 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2178 const size_type size,
2179 const Number a,
2180 const Number b,
2181 const ::MemorySpace::
2182 MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2183 const ::MemorySpace::
2184 MemorySpaceData<Number, ::MemorySpace::CUDA> &w_data,
2187 &data)
2188 {
2189 const int n_blocks = 1 + size / (chunk_size * block_size);
2190 ::LinearAlgebra::CUDAWrappers::kernel::add_aVbW<Number>
2191 <<<dim3(n_blocks, 1), dim3(block_size)>>>(data.values_dev.get(),
2192 a,
2193 v_data.values_dev.get(),
2194 b,
2195 w_data.values_dev.get(),
2196 size);
2198 }
2199
2200 static void
2202 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2203 const size_type size,
2204 const Number x,
2205 const ::MemorySpace::
2206 MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2209 &data)
2210 {
2211 const int n_blocks = 1 + size / (chunk_size * block_size);
2212 ::LinearAlgebra::CUDAWrappers::kernel::sadd<Number>
2213 <<<dim3(n_blocks, 1), dim3(block_size)>>>(
2214 x, data.values_dev.get(), 1., v_data.values_dev.get(), size);
2216 }
2217
2218 static void
2220 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2221 const size_type size,
2222 const Number x,
2223 const Number a,
2224 const ::MemorySpace::
2225 MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2228 &data)
2229 {
2230 const int n_blocks = 1 + size / (chunk_size * block_size);
2231 ::LinearAlgebra::CUDAWrappers::kernel::sadd<Number>
2232 <<<dim3(n_blocks, 1), dim3(block_size)>>>(
2233 x, data.values_dev.get(), a, v_data.values_dev.get(), size);
2235 }
2236
2237 static void
2239 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2240 const size_type size,
2241 const Number x,
2242 const Number a,
2243 const Number b,
2244 const ::MemorySpace::
2245 MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2246 const ::MemorySpace::
2247 MemorySpaceData<Number, ::MemorySpace::CUDA> &w_data,
2250 &data)
2251 {
2252 const int n_blocks = 1 + size / (chunk_size * block_size);
2253 ::LinearAlgebra::CUDAWrappers::kernel::sadd<Number>
2254 <<<dim3(n_blocks, 1), dim3(block_size)>>>(x,
2255 data.values_dev.get(),
2256 a,
2257 v_data.values_dev.get(),
2258 b,
2259 w_data.values_dev.get(),
2260 size);
2262 }
2263
2264 static void
2266 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2267 const size_type size,
2268 const Number factor,
2271 &data)
2272 {
2273 const int n_blocks = 1 + size / (chunk_size * block_size);
2274 ::LinearAlgebra::CUDAWrappers::kernel::vec_scale<Number>
2275 <<<n_blocks, block_size>>>(data.values_dev.get(), factor, size);
2277 }
2278
2279 static void
2281 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2282 const size_type size,
2283 const ::MemorySpace::
2284 MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2287 &data)
2288 {
2289 const int n_blocks = 1 + size / (chunk_size * block_size);
2290 ::LinearAlgebra::CUDAWrappers::kernel::scale<Number>
2291 <<<dim3(n_blocks, 1), dim3(block_size)>>>(data.values_dev.get(),
2292 v_data.values_dev.get(),
2293 size);
2295 }
2296
2297 static void
2299 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2300 const size_type size,
2301 const Number a,
2302 const ::MemorySpace::
2303 MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2306 &data)
2307 {
2308 const int n_blocks = 1 + size / (chunk_size * block_size);
2309 ::LinearAlgebra::CUDAWrappers::kernel::equ<Number>
2310 <<<dim3(n_blocks, 1), dim3(block_size)>>>(data.values_dev.get(),
2311 a,
2312 v_data.values_dev.get(),
2313 size);
2315 }
2316
2317 static void
2319 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2320 const size_type size,
2321 const Number a,
2322 const Number b,
2323 const ::MemorySpace::
2324 MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2325 const ::MemorySpace::
2326 MemorySpaceData<Number, ::MemorySpace::CUDA> &w_data,
2329 &data)
2330 {
2331 const int n_blocks = 1 + size / (chunk_size * block_size);
2332 ::LinearAlgebra::CUDAWrappers::kernel::equ<Number>
2333 <<<dim3(n_blocks, 1), dim3(block_size)>>>(data.values_dev.get(),
2334 a,
2335 v_data.values_dev.get(),
2336 b,
2337 w_data.values_dev.get(),
2338 size);
2340 }
2341
2342 static Number
2343 dot(const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2344 const size_type size,
2345 const ::MemorySpace::
2346 MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2349 &data)
2350 {
2351 Number * result_device;
2352 cudaError_t error_code = cudaMalloc(&result_device, sizeof(Number));
2353 AssertCuda(error_code);
2354 error_code = cudaMemset(result_device, 0, sizeof(Number));
2355 AssertCuda(error_code);
2356
2357 const int n_blocks = 1 + size / (chunk_size * block_size);
2359 Number,
2361 <<<dim3(n_blocks, 1), dim3(block_size)>>>(result_device,
2362 data.values_dev.get(),
2363 v_data.values_dev.get(),
2364 static_cast<unsigned int>(
2365 size));
2367
2368 // Copy the result back to the host
2369 Number result;
2370 error_code = cudaMemcpy(&result,
2371 result_device,
2372 sizeof(Number),
2373 cudaMemcpyDeviceToHost);
2374 AssertCuda(error_code);
2375 // Free the memory on the device
2376 error_code = cudaFree(result_device);
2377 AssertCuda(error_code);
2378
2379 AssertIsFinite(result);
2380
2381 return result;
2382 }
2383
2384 template <typename real_type>
2385 static void
2386 norm_2(const std::shared_ptr<::parallel::internal::TBBPartitioner>
2387 & thread_loop_partitioner,
2388 const size_type size,
2389 real_type & sum,
2392 &data)
2393 {
2394 sum = dot(thread_loop_partitioner, size, data, data);
2395 }
2396
2397 static Number
2399 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2400 const size_type size,
2401 const ::MemorySpace::
2402 MemorySpaceData<Number, ::MemorySpace::CUDA> &data)
2403 {
2404 Number * result_device;
2405 cudaError_t error_code = cudaMalloc(&result_device, sizeof(Number));
2406 AssertCuda(error_code);
2407 error_code = cudaMemset(result_device, 0, sizeof(Number));
2408
2409 const int n_blocks = 1 + size / (chunk_size * block_size);
2411 Number,
2413 <<<dim3(n_blocks, 1), dim3(block_size)>>>(result_device,
2414 data.values_dev.get(),
2415 size);
2416
2417 // Copy the result back to the host
2418 Number result;
2419 error_code = cudaMemcpy(&result,
2420 result_device,
2421 sizeof(Number),
2422 cudaMemcpyDeviceToHost);
2423 AssertCuda(error_code);
2424 // Free the memory on the device
2425 error_code = cudaFree(result_device);
2426 AssertCuda(error_code);
2427
2428 return result;
2429 }
2430
2431 template <typename real_type>
2432 static void
2434 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2435 const size_type size,
2436 real_type & sum,
2439 &data)
2440 {
2441 Number * result_device;
2442 cudaError_t error_code = cudaMalloc(&result_device, sizeof(Number));
2443 AssertCuda(error_code);
2444 error_code = cudaMemset(result_device, 0, sizeof(Number));
2445
2446 const int n_blocks = 1 + size / (chunk_size * block_size);
2448 Number,
2450 <<<dim3(n_blocks, 1), dim3(block_size)>>>(result_device,
2451 data.values_dev.get(),
2452 size);
2453
2454 // Copy the result back to the host
2455 error_code = cudaMemcpy(&sum,
2456 result_device,
2457 sizeof(Number),
2458 cudaMemcpyDeviceToHost);
2459 AssertCuda(error_code);
2460 // Free the memory on the device
2461 error_code = cudaFree(result_device);
2462 AssertCuda(error_code);
2463 }
2464
2465 template <typename real_type>
2466 static void
2468 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2469 const size_type,
2470 real_type &,
2471 real_type,
2474 {
2475 Assert(false, ExcNotImplemented());
2476 }
2477
2478 static Number
2480 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2481 const size_type size,
2482 const Number a,
2483 const ::MemorySpace::
2484 MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2485 const ::MemorySpace::
2486 MemorySpaceData<Number, ::MemorySpace::CUDA> &w_data,
2489 &data)
2490 {
2491 Number * res_d;
2492 cudaError_t error_code = cudaMalloc(&res_d, sizeof(Number));
2493 AssertCuda(error_code);
2494 error_code = cudaMemset(res_d, 0, sizeof(Number));
2495 AssertCuda(error_code);
2496
2497 const int n_blocks = 1 + size / (chunk_size * block_size);
2498 ::LinearAlgebra::CUDAWrappers::kernel::add_and_dot<Number>
2499 <<<dim3(n_blocks, 1), dim3(block_size)>>>(res_d,
2500 data.values_dev.get(),
2501 v_data.values_dev.get(),
2502 w_data.values_dev.get(),
2503 a,
2504 size);
2505
2506 Number res;
2507 error_code =
2508 cudaMemcpy(&res, res_d, sizeof(Number), cudaMemcpyDeviceToHost);
2509 AssertCuda(error_code);
2510 error_code = cudaFree(res_d);
2511
2512 return res;
2513 }
2514
2515 template <typename MemorySpace2>
2516 static void
2518 const std::shared_ptr<::parallel::internal::TBBPartitioner>
2519 & thread_loop_partitioner,
2520 const size_type size,
2521 VectorOperation::values operation,
2522 const ::MemorySpace::MemorySpaceData<Number, MemorySpace2>
2523 &v_data,
2526 &data,
2527 typename std::enable_if<
2528 std::is_same<MemorySpace2, ::MemorySpace::CUDA>::value,
2529 int>::type = 0)
2530 {
2531 if (operation == VectorOperation::insert)
2532 {
2533 copy(thread_loop_partitioner, size, v_data, data);
2534 }
2535 else if (operation == VectorOperation::add)
2536 {
2537 add_vector(thread_loop_partitioner, size, v_data, data);
2538 }
2539 else
2540 {
2542 }
2543 }
2544
2545 template <typename MemorySpace2>
2546 static void
2548 const std::shared_ptr<::parallel::internal::TBBPartitioner>
2549 & /*thread_loop_partitioner*/,
2550 const size_type size,
2551 VectorOperation::values operation,
2552 const ::MemorySpace::MemorySpaceData<Number, MemorySpace2>
2553 &v_data,
2556 &data,
2557 typename std::enable_if<
2558 std::is_same<MemorySpace2, ::MemorySpace::Host>::value,
2559 int>::type = 0)
2560 {
2561 if (operation == VectorOperation::insert)
2562 {
2563 cudaError_t cuda_error_code = cudaMemcpy(data.values_dev.get(),
2564 v_data.values.get(),
2565 size * sizeof(Number),
2566 cudaMemcpyHostToDevice);
2567 AssertCuda(cuda_error_code);
2568 }
2569 else
2570 {
2572 }
2573 }
2574 };
2575#endif
2576 } // namespace VectorOperations
2577} // namespace internal
2578
2580
2581#endif
static unsigned int n_threads()
void store(Number *ptr) const
void load(const Number *ptr)
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition: config.h:140
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
#define DEAL_II_FALLTHROUGH
Definition: config.h:171
Point< 2 > first
Definition: grid_out.cc:4587
__global__ void double_vector_reduction(Number *result, const Number *v1, const Number *v2, const size_type N)
__global__ void reduction(Number *result, const Number *v, const size_type N)
#define AssertCudaKernel()
Definition: exceptions.h:1827
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1465
#define AssertIsFinite(number)
Definition: exceptions.h:1721
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
static ::ExceptionBase & ExcInternalError()
#define AssertCuda(error_code)
Definition: exceptions.h:1772
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
constexpr int chunk_size
Definition: cuda_size.h:35
constexpr int block_size
Definition: cuda_size.h:29
static const char U
static const char T
std::enable_if< IsBlockVector< VectorType >::value, unsignedint >::type n_blocks(const VectorType &vector)
Definition: operators.h:49
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm &mpi_communicator)
unsigned int minimum_parallel_grain_size
Definition: parallel.cc:34
void accumulate_recursive(const Operation &op, const size_type first, const size_type last, ResultType &result)
void parallel_reduce(const Operation &op, const size_type start, const size_type end, ResultType &result, const std::shared_ptr<::parallel::internal::TBBPartitioner > &partitioner)
void copy(const T *begin, const T *end, U *dest)
void parallel_for(Functor &functor, const size_type start, const size_type end, const std::shared_ptr<::parallel::internal::TBBPartitioner > &partitioner)
void copy(const std::complex< T > *, const std::complex< T > *, U *)
const unsigned int vector_accumulation_recursion_threshold
void accumulate_regular(const Operation &op, const size_type &n_chunks, size_type &index, ResultType(&outer_results)[vector_accumulation_recursion_threshold], std::integral_constant< bool, false >)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index
Definition: types.h:76
AddAndDot(Number *const X, const Number *const V, const Number *const W, const Number a)
Number operator()(const size_type i) const
VectorizedArray< Number > do_vectorized(const size_type i) const
Dot(const Number *const X, const Number2 *const Y)
Number operator()(const size_type i) const
VectorizedArray< Number > do_vectorized(const size_type i) const
Number operator()(const size_type i) const
VectorizedArray< Number > do_vectorized(const size_type i) const
RealType operator()(const size_type i) const
VectorizedArray< Number > do_vectorized(const size_type i) const
RealType operator()(const size_type i) const
VectorizedArray< Number > do_vectorized(const size_type i) const
RealType operator()(const size_type i) const
VectorizedArray< Number > do_vectorized(const size_type i) const
void operator()(const tbb::blocked_range< size_type > &range) const
TBBForFunctor(Functor &functor, const size_type start, const size_type end)
TBBReduceFunctor(const Operation &op, const size_type start, const size_type end)
ResultType small_array[threshold_array_allocate]
void operator()(const tbb::blocked_range< size_type > &range) const
Vector_copy(const OtherNumber *const src, Number *const dst)
void operator()(const size_type begin, const size_type end) const
Vector_set(const Number value, Number *const dst)
void operator()(const size_type begin, const size_type end) const
Vectorization_add_av(Number *const val, const Number *const v_val, const Number factor)
void operator()(const size_type begin, const size_type end) const
Vectorization_add_avpbw(Number *const val, const Number *const v_val, const Number *const w_val, const Number a, const Number b)
void operator()(const size_type begin, const size_type end) const
void operator()(const size_type begin, const size_type end) const
Vectorization_add_factor(Number *const val, const Number factor)
void operator()(const size_type begin, const size_type end) const
Vectorization_add_v(Number *const val, const Number *const v_val)
Vectorization_equ_au(Number *const val, const Number *const u_val, const Number a)
void operator()(const size_type begin, const size_type end) const
Vectorization_equ_aubv(Number *const val, const Number *const u_val, const Number *const v_val, const Number a, const Number b)
void operator()(const size_type begin, const size_type end) const
Vectorization_equ_aubvcw(Number *val, const Number *u_val, const Number *v_val, const Number *w_val, const Number a, const Number b, const Number c)
void operator()(const size_type begin, const size_type end) const
Vectorization_multiply_factor(Number *const val, const Number factor)
void operator()(const size_type begin, const size_type end) const
void operator()(const size_type begin, const size_type end) const
Vectorization_ratio(Number *val, const Number *a_val, const Number *b_val)
Vectorization_sadd_xav(Number *val, const Number *const v_val, const Number a, const Number x)
void operator()(const size_type begin, const size_type end) const
void operator()(const size_type begin, const size_type end) const
Vectorization_sadd_xavbw(Number *val, const Number *v_val, const Number *w_val, Number x, Number a, Number b)
Vectorization_sadd_xv(Number *const val, const Number *const v_val, const Number x)
void operator()(const size_type begin, const size_type end) const
void operator()(const size_type begin, const size_type end) const
Vectorization_scale(Number *const val, const Number *const v_val)
Vectorization_subtract_v(Number *val, const Number *const v_val)
void operator()(const size_type begin, const size_type end) const
static void set(const std::shared_ptr<::parallel::internal::TBBPartitioner > &thread_loop_partitioner, const size_type size, const Number s, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data)
static void norm_2(const std::shared_ptr<::parallel::internal::TBBPartitioner > &thread_loop_partitioner, const size_type size, real_type &sum, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data)
static void norm_p(const std::shared_ptr<::parallel::internal::TBBPartitioner > &thread_loop_partitioner, const size_type size, real_type &sum, const real_type p, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data)
static void import_elements(const std::shared_ptr<::parallel::internal::TBBPartitioner > &thread_loop_partitioner, const size_type size, VectorOperation::values operation, const ::MemorySpace::MemorySpaceData< Number, MemorySpace2 > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data, typename std::enable_if< std::is_same< MemorySpace2, ::MemorySpace::Host >::value, int >::type=0)
static void add_avpbw(const std::shared_ptr<::parallel::internal::TBBPartitioner > &thread_loop_partitioner, const size_type size, const Number a, const Number b, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &v_data, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &w_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data)
static void sadd_xav(const std::shared_ptr<::parallel::internal::TBBPartitioner > &thread_loop_partitioner, const size_type size, const Number x, const Number a, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data)
static void equ_au(const std::shared_ptr<::parallel::internal::TBBPartitioner > &thread_loop_partitioner, const size_type size, const Number a, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data)
static void scale(const std::shared_ptr<::parallel::internal::TBBPartitioner > &thread_loop_partitioner, const size_type size, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data)
static Number mean_value(const std::shared_ptr<::parallel::internal::TBBPartitioner > &thread_loop_partitioner, const size_type size, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data)
static void add_factor(const std::shared_ptr<::parallel::internal::TBBPartitioner > &thread_loop_partitioner, const size_type size, Number a, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data)
static void sadd_xv(const std::shared_ptr<::parallel::internal::TBBPartitioner > &thread_loop_partitioner, const size_type size, const Number x, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data)
static void norm_1(const std::shared_ptr<::parallel::internal::TBBPartitioner > &thread_loop_partitioner, const size_type size, real_type &sum, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data)
static void copy(const std::shared_ptr<::parallel::internal::TBBPartitioner > &thread_loop_partitioner, const size_type size, const ::MemorySpace::MemorySpaceData< Number2, ::MemorySpace::Host > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data)
static void import_elements(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, VectorOperation::values operation, const ::MemorySpace::MemorySpaceData< Number, MemorySpace2 > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data, typename std::enable_if< std::is_same< MemorySpace2, ::MemorySpace::CUDA >::value, int >::type=0)
static void add_av(const std::shared_ptr<::parallel::internal::TBBPartitioner > &thread_loop_partitioner, const size_type size, const Number a, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data)
static void add_vector(const std::shared_ptr<::parallel::internal::TBBPartitioner > &thread_loop_partitioner, const size_type size, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data)
static void sadd_xavbw(const std::shared_ptr<::parallel::internal::TBBPartitioner > &thread_loop_partitioner, const size_type size, const Number x, const Number a, const Number b, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &v_data, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &w_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data)
static void subtract_vector(const std::shared_ptr<::parallel::internal::TBBPartitioner > &thread_loop_partitioner, const size_type size, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data)
static Number add_and_dot(const std::shared_ptr<::parallel::internal::TBBPartitioner > &thread_loop_partitioner, const size_type size, const Number a, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &v_data, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &w_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data)
static Number dot(const std::shared_ptr<::parallel::internal::TBBPartitioner > &thread_loop_partitioner, const size_type size, const ::MemorySpace::MemorySpaceData< Number2, ::MemorySpace::Host > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data)
static void equ_aubv(const std::shared_ptr<::parallel::internal::TBBPartitioner > &thread_loop_partitioner, const size_type size, const Number a, const Number b, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &v_data, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &w_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data)
static void multiply_factor(const std::shared_ptr<::parallel::internal::TBBPartitioner > &thread_loop_partitioner, const size_type size, const Number factor, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data)
static void add_factor(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, Number a, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &data)
static Number add_and_dot(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const Number a, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &v_data, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &w_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &data)
static void multiply_factor(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const Number factor, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &data)
static void subtract_vector(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &data)
static void scale(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &data)
static void copy(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &data)
static Number mean_value(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &data)
static void norm_p(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type, real_type &, real_type, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &)
static void sadd_xavbw(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const Number x, const Number a, const Number b, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &v_data, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &w_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &data)
static void equ_aubv(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const Number a, const Number b, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &v_data, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &w_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &data)
static void sadd_xav(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const Number x, const Number a, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &data)
static Number dot(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &data)
static void norm_1(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, real_type &sum, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &data)
static void sadd_xv(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const Number x, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &data)
static void import_elements(const std::shared_ptr<::parallel::internal::TBBPartitioner > &thread_loop_partitioner, const size_type size, VectorOperation::values operation, const ::MemorySpace::MemorySpaceData< Number, MemorySpace2 > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &data, typename std::enable_if< std::is_same< MemorySpace2, ::MemorySpace::CUDA >::value, int >::type=0)
static void equ_au(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const Number a, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &data)
static void norm_2(const std::shared_ptr<::parallel::internal::TBBPartitioner > &thread_loop_partitioner, const size_type size, real_type &sum, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &data)
static void import_elements(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, VectorOperation::values operation, const ::MemorySpace::MemorySpaceData< Number, MemorySpace2 > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &data, typename std::enable_if< std::is_same< MemorySpace2, ::MemorySpace::Host >::value, int >::type=0)
static void add_vector(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &data)
static void set(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const Number s, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &data)
static void add_avpbw(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const Number a, const Number b, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &v_data, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &w_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &data)
static void add_av(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const Number a, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::CUDA > &data)
static Number mean_value(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type, const ::MemorySpace::MemorySpaceData< Number, MemorySpace > &)
static void equ_au(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type, const Number, const ::MemorySpace::MemorySpaceData< Number, MemorySpace > &, ::MemorySpace::MemorySpaceData< Number, MemorySpace > &)
static void add_avpbw(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type, const Number, const Number, const ::MemorySpace::MemorySpaceData< Number, MemorySpace > &, const ::MemorySpace::MemorySpaceData< Number, MemorySpace > &, ::MemorySpace::MemorySpaceData< Number, MemorySpace > &)
static void sadd_xv(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type, const Number, const ::MemorySpace::MemorySpaceData< Number, MemorySpace > &, ::MemorySpace::MemorySpaceData< Number, MemorySpace > &)
static void add_factor(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type, Number, ::MemorySpace::MemorySpaceData< Number, MemorySpace > &)
static void norm_2(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type, real_type &, const ::MemorySpace::MemorySpaceData< Number, MemorySpace > &, ::MemorySpace::MemorySpaceData< Number, MemorySpace > &)
static void add_vector(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type, const ::MemorySpace::MemorySpaceData< Number, MemorySpace > &, ::MemorySpace::MemorySpaceData< Number, MemorySpace > &)
static void scale(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type, const ::MemorySpace::MemorySpaceData< Number, MemorySpace > &, ::MemorySpace::MemorySpaceData< Number, MemorySpace > &)
static void subtract_vector(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type, const ::MemorySpace::MemorySpaceData< Number, MemorySpace > &, ::MemorySpace::MemorySpaceData< Number, MemorySpace > &)
static void sadd_xavbw(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type, const Number, const Number, const Number, const ::MemorySpace::MemorySpaceData< Number, MemorySpace > &, const ::MemorySpace::MemorySpaceData< Number, MemorySpace > &, ::MemorySpace::MemorySpaceData< Number, MemorySpace > &)
static Number dot(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type, const ::MemorySpace::MemorySpaceData< Number2, MemorySpace > &, ::MemorySpace::MemorySpaceData< Number, MemorySpace > &)
static Number add_and_dot(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type, const Number, const ::MemorySpace::MemorySpaceData< Number, MemorySpace > &, const ::MemorySpace::MemorySpaceData< Number, MemorySpace > &, ::MemorySpace::MemorySpaceData< Number, MemorySpace > &)
static void sadd_xav(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type, const Number, const Number, const ::MemorySpace::MemorySpaceData< Number, MemorySpace > &, ::MemorySpace::MemorySpaceData< Number, MemorySpace > &)
static void import_elements(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type, VectorOperation::values, const ::MemorySpace::MemorySpaceData< Number, MemorySpace2 > &, ::MemorySpace::MemorySpaceData< Number, MemorySpace > &)
static void copy(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type, const ::MemorySpace::MemorySpaceData< Number2, MemorySpace > &, ::MemorySpace::MemorySpaceData< Number, MemorySpace > &)
static void norm_1(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type, real_type &, Number *, Number *)
static void add_av(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type, const Number, const ::MemorySpace::MemorySpaceData< Number, MemorySpace > &, ::MemorySpace::MemorySpaceData< Number, MemorySpace > &)
static void multiply_factor(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type, const Number, ::MemorySpace::MemorySpaceData< Number, MemorySpace > &)
static void norm_p(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type, real_type &, real_type, ::MemorySpace::MemorySpaceData< Number, MemorySpace > &)
static void equ_aubv(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type, const Number, const Number, const ::MemorySpace::MemorySpaceData< Number, MemorySpace > &, const ::MemorySpace::MemorySpaceData< Number, MemorySpace > &, ::MemorySpace::MemorySpaceData< Number, MemorySpace > &)
static void set(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type, const Number, ::MemorySpace::MemorySpaceData< Number, MemorySpace > &)
static constexpr std::enable_if< std::is_same< Dummy, number >::value &&is_cuda_compatible< Dummy >::value, real_type >::type abs_square(const number &x)
Definition: numbers.h:577
static real_type abs(const number &x)
Definition: numbers.h:599