deal.II version GIT relicensing-2652-g41e7496fd4 2025-02-17 19:10:00+00:00
\(\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\}}\)
Loading...
Searching...
No Matches
vector_operations_internal.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2016 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
16#ifndef dealii_vector_operations_internal_h
17#define dealii_vector_operations_internal_h
18
19#include <deal.II/base/config.h>
20
25#include <deal.II/base/types.h>
27
29
30#include <cstdio>
31#include <cstring>
32
33#ifdef DEAL_II_WITH_TBB
34# include <tbb/blocked_range.h>
35# include <tbb/partitioner.h>
36#endif
37
38
40
41namespace internal
42{
43 namespace VectorOperations
44 {
46
47 template <typename T>
48 bool
49 is_non_negative(const T &t)
50 {
51 return t >= 0;
52 }
53
54
55 template <typename T>
56 bool
57 is_non_negative(const std::complex<T> &)
58 {
59 Assert(false, ExcMessage("Complex numbers do not have an ordering."));
60
61 return false;
62 }
63
64
65 // call std::copy, except for in
66 // the case where we want to copy
67 // from std::complex to a
68 // non-complex type
69 template <typename T, typename U>
70 void
71 copy(const T *begin, const T *end, U *dest)
72 {
73 std::copy(begin, end, dest);
74 }
75
76 template <typename T, typename U>
77 void
78 copy(const std::complex<T> *begin,
79 const std::complex<T> *end,
80 std::complex<U> *dest)
81 {
82 std::copy(begin, end, dest);
83 }
84
85 template <typename T, typename U>
86 void
87 copy(const std::complex<T> *, const std::complex<T> *, U *)
88 {
89 Assert(false,
90 ExcMessage("Can't convert a vector of complex numbers "
91 "into a vector of reals/doubles"));
92 }
93
94
95
96#ifdef DEAL_II_WITH_TBB
105 template <typename Functor>
107 {
109 const size_type start,
110 const size_type end)
112 , start(start)
113 , end(end)
114 {
115 const size_type vec_size = end - start;
116 // set chunk size for sub-tasks
117 const unsigned int gs =
119 n_chunks =
121 vec_size / gs);
123
124 // round to next multiple of 512 (or minimum grain size if that happens
125 // to be smaller). this is advantageous because our accumulation
126 // algorithms favor lengths of a power of 2 due to pairwise summation ->
127 // at most one 'oddly' sized chunk
128 if (chunk_size > 512)
129 chunk_size = ((chunk_size + 511) / 512) * 512;
133 }
134
135 void
136 operator()(const tbb::blocked_range<size_type> &range) const
137 {
138 const size_type r_begin = start + range.begin() * chunk_size;
139 const size_type r_end = std::min(start + range.end() * chunk_size, end);
141 }
142
146 unsigned int n_chunks;
148 };
149#endif
150
151 template <typename Functor>
152 void
154 Functor &functor,
155 const size_type start,
156 const size_type end,
157 const std::shared_ptr<::parallel::internal::TBBPartitioner>
158 &partitioner)
159 {
160#ifdef DEAL_II_WITH_TBB
161 const size_type vec_size = end - start;
162 // only go to the parallel function in case there are at least 4 parallel
163 // items, otherwise the overhead is too large
164 if (vec_size >=
167 {
168 Assert(partitioner.get() != nullptr,
170 "Unexpected initialization of Vector that does "
171 "not set the TBB partitioner to a usable state."));
172 std::shared_ptr<tbb::affinity_partitioner> tbb_partitioner =
173 partitioner->acquire_one_partitioner();
174
175 TBBForFunctor<Functor> generic_functor(functor, start, end);
176 // We use a minimum grain size of 1 here since the grains at this
177 // stage of dividing the work refer to the number of vector chunks
178 // that are processed by (possibly different) threads in the
179 // parallelized for loop (i.e., they do not refer to individual
180 // vector entries). The number of chunks here is calculated inside
181 // TBBForFunctor. See also GitHub issue #2496 for further discussion
182 // of this strategy.
184 static_cast<size_type>(0),
185 static_cast<size_type>(generic_functor.n_chunks),
187 1,
189 partitioner->release_one_partitioner(tbb_partitioner);
190 }
191 else if (vec_size > 0)
192 functor(start, end);
193#else
194 functor(start, end);
195 (void)partitioner;
196#endif
197 }
198
199
200 // Define the functors necessary to use SIMD with TBB. we also include the
201 // simple copy and set operations
202
203 template <typename Number>
205 {
206 Vector_set(const Number value, Number *const dst)
207 : value(value)
208 , dst(dst)
209 {
210 Assert(dst != nullptr, ExcInternalError());
211 }
212
213 void
214 operator()(const size_type begin, const size_type end) const
215 {
216 Assert(end >= begin, ExcInternalError());
217
218 if (value == Number())
219 std::fill(dst + begin, dst + end, Number());
220 else
221 std::fill(dst + begin, dst + end, value);
222 }
223
224 const Number value;
225 Number *const dst;
226 };
227
228 template <typename Number, typename OtherNumber>
230 {
231 Vector_copy(const OtherNumber *const src, Number *const dst)
232 : src(src)
233 , dst(dst)
234 {
235 Assert(src != nullptr, ExcInternalError());
236 Assert(dst != nullptr, ExcInternalError());
237 }
238
239 void
240 operator()(const size_type begin, const size_type end) const
241 {
242 Assert(end >= begin, ExcInternalError());
243
244 if constexpr (std::is_trivially_copyable<Number>() &&
245 std::is_same_v<Number, OtherNumber>)
246 std::memcpy(dst + begin, src + begin, (end - begin) * sizeof(Number));
247 else
248 {
250 for (size_type i = begin; i < end; ++i)
251 dst[i] = src[i];
252 }
253 }
254
255 const OtherNumber *const src;
256 Number *const dst;
257 };
258
259 template <typename Number>
261 {
262 Vectorization_multiply_factor(Number *const val, const Number factor)
263 : val(val)
264 , stored_factor(factor)
265 {}
266
267 void
268 operator()(const size_type begin, const size_type end) const
269 {
270 // create a local copy of the variable to help the compiler with the
271 // aliasing analysis
272 const Number factor = stored_factor;
273
275 {
277 for (size_type i = begin; i < end; ++i)
278 val[i] *= factor;
279 }
280 else
281 {
282 for (size_type i = begin; i < end; ++i)
283 val[i] *= factor;
284 }
285 }
286
287 Number *const val;
288 const Number stored_factor;
289 };
290
291 template <typename Number>
293 {
295 const Number *const v_val,
296 const Number factor)
297 : val(val)
298 , v_val(v_val)
299 , stored_factor(factor)
300 {}
301
302 void
303 operator()(const size_type begin, const size_type end) const
304 {
305 // create a local copy of the variable to help the compiler with the
306 // aliasing analysis
307 const Number factor = stored_factor;
309 {
311 for (size_type i = begin; i < end; ++i)
312 val[i] += factor * v_val[i];
313 }
314 else
315 {
316 for (size_type i = begin; i < end; ++i)
317 val[i] += factor * v_val[i];
318 }
319 }
320
321 Number *const val;
322 const Number *const v_val;
323 const Number stored_factor;
324 };
325
326 template <typename Number>
328 {
330 const Number *const v_val,
331 const Number a,
332 const Number x)
333 : val(val)
334 , v_val(v_val)
335 , stored_a(a)
336 , stored_x(x)
337 {}
338
339 void
340 operator()(const size_type begin, const size_type end) const
341 {
342 // create a local copy of the variable to help the compiler with the
343 // aliasing analysis
344 const Number x = stored_x, a = stored_a;
345
347 {
349 for (size_type i = begin; i < end; ++i)
350 val[i] = x * val[i] + a * v_val[i];
351 }
352 else
353 {
354 for (size_type i = begin; i < end; ++i)
355 val[i] = x * val[i] + a * v_val[i];
356 }
357 }
358
359 Number *const val;
360 const Number *const v_val;
361 const Number stored_a;
362 const Number stored_x;
363 };
364
365 template <typename Number>
367 {
368 Vectorization_subtract_v(Number *val, const Number *const v_val)
369 : val(val)
370 , v_val(v_val)
371 {}
372
373 void
374 operator()(const size_type begin, const size_type end) const
375 {
377 {
379 for (size_type i = begin; i < end; ++i)
380 val[i] -= v_val[i];
381 }
382 else
383 {
384 for (size_type i = begin; i < end; ++i)
385 val[i] -= v_val[i];
386 }
387 }
388
389 Number *const val;
390 const Number *const v_val;
391 };
392
393 template <typename Number>
395 {
396 Vectorization_add_factor(Number *const val, const Number factor)
397 : val(val)
398 , stored_factor(factor)
399 {}
400
401 void
402 operator()(const size_type begin, const size_type end) const
403 {
404 const Number factor = stored_factor;
405
407 {
409 for (size_type i = begin; i < end; ++i)
410 val[i] += factor;
411 }
412 else
413 {
414 for (size_type i = begin; i < end; ++i)
415 val[i] += factor;
416 }
417 }
418
419 Number *const val;
420 const Number stored_factor;
421 };
422
423 template <typename Number>
425 {
426 Vectorization_add_v(Number *const val, const Number *const v_val)
427 : val(val)
428 , v_val(v_val)
429 {}
430
431 void
432 operator()(const size_type begin, const size_type end) const
433 {
435 {
437 for (size_type i = begin; i < end; ++i)
438 val[i] += v_val[i];
439 }
440 else
441 {
442 for (size_type i = begin; i < end; ++i)
443 val[i] += v_val[i];
444 }
445 }
446
447 Number *const val;
448 const Number *const v_val;
449 };
450
451 template <typename Number>
453 {
455 const Number *const v_val,
456 const Number *const w_val,
457 const Number a,
458 const Number b)
459 : val(val)
460 , v_val(v_val)
461 , w_val(w_val)
462 , stored_a(a)
463 , stored_b(b)
464 {}
465
466 void
467 operator()(const size_type begin, const size_type end) const
468 {
469 const Number a = stored_a, b = stored_b;
470
472 {
474 for (size_type i = begin; i < end; ++i)
475 val[i] = val[i] + a * v_val[i] + b * w_val[i];
476 }
477 else
478 {
479 for (size_type i = begin; i < end; ++i)
480 val[i] = val[i] + a * v_val[i] + b * w_val[i];
481 }
482 }
483
484 Number *const val;
485 const Number *const v_val;
486 const Number *const w_val;
487 const Number stored_a;
488 const Number stored_b;
489 };
490
491 template <typename Number>
493 {
495 const Number *const v_val,
496 const Number x)
497 : val(val)
498 , v_val(v_val)
499 , stored_x(x)
500 {}
501
502 void
503 operator()(const size_type begin, const size_type end) const
504 {
505 const Number x = stored_x;
506
508 {
510 for (size_type i = begin; i < end; ++i)
511 val[i] = x * val[i] + v_val[i];
512 }
513 else
514 {
515 for (size_type i = begin; i < end; ++i)
516 val[i] = x * val[i] + v_val[i];
517 }
518 }
519
520 Number *const val;
521 const Number *const v_val;
522 const Number stored_x;
523 };
524
525 template <typename Number>
527 {
529 const Number *v_val,
530 const Number *w_val,
531 Number x,
532 Number a,
533 Number b)
534 : val(val)
535 , v_val(v_val)
536 , w_val(w_val)
537 , stored_x(x)
538 , stored_a(a)
539 , stored_b(b)
540 {}
541
542 void
543 operator()(const size_type begin, const size_type end) const
544 {
545 const Number x = stored_x, a = stored_a, b = stored_b;
546
548 {
550 for (size_type i = begin; i < end; ++i)
551 val[i] = x * val[i] + a * v_val[i] + b * w_val[i];
552 }
553 else
554 {
555 for (size_type i = begin; i < end; ++i)
556 val[i] = x * val[i] + a * v_val[i] + b * w_val[i];
557 }
558 }
559
560 Number *const val;
561 const Number *const v_val;
562 const Number *const w_val;
563 const Number stored_x;
564 const Number stored_a;
565 const Number stored_b;
566 };
567
568 template <typename Number>
570 {
571 Vectorization_scale(Number *const val, const Number *const v_val)
572 : val(val)
573 , v_val(v_val)
574 {}
575
576 void
577 operator()(const size_type begin, const size_type end) const
578 {
580 {
582 for (size_type i = begin; i < end; ++i)
583 val[i] *= v_val[i];
584 }
585 else
586 {
587 for (size_type i = begin; i < end; ++i)
588 val[i] *= v_val[i];
589 }
590 }
591
592 Number *const val;
593 const Number *const v_val;
594 };
595
596 template <typename Number>
598 {
600 const Number *const u_val,
601 const Number a)
602 : val(val)
603 , u_val(u_val)
604 , stored_a(a)
605 {}
606
607 void
608 operator()(const size_type begin, const size_type end) const
609 {
610 const Number a = stored_a;
611
613 {
615 for (size_type i = begin; i < end; ++i)
616 val[i] = a * u_val[i];
617 }
618 else
619 {
620 for (size_type i = begin; i < end; ++i)
621 val[i] = a * u_val[i];
622 }
623 }
624
625 Number *const val;
626 const Number *const u_val;
627 const Number stored_a;
628 };
629
630 template <typename Number>
632 {
634 const Number *const u_val,
635 const Number *const v_val,
636 const Number a,
637 const Number b)
638 : val(val)
639 , u_val(u_val)
640 , v_val(v_val)
641 , stored_a(a)
642 , stored_b(b)
643 {}
644
645 void
646 operator()(const size_type begin, const size_type end) const
647 {
648 const Number a = stored_a, b = stored_b;
649
651 {
653 for (size_type i = begin; i < end; ++i)
654 val[i] = a * u_val[i] + b * v_val[i];
655 }
656 else
657 {
658 for (size_type i = begin; i < end; ++i)
659 val[i] = a * u_val[i] + b * v_val[i];
660 }
661 }
662
663 Number *const val;
664 const Number *const u_val;
665 const Number *const v_val;
666 const Number stored_a;
667 const Number stored_b;
668 };
669
670 template <typename Number>
672 {
674 const Number *u_val,
675 const Number *v_val,
676 const Number *w_val,
677 const Number a,
678 const Number b,
679 const Number c)
680 : val(val)
681 , u_val(u_val)
682 , v_val(v_val)
683 , w_val(w_val)
684 , stored_a(a)
685 , stored_b(b)
686 , stored_c(c)
687 {}
688
689 void
690 operator()(const size_type begin, const size_type end) const
691 {
692 const Number a = stored_a, b = stored_b, c = stored_c;
693
695 {
697 for (size_type i = begin; i < end; ++i)
698 val[i] = a * u_val[i] + b * v_val[i] + c * w_val[i];
699 }
700 else
701 {
702 for (size_type i = begin; i < end; ++i)
703 val[i] = a * u_val[i] + b * v_val[i] + c * w_val[i];
704 }
705 }
706
707 Number *const val;
708 const Number *const u_val;
709 const Number *const v_val;
710 const Number *const w_val;
711 const Number stored_a;
712 const Number stored_b;
713 const Number stored_c;
714 };
715
716 template <typename Number>
718 {
719 Vectorization_ratio(Number *val, const Number *a_val, const Number *b_val)
720 : val(val)
721 , a_val(a_val)
722 , b_val(b_val)
723 {}
724
725 void
726 operator()(const size_type begin, const size_type end) const
727 {
729 {
731 for (size_type i = begin; i < end; ++i)
732 val[i] = a_val[i] / b_val[i];
733 }
734 else
735 {
736 for (size_type i = begin; i < end; ++i)
737 val[i] = a_val[i] / b_val[i];
738 }
739 }
740
741 Number *const val;
742 const Number *const a_val;
743 const Number *const b_val;
744 };
745
746
747
748 // All sums over all the vector entries (l2-norm, inner product, etc.) are
749 // performed with the same code, using a templated operation defined
750 // here. There are always two versions defined, a standard one that covers
751 // most cases and a vectorized one which is only for equal types and float
752 // and double.
753 template <typename Number, typename Number2>
754 struct Dot
755 {
756 static constexpr bool vectorizes = std::is_same_v<Number, Number2> &&
758
759 Dot(const Number *const X, const Number2 *const Y)
760 : X(X)
761 , Y(Y)
762 {}
763
764 Number
765 operator()(const size_type i) const
766 {
767 return X[i] * Number(numbers::NumberTraits<Number2>::conjugate(Y[i]));
768 }
769
772 {
774 x.load(X + i);
775 y.load(Y + i);
776
777 // the following operation in VectorizedArray does an element-wise
778 // scalar product without taking into account complex values and
779 // the need to take the complex-conjugate of one argument. this
780 // may be a bug, but because all VectorizedArray classes only
781 // work on real scalars, it doesn't really matter very much.
782 // in any case, assert that we really don't get here for
783 // complex-valued objects
784 static_assert(numbers::NumberTraits<Number>::is_complex == false,
785 "This operation is not correctly implemented for "
786 "complex-valued objects.");
787 return x * y;
788 }
789
790 const Number *const X;
791 const Number2 *const Y;
792 };
793
794 template <typename Number, typename RealType>
795 struct Norm2
796 {
797 static const bool vectorizes = VectorizedArray<Number>::size() > 1;
798
799 Norm2(const Number *const X)
800 : X(X)
801 {}
802
804 operator()(const size_type i) const
805 {
807 }
808
811 {
813 x.load(X + i);
814 return x * x;
815 }
816
817 const Number *const X;
818 };
819
820 template <typename Number, typename RealType>
821 struct Norm1
822 {
823 static const bool vectorizes = VectorizedArray<Number>::size() > 1;
824
825 Norm1(const Number *X)
826 : X(X)
827 {}
828
830 operator()(const size_type i) const
831 {
833 }
834
837 {
839 x.load(X + i);
840 return std::abs(x);
841 }
842
843 const Number *X;
844 };
845
846 template <typename Number, typename RealType>
847 struct NormP
848 {
849 static const bool vectorizes = VectorizedArray<Number>::size() > 1;
850
851 NormP(const Number *X, RealType p)
852 : X(X)
853 , p(p)
854 {}
855
857 operator()(const size_type i) const
858 {
860 }
861
864 {
866 x.load(X + i);
867 return std::pow(std::abs(x), p);
868 }
869
870 const Number *X;
871 const RealType p;
872 };
873
874 template <typename Number>
876 {
877 static const bool vectorizes = VectorizedArray<Number>::size() > 1;
878
879 MeanValue(const Number *X)
880 : X(X)
881 {}
882
883 Number
884 operator()(const size_type i) const
885 {
886 return X[i];
887 }
888
891 {
893 x.load(X + i);
894 return x;
895 }
896
897 const Number *X;
898 };
899
900 template <typename Number>
902 {
903 static const bool vectorizes = VectorizedArray<Number>::size() > 1;
904
905 AddAndDot(Number *const X,
906 const Number *const V,
907 const Number *const W,
908 const Number a)
909 : X(X)
910 , V(V)
911 , W(W)
912 , a(a)
913 {}
914
915 Number
916 operator()(const size_type i) const
917 {
918 X[i] += a * V[i];
919 return X[i] * Number(numbers::NumberTraits<Number>::conjugate(W[i]));
920 }
921
924 {
926 x.load(X + i);
927 v.load(V + i);
928 x += a * v;
929 x.store(X + i);
930 // may only load from W after storing in X because the pointers might
931 // point to the same memory
932 w.load(W + i);
933
934 // the following operation in VectorizedArray does an element-wise
935 // scalar product without taking into account complex values and
936 // the need to take the complex-conjugate of one argument. this
937 // may be a bug, but because all VectorizedArray classes only
938 // work on real scalars, it doesn't really matter very much.
939 // in any case, assert that we really don't get here for
940 // complex-valued objects
941 static_assert(numbers::NumberTraits<Number>::is_complex == false,
942 "This operation is not correctly implemented for "
943 "complex-valued objects.");
944 return x * w;
945 }
946
947 Number *const X;
948 const Number *const V;
949 const Number *const W;
950 const Number a;
951 };
952
953
954
955 // this is the main working loop for all vector sums using the templated
956 // operation above. it accumulates the sums using a block-wise summation
957 // algorithm with post-update. this blocked algorithm has been proposed in
958 // a similar form by Castaldo, Whaley and Chronopoulos (SIAM
959 // J. Sci. Comput. 31, 1156-1174, 2008) and we use the smallest possible
960 // block size, 2. Sometimes it is referred to as pairwise summation. The
961 // worst case error made by this algorithm is on the order O(eps *
962 // log2(vec_size)), whereas a naive summation is O(eps * vec_size). Even
963 // though the Kahan summation is even more accurate with an error O(eps)
964 // by carrying along remainders not captured by the main sum, that involves
965 // additional costs which are not worthwhile. See the Wikipedia article on
966 // the Kahan summation algorithm.
967
968 // The algorithm implemented here has the additional benefit that it is
969 // easily parallelized without changing the order of how the elements are
970 // added (floating point addition is not associative). For the same vector
971 // size and minimum_parallel_grainsize, the blocks are always the
972 // same and added pairwise.
973
974 // The depth of recursion is controlled by the 'magic' parameter
975 // vector_accumulation_recursion_threshold: If the length is below
976 // vector_accumulation_recursion_threshold * 32 (32 is the part of code we
977 // unroll), a straight loop instead of recursion will be used. At the
978 // innermost level, eight values are added consecutively in order to better
979 // balance multiplications and additions.
980
981 // Loops are unrolled as follows: the range [first,last) is broken into
982 // @p n_chunks each of size 32 plus the @p remainder.
983 // accumulate_regular() does the work on 32*n_chunks elements employing SIMD
984 // if possible and stores the result of the operation for each chunk in @p outer_results.
985
986 // The code returns the result as the last argument in order to make
987 // spawning tasks simpler and use automatic template deduction.
988
989
996
997 template <typename Operation, typename ResultType>
998 void
1000 const size_type first,
1001 const size_type last,
1002 ResultType &result)
1003 {
1004 if (first == last)
1005 {
1006 result = ResultType();
1007 return;
1008 }
1009
1010 const size_type vec_size = last - first;
1012 {
1013 // The vector is short enough so we perform the summation. We store
1014 // the number of chunks (each 32 indices) for the given vector
1015 // length; all results are stored in outer_results[0,n_chunks). We
1016 // keep twice the number around to be able to do the pairwise
1017 // summation with a single for loop (see the loop over j below)
1019
1020 // Select between the regular version and vectorized version based
1021 // on the number types we are given. To choose the vectorized
1022 // version often enough, we need to have all tasks but the last one
1023 // to be divisible by the vectorization length
1024 size_type n_chunks =
1025 do_accumulate(op,
1026 vec_size,
1027 first,
1029 std::bool_constant<Operation::vectorizes>());
1030
1031 AssertIndexRange(n_chunks,
1033
1034 // now sum the results from the chunks stored in
1035 // outer_results[0,n_chunks) recursively
1036 unsigned int j = 0;
1037 constexpr unsigned int n_lanes = VectorizedArray<ResultType>::size();
1038 for (; j + 2 * n_lanes - 1 < n_chunks;
1039 j += 2 * n_lanes, n_chunks += n_lanes)
1040 {
1042 a.load(outer_results + j);
1043 b.load(outer_results + j + n_lanes);
1044 a += b;
1045 a.store(outer_results + n_chunks);
1046 }
1047
1048 // In the vectorized case, we know the loop bounds and can do things
1049 // more efficiently
1050 if (Operation::vectorizes)
1051 {
1052 AssertDimension(j + n_lanes, n_chunks);
1053 AssertIndexRange(n_chunks,
1055 ResultType *result_ptr = outer_results + j;
1056 if (n_lanes >= 16)
1057 for (unsigned int i = 0; i < 8; ++i)
1058 result_ptr[i] = result_ptr[i] + result_ptr[i + 8];
1059 if (n_lanes >= 8)
1060 for (unsigned int i = 0; i < 4; ++i)
1061 result_ptr[i] = result_ptr[i] + result_ptr[i + 4];
1062 if (n_lanes >= 4)
1063 for (unsigned int i = 0; i < 2; ++i)
1064 result_ptr[i] = result_ptr[i] + result_ptr[i + 2];
1065 result = result_ptr[0] + result_ptr[1];
1066 }
1067 else
1068 {
1069 // Without vectorization, we do not know the exact bounds, so we
1070 // need to continue the variable-length pairwise summation loop
1071 // from above
1072 for (; j + 1 < n_chunks; j += 2, ++n_chunks)
1073 outer_results[n_chunks] =
1075
1076 AssertIndexRange(n_chunks,
1078 Assert(n_chunks > 0, ExcInternalError());
1079 result = outer_results[n_chunks - 1];
1080 }
1081 }
1082 else
1083 {
1084 // split vector into four pieces and work on the pieces
1085 // recursively. Make pieces (except last) divisible by one fourth the
1086 // recursion threshold.
1087 const size_type new_size =
1090 Assert(first + 3 * new_size < last, ExcInternalError());
1091 ResultType r0, r1, r2, r3;
1095 first + 2 * new_size,
1096 first + 3 * new_size,
1097 r2);
1098 accumulate_recursive(op, first + 3 * new_size, last, r3);
1099 result = (r0 + r1) + (r2 + r3);
1100 }
1101 }
1102
1103
1104 // this is the inner working routine for the accumulation loops below. We
1105 // pulled this part out of the regular accumulate routine because we might
1106 // do this thing vectorized (see specialized function below; this is the
1107 // un-vectorized version). As opposed to the vector add functions above,
1108 // we here pass the functor 'op' by value, because we cannot create a copy
1109 // of the scalar inline, and instead make sure that the numbers get local
1110 // (and thus definitely not aliased) for the compiler
1111 template <typename Operation, typename ResultType>
1112 size_type
1114 const size_type vec_size,
1115 const size_type start_index,
1116 ResultType *outer_results,
1117 std::bool_constant<false>)
1118 {
1119 // Create local copy to indicate no aliasing to the compiler
1121
1122 // choose each chunk to have a width of 32, thereby the index
1123 // is incremented by 4*8 for each @p i.
1124 size_type n_chunks = vec_size / 32;
1125 for (size_type i = 0; i < n_chunks; ++i)
1126 {
1127 ResultType r = {};
1128 for (unsigned int k = 0; k < 2; ++k)
1129 {
1130 ResultType r0 = op(index);
1131 ResultType r1 = op(index + 1);
1132 ResultType r2 = op(index + 2);
1133 ResultType r3 = op(index + 3);
1134 index += 4;
1135 for (size_type j = 1; j < 4; ++j, index += 4)
1136 {
1137 r0 += op(index);
1138 r1 += op(index + 1);
1139 r2 += op(index + 2);
1140 r3 += op(index + 3);
1141 }
1142 r += (r0 + r1) + (r2 + r3);
1143 }
1144 outer_results[i] = r;
1145 }
1146
1147 if (n_chunks * 32 < vec_size)
1148 {
1149 const size_type remainder = vec_size - n_chunks * 32;
1150 const size_type inner_chunks = remainder / 8;
1152 ResultType r0 = ResultType(), r1 = ResultType(), r2 = ResultType();
1153 switch (inner_chunks)
1154 {
1155 case 3:
1156 r2 = op(index++);
1157 for (size_type j = 1; j < 8; ++j)
1158 r2 += op(index++);
1160 case 2:
1161 r1 = op(index++);
1162 for (size_type j = 1; j < 8; ++j)
1163 r1 += op(index++);
1164 r1 += r2;
1166 case 1:
1167 r2 = op(index++);
1168 for (size_type j = 1; j < 8; ++j)
1169 r2 += op(index++);
1171 default:
1172 for (size_type j = 0; j < remainder_inner; ++j)
1173 r0 += op(index++);
1174 outer_results[n_chunks++] = (r0 + r2) + r1;
1175 break;
1176 }
1177 }
1178
1179 // make sure we worked through all indices
1181
1182 return n_chunks;
1183 }
1184
1185
1186
1187 // this is the inner working routine for the accumulation loops
1188 // below. This is the specialized case where we can vectorize. We request
1189 // the 'do_vectorized' routine of the operation instead of the regular one
1190 // which does several operations at once. As above, pass in the functor by
1191 // value to create a local copy of the scalar factors in the function (if
1192 // there are any).
1193 template <typename Operation, typename Number>
1194 size_type
1196 const size_type vec_size,
1197 const size_type start_index,
1198 Number *outer_results,
1199 std::bool_constant<true>)
1200 {
1201 // Create local copy to indicate no aliasing to the compiler
1203
1204 // we start from @p index and workout @p n_chunks each of size 32.
1205 // in order employ SIMD and work on @p nvecs at a time, we split this
1206 // loop yet again:
1207 // First we work on (n_chunks/nvecs) chunks, where each chunk processes
1208 // nvecs*(4*8) elements.
1209
1210 constexpr size_type n_lanes = VectorizedArray<Number>::size();
1211 const size_type regular_chunks = vec_size / (32 * n_lanes);
1212 for (size_type i = 0; i < regular_chunks; ++i)
1213 {
1215 for (unsigned int k = 0; k < 2; ++k)
1216 {
1217 VectorizedArray<Number> r0 = op.do_vectorized(index);
1218 VectorizedArray<Number> r1 = op.do_vectorized(index + n_lanes);
1220 op.do_vectorized(index + 2 * n_lanes);
1222 op.do_vectorized(index + 3 * n_lanes);
1223 index += n_lanes * 4;
1224 for (size_type j = 1; j < 4; ++j, index += n_lanes * 4)
1225 {
1226 r0 += op.do_vectorized(index);
1227 r1 += op.do_vectorized(index + n_lanes);
1228 r2 += op.do_vectorized(index + 2 * n_lanes);
1229 r3 += op.do_vectorized(index + 3 * n_lanes);
1230 }
1231 r += (r0 + r1) + (r2 + r3);
1232 }
1233 r.store(&outer_results[i * n_lanes]);
1234 }
1235
1236 // If we are treating a case where the vector length is not divisible by
1237 // the vectorization length, need a cleanup loop
1238 // The remaining chunks are processed one by one starting from
1239 // regular_chunks * n_lanes; We do as much as possible with 2 SIMD
1240 // operations within each chunk. Here we assume that n_lanes < 32/2 = 16
1241 // as well as 16 % n_lanes == 0.
1242 static_assert(n_lanes <= 16 && 16 % n_lanes == 0,
1243 "VectorizedArray::size() must be 1, 2, 4, 8, or 16");
1244 size_type n_chunks = regular_chunks * n_lanes;
1245 const size_type start_irregular = regular_chunks * n_lanes * 32;
1247 {
1251 const size_type loop_length = remainder / (2 * n_lanes);
1252 for (size_type j = 0; j < loop_length; ++j, index += 2 * n_lanes)
1253 {
1254 r0 += op.do_vectorized(index);
1255 r1 += op.do_vectorized(index + n_lanes);
1256 }
1257 Number scalar_part = Number();
1258 size_type last = remainder % (2 * n_lanes);
1259 if (last > 0)
1260 {
1261 if (last >= n_lanes)
1262 {
1263 r0 += op.do_vectorized(index);
1264 index += n_lanes;
1265 last -= n_lanes;
1266 }
1267 for (unsigned int i = 0; i < last; ++i)
1268 scalar_part += op(index++);
1269 }
1270
1271 r0 += r1;
1272 r0.store(&outer_results[n_chunks]);
1273 outer_results[n_chunks] += scalar_part;
1274
1275 // update n_chunks to denote range of entries to sum up in
1276 // outer_results[].
1277 n_chunks += n_lanes;
1278 }
1279
1280 // make sure we worked through all indices
1282
1283 return n_chunks;
1284 }
1285
1286
1287
1288#ifdef DEAL_II_WITH_TBB
1317 template <typename Operation, typename ResultType>
1319 {
1320 static const unsigned int threshold_array_allocate = 512;
1321
1323 const size_type start,
1324 const size_type end)
1325 : op(op)
1326 , start(start)
1327 , end(end)
1328 {
1329 const size_type vec_size = end - start;
1330 // set chunk size for sub-tasks
1331 const unsigned int gs =
1333 n_chunks =
1334 std::min(static_cast<size_type>(4 * MultithreadInfo::n_threads()),
1335 vec_size / gs);
1337
1338 // round to next multiple of 512 (or leave it at the minimum grain size
1339 // if that happens to be smaller). this is advantageous because our
1340 // algorithm favors lengths of a power of 2 due to pairwise summation ->
1341 // at most one 'oddly' sized chunk
1342 if (chunk_size > 512)
1343 chunk_size = ((chunk_size + 511) / 512) * 512;
1347
1349 {
1350 // make sure we allocate an even number of elements,
1351 // access to the new last element is needed in do_sum()
1352 large_array.resize(2 * ((n_chunks + 1) / 2));
1353 array_ptr = large_array.data();
1354 }
1355 else
1356 array_ptr = &small_array[0];
1357 }
1358
1363 void
1364 operator()(const tbb::blocked_range<size_type> &range) const
1365 {
1366 for (size_type i = range.begin(); i < range.end(); ++i)
1368 start + i * chunk_size,
1369 std::min(start + (i + 1) * chunk_size, end),
1370 array_ptr[i]);
1371 }
1372
1373 ResultType
1374 do_sum() const
1375 {
1376 while (n_chunks > 1)
1377 {
1378 if (n_chunks % 2 == 1)
1379 array_ptr[n_chunks++] = ResultType();
1380 for (size_type i = 0; i < n_chunks; i += 2)
1381 array_ptr[i / 2] = array_ptr[i] + array_ptr[i + 1];
1382 n_chunks /= 2;
1383 }
1384 return array_ptr[0];
1385 }
1386
1390
1391 mutable unsigned int n_chunks;
1392 unsigned int chunk_size;
1394 std::vector<ResultType> large_array;
1395 // this variable either points to small_array or large_array depending on
1396 // the number of threads we want to feed
1397 mutable ResultType *array_ptr;
1398 };
1399#endif
1400
1401
1402
1407 template <typename Operation, typename ResultType>
1408#ifndef DEBUG
1410#endif
1411 inline void
1413 const Operation &op,
1414 const size_type start,
1415 const size_type end,
1416 ResultType &result,
1417 const std::shared_ptr<::parallel::internal::TBBPartitioner>
1418 &partitioner)
1419 {
1420#ifdef DEAL_II_WITH_TBB
1421 const size_type vec_size = end - start;
1422 // only go to the parallel function in case there are at least 4 parallel
1423 // items, otherwise the overhead is too large
1424 if (vec_size >=
1427 {
1428 Assert(partitioner.get() != nullptr,
1430 "Unexpected initialization of Vector that does "
1431 "not set the TBB partitioner to a usable state."));
1432 std::shared_ptr<tbb::affinity_partitioner> tbb_partitioner =
1433 partitioner->acquire_one_partitioner();
1434
1436 start,
1437 end);
1438 // We use a minimum grain size of 1 here since the grains at this
1439 // stage of dividing the work refer to the number of vector chunks
1440 // that are processed by (possibly different) threads in the
1441 // parallelized for loop (i.e., they do not refer to individual
1442 // vector entries). The number of chunks here is calculated inside
1443 // TBBForFunctor. See also GitHub issue #2496 for further discussion
1444 // of this strategy.
1446 static_cast<size_type>(0),
1447 static_cast<size_type>(generic_functor.n_chunks),
1449 1,
1451 partitioner->release_one_partitioner(tbb_partitioner);
1452 result = generic_functor.do_sum();
1453 }
1454 else
1455 accumulate_recursive(op, start, end, result);
1456#else
1457 accumulate_recursive(op, start, end, result);
1458 (void)partitioner;
1459#endif
1460 }
1461
1462
1463 template <typename Number, typename Number2, typename MemorySpace>
1465 {
1466 static void
1468 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1469 /*thread_loop_partitioner*/,
1470 const size_type /*size*/,
1471 const ::MemorySpace::MemorySpaceData<Number2, MemorySpace>
1472 & /*v_data*/,
1474 {
1475 static_assert(
1476 std::is_same_v<MemorySpace, ::MemorySpace::Default> &&
1477 std::is_same_v<Number, Number2>,
1478 "For the Default MemorySpace Number and Number2 should be the same type");
1479 }
1480
1481 static void
1483 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1484 /*thread_loop_partitioner*/,
1485 const size_type /*size*/,
1486 const Number /*s*/,
1488 {}
1489
1490 static void
1492 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1493 /*thread_loop_partitioner*/,
1494 const size_type /*size*/,
1495 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1496 & /*v_data*/,
1498 {}
1499
1500 static void
1502 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1503 /*thread_loop_partitioner*/,
1504 const size_type /*size*/,
1505 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1506 & /*v_data*/,
1508 {}
1509
1510 static void
1512 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1513 /*thread_loop_partitioner*/,
1514 const size_type /*size*/,
1515 Number /*a*/,
1517 {}
1518
1519 static void
1521 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1522 /*thread_loop_partitioner*/,
1523 const size_type /*size*/,
1524 const Number /*a*/,
1525 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1526 & /*v_data*/,
1528 {}
1529
1530 static void
1532 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1533 /*thread_loop_partitioner*/,
1534 const size_type /*size*/,
1535 const Number /*a*/,
1536 const Number /*b*/,
1537 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1538 & /*v_data*/,
1539 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1540 & /*w_data*/,
1542 {}
1543
1544 static void
1546 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1547 /*thread_loop_partitioner*/,
1548 const size_type /*size*/,
1549 const Number /*x*/,
1550 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1551 & /*v_data*/,
1553 {}
1554
1555 static void
1557 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1558 /*thread_loop_partitioner*/,
1559 const size_type /*size*/,
1560 const Number /*x*/,
1561 const Number /*a*/,
1562 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1563 & /*v_data*/,
1565 {}
1566
1567 static void
1569 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1570 /*thread_loop_partitioner*/,
1571 const size_type /*size*/,
1572 const Number /*x*/,
1573 const Number /*a*/,
1574 const Number /*b*/,
1575 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1576 & /*v_data*/,
1577 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1578 & /*w_data*/,
1580 {}
1581
1582 static void
1584 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1585 /*thread_loop_partitioner*/,
1586 const size_type /*size*/,
1587 const Number /*factor*/,
1589 {}
1590
1591 static void
1593 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1594 /*thread_loop_partitioner*/,
1595 const size_type /*size*/,
1596 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1597 & /*v_data*/,
1599 {}
1600
1601 static void
1603 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1604 /*thread_loop_partitioner*/,
1605 const size_type /*size*/,
1606 const Number /*a*/,
1607 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1608 & /*v_data*/,
1610 {}
1611
1612 static void
1614 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1615 /*thread_loop_partitioner*/,
1616 const size_type /*size*/,
1617 const Number /*a*/,
1618 const Number /*b*/,
1619 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1620 & /*v_data*/,
1621 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1622 & /*w_data*/,
1624 {}
1625
1626 static Number
1628 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1629 /*thread_loop_partitioner*/,
1630 const size_type /*size*/,
1631 const ::MemorySpace::MemorySpaceData<Number2, MemorySpace>
1632 & /*v_data*/,
1634 {
1635 return Number();
1636 }
1637
1638 template <typename real_type>
1639 static void
1641 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1642 /*thread_loop_partitioner*/,
1643 const size_type /*size*/,
1644 real_type & /*sum*/,
1645 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1646 & /*v_data*/,
1648 {}
1649
1650 static Number
1652 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1653 /*thread_loop_partitioner*/,
1654 const size_type /*size*/,
1655 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1656 & /*data*/)
1657 {
1658 return Number();
1659 }
1660
1661 template <typename real_type>
1662 static void
1664 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1665 /*thread_loop_partitioner*/,
1666 const size_type /*size*/,
1667 real_type & /*sum*/,
1668 Number * /*values*/,
1669 Number * /*values*/)
1670 {}
1671
1672 template <typename real_type>
1673 static void
1675 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1676 /*thread_loop_partitioner*/,
1677 const size_type /*size*/,
1678 real_type & /*sum*/,
1679 real_type /*p*/,
1681 {}
1682
1683 static Number
1685 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1686 /*thread_loop_partitioner*/,
1687 const size_type /*size*/,
1688 const Number /*a*/,
1689 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1690 & /*v_data*/,
1691 const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1692 & /*w_data*/,
1694 {
1695 return Number();
1696 }
1697
1698 template <typename MemorySpace2>
1699 static void
1701 const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1702 /*thread_loop_partitioner*/,
1703 const size_type /*size*/,
1704 VectorOperation::values /*operation*/,
1705 const ::MemorySpace::MemorySpaceData<Number, MemorySpace2>
1706 & /*v_data*/,
1708 {}
1709 };
1710
1711
1712
1713 template <typename Number, typename Number2>
1714 struct functions<Number, Number2, ::MemorySpace::Host>
1715 {
1716 static void
1717 copy(const std::shared_ptr<::parallel::internal::TBBPartitioner>
1718 &thread_loop_partitioner,
1719 const size_type size,
1720 const ::MemorySpace::
1721 MemorySpaceData<Number2, ::MemorySpace::Host> &v_data,
1722 ::MemorySpace::MemorySpaceData<Number,
1723 ::MemorySpace::Host>
1724 &data)
1725 {
1727 data.values.data());
1728 parallel_for(copier, 0, size, thread_loop_partitioner);
1729 }
1730
1731 static void
1732 set(const std::shared_ptr<::parallel::internal::TBBPartitioner>
1733 &thread_loop_partitioner,
1734 const size_type size,
1735 const Number s,
1738 &data)
1739 {
1740 Vector_set<Number> setter(s, data.values.data());
1741 parallel_for(setter, 0, size, thread_loop_partitioner);
1742 }
1743
1744 static void
1746 const std::shared_ptr<::parallel::internal::TBBPartitioner>
1747 &thread_loop_partitioner,
1748 const size_type size,
1749 const ::MemorySpace::
1750 MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1751 ::MemorySpace::MemorySpaceData<Number,
1752 ::MemorySpace::Host>
1753 &data)
1754 {
1756 v_data.values.data());
1757 parallel_for(vector_add, 0, size, thread_loop_partitioner);
1758 }
1759
1760 static void
1762 const std::shared_ptr<::parallel::internal::TBBPartitioner>
1763 &thread_loop_partitioner,
1764 const size_type size,
1765 const ::MemorySpace::
1766 MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1767 ::MemorySpace::MemorySpaceData<Number,
1768 ::MemorySpace::Host>
1769 &data)
1770 {
1772 v_data.values.data());
1773 parallel_for(vector_subtract, 0, size, thread_loop_partitioner);
1774 }
1775
1776 static void
1778 const std::shared_ptr<::parallel::internal::TBBPartitioner>
1779 &thread_loop_partitioner,
1780 const size_type size,
1781 Number a,
1784 &data)
1785 {
1787 parallel_for(vector_add, 0, size, thread_loop_partitioner);
1788 }
1789
1790 static void
1791 add_av(const std::shared_ptr<::parallel::internal::TBBPartitioner>
1792 &thread_loop_partitioner,
1793 const size_type size,
1794 const Number a,
1795 const ::MemorySpace::
1796 MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1797 ::MemorySpace::MemorySpaceData<Number,
1798 ::MemorySpace::Host>
1799 &data)
1800 {
1802 v_data.values.data(),
1803 a);
1804 parallel_for(vector_add, 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 a,
1813 const Number b,
1814 const ::MemorySpace::
1815 MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1816 const ::MemorySpace::
1817 MemorySpaceData<Number, ::MemorySpace::Host> &w_data,
1818 ::MemorySpace::MemorySpaceData<Number,
1819 ::MemorySpace::Host>
1820 &data)
1821 {
1823 data.values.data(), v_data.values.data(), w_data.values.data(), a, b);
1824 parallel_for(vector_add, 0, size, thread_loop_partitioner);
1825 }
1826
1827 static void
1829 const std::shared_ptr<::parallel::internal::TBBPartitioner>
1830 &thread_loop_partitioner,
1831 const size_type size,
1832 const Number x,
1833 const ::MemorySpace::
1834 MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1835 ::MemorySpace::MemorySpaceData<Number,
1836 ::MemorySpace::Host>
1837 &data)
1838 {
1840 v_data.values.data(),
1841 x);
1842 parallel_for(vector_sadd, 0, size, thread_loop_partitioner);
1843 }
1844
1845 static void
1847 const std::shared_ptr<::parallel::internal::TBBPartitioner>
1848 &thread_loop_partitioner,
1849 const size_type size,
1850 const Number x,
1851 const Number a,
1852 const ::MemorySpace::
1853 MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1854 ::MemorySpace::MemorySpaceData<Number,
1855 ::MemorySpace::Host>
1856 &data)
1857 {
1859 v_data.values.data(),
1860 a,
1861 x);
1862 parallel_for(vector_sadd, 0, size, thread_loop_partitioner);
1863 }
1864
1865 static void
1867 const std::shared_ptr<::parallel::internal::TBBPartitioner>
1868 &thread_loop_partitioner,
1869 const size_type size,
1870 const Number x,
1871 const Number a,
1872 const Number b,
1873 const ::MemorySpace::
1874 MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1875 const ::MemorySpace::
1876 MemorySpaceData<Number, ::MemorySpace::Host> &w_data,
1877 ::MemorySpace::MemorySpaceData<Number,
1878 ::MemorySpace::Host>
1879 &data)
1880 {
1882 v_data.values.data(),
1883 w_data.values.data(),
1884 x,
1885 a,
1886 b);
1887 parallel_for(vector_sadd, 0, size, thread_loop_partitioner);
1888 }
1889
1890 static void
1892 const std::shared_ptr<::parallel::internal::TBBPartitioner>
1893 &thread_loop_partitioner,
1894 const size_type size,
1895 const Number factor,
1898 &data)
1899 {
1901 data.values.data(), factor);
1902 parallel_for(vector_multiply, 0, size, thread_loop_partitioner);
1903 }
1904
1905 static void
1906 scale(const std::shared_ptr<::parallel::internal::TBBPartitioner>
1907 &thread_loop_partitioner,
1908 const size_type size,
1909 const ::MemorySpace::
1910 MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1911 ::MemorySpace::MemorySpaceData<Number,
1912 ::MemorySpace::Host>
1913 &data)
1914 {
1916 v_data.values.data());
1917 parallel_for(vector_scale, 0, size, thread_loop_partitioner);
1918 }
1919
1920 static void
1921 equ_au(const std::shared_ptr<::parallel::internal::TBBPartitioner>
1922 &thread_loop_partitioner,
1923 const size_type size,
1924 const Number a,
1925 const ::MemorySpace::
1926 MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1927 ::MemorySpace::MemorySpaceData<Number,
1928 ::MemorySpace::Host>
1929 &data)
1930 {
1932 v_data.values.data(),
1933 a);
1934 parallel_for(vector_equ, 0, size, thread_loop_partitioner);
1935 }
1936
1937 static void
1939 const std::shared_ptr<::parallel::internal::TBBPartitioner>
1940 &thread_loop_partitioner,
1941 const size_type size,
1942 const Number a,
1943 const Number b,
1944 const ::MemorySpace::
1945 MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1946 const ::MemorySpace::
1947 MemorySpaceData<Number, ::MemorySpace::Host> &w_data,
1948 ::MemorySpace::MemorySpaceData<Number,
1949 ::MemorySpace::Host>
1950 &data)
1951 {
1953 data.values.data(), v_data.values.data(), w_data.values.data(), a, b);
1954 parallel_for(vector_equ, 0, size, thread_loop_partitioner);
1955 }
1956
1957 static Number
1958 dot(const std::shared_ptr<::parallel::internal::TBBPartitioner>
1959 &thread_loop_partitioner,
1960 const size_type size,
1961 const ::MemorySpace::
1962 MemorySpaceData<Number2, ::MemorySpace::Host> &v_data,
1963 ::MemorySpace::MemorySpaceData<Number,
1964 ::MemorySpace::Host>
1965 &data)
1966 {
1967 Number sum;
1969 data.values.data(), v_data.values.data());
1971 dot, 0, size, sum, thread_loop_partitioner);
1972 AssertIsFinite(sum);
1973
1974 return sum;
1975 }
1976
1977 template <typename real_type>
1978 static void
1979 norm_2(const std::shared_ptr<::parallel::internal::TBBPartitioner>
1980 &thread_loop_partitioner,
1981 const size_type size,
1982 real_type &sum,
1985 &data)
1986 {
1987 Norm2<Number, real_type> norm2(data.values.data());
1988 parallel_reduce(norm2, 0, size, sum, thread_loop_partitioner);
1989 }
1990
1991 static Number
1993 const std::shared_ptr<::parallel::internal::TBBPartitioner>
1994 &thread_loop_partitioner,
1995 const size_type size,
1996 const ::MemorySpace::
1997 MemorySpaceData<Number, ::MemorySpace::Host> &data)
1998 {
1999 Number sum;
2000 MeanValue<Number> mean(data.values.data());
2001 parallel_reduce(mean, 0, size, sum, thread_loop_partitioner);
2002
2003 return sum;
2004 }
2005
2006 template <typename real_type>
2007 static void
2008 norm_1(const std::shared_ptr<::parallel::internal::TBBPartitioner>
2009 &thread_loop_partitioner,
2010 const size_type size,
2011 real_type &sum,
2014 &data,
2015 const size_type optional_offset = 0)
2016 {
2017 Norm1<Number, real_type> norm1(data.values.data());
2021 sum,
2022 thread_loop_partitioner);
2023 }
2024
2025 template <typename real_type>
2026 static void
2027 norm_p(const std::shared_ptr<::parallel::internal::TBBPartitioner>
2028 &thread_loop_partitioner,
2029 const size_type size,
2030 real_type &sum,
2031 const real_type p,
2034 &data)
2035 {
2036 NormP<Number, real_type> normp(data.values.data(), p);
2037 parallel_reduce(normp, 0, size, sum, thread_loop_partitioner);
2038 }
2039
2040 static Number
2042 const std::shared_ptr<::parallel::internal::TBBPartitioner>
2043 &thread_loop_partitioner,
2044 const size_type size,
2045 const Number a,
2046 const ::MemorySpace::
2047 MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
2048 const ::MemorySpace::
2049 MemorySpaceData<Number, ::MemorySpace::Host> &w_data,
2050 ::MemorySpace::MemorySpaceData<Number,
2051 ::MemorySpace::Host>
2052 &data)
2053 {
2054 Number sum;
2055 AddAndDot<Number> adder(data.values.data(),
2056 v_data.values.data(),
2057 w_data.values.data(),
2058 a);
2059 parallel_reduce(adder, 0, size, sum, thread_loop_partitioner);
2060
2061 return sum;
2062 }
2063
2064 template <typename MemorySpace2>
2065 static void
2067 const std::shared_ptr<::parallel::internal::TBBPartitioner>
2068 &thread_loop_partitioner,
2069 const size_type size,
2071 const ::MemorySpace::MemorySpaceData<Number, MemorySpace2>
2072 &v_data,
2075 &data,
2076 std::enable_if_t<
2077 std::is_same_v<MemorySpace2, ::MemorySpace::Host>,
2078 int> = 0)
2079 {
2081 {
2082 copy(thread_loop_partitioner, size, v_data, data);
2083 }
2084 else if (operation == VectorOperation::add)
2085 {
2086 add_vector(thread_loop_partitioner, size, v_data, data);
2087 }
2088 else
2089 {
2091 }
2092 }
2093
2094 template <typename MemorySpace2>
2095 static void
2097 const std::shared_ptr<::parallel::internal::TBBPartitioner>
2098 & /*thread_loop_partitioner*/,
2099 const size_type size,
2101 const ::MemorySpace::MemorySpaceData<Number, MemorySpace2>
2102 &v_data,
2105 &data,
2106 std::enable_if_t<
2107 std::is_same_v<MemorySpace2, ::MemorySpace::Default>,
2108 int> = 0)
2109 {
2111 {
2112 Kokkos::deep_copy(
2113 Kokkos::subview(data.values,
2114 Kokkos::pair<size_type, size_type>(0, size)),
2115 Kokkos::subview(v_data.values,
2116 Kokkos::pair<size_type, size_type>(0, size)));
2117 }
2118 else
2119 {
2121 }
2122 }
2123 };
2124
2125
2126
2127 template <typename Number>
2128 struct functions<Number, Number, ::MemorySpace::Default>
2129 {
2130 static void
2132 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2133 const size_type size,
2134 const ::MemorySpace::
2135 MemorySpaceData<Number, ::MemorySpace::Default> &v_data,
2136 ::MemorySpace::MemorySpaceData<Number,
2137 ::MemorySpace::Default>
2138 &data)
2139 {
2140 Kokkos::deep_copy(
2141 Kokkos::subview(data.values,
2142 Kokkos::pair<size_type, size_type>(0, size)),
2143 Kokkos::subview(v_data.values,
2144 Kokkos::pair<size_type, size_type>(0, size)));
2145 }
2146
2147 static void
2148 set(const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2149 const size_type size,
2150 const Number s,
2153 &data)
2154 {
2155 Kokkos::deep_copy(
2156 Kokkos::subview(data.values,
2157 Kokkos::pair<size_type, size_type>(0, size)),
2158 s);
2159 }
2160
2161 static void
2163 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2164 const size_type size,
2165 const ::MemorySpace::
2166 MemorySpaceData<Number, ::MemorySpace::Default> &v_data,
2167 ::MemorySpace::MemorySpaceData<Number,
2168 ::MemorySpace::Default>
2169 &data)
2170 {
2171 auto exec = typename ::MemorySpace::Default::kokkos_space::
2172 execution_space{};
2173 Kokkos::parallel_for(
2174 "::add_vector",
2175 Kokkos::RangePolicy<
2176 ::MemorySpace::Default::kokkos_space::execution_space>(
2177 exec, 0, size),
2178 KOKKOS_LAMBDA(int i) { data.values(i) += v_data.values(i); });
2179 exec.fence();
2180 }
2181
2182 static void
2184 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2185 const size_type size,
2186 const ::MemorySpace::
2187 MemorySpaceData<Number, ::MemorySpace::Default> &v_data,
2188 ::MemorySpace::MemorySpaceData<Number,
2189 ::MemorySpace::Default>
2190 &data)
2191 {
2192 auto exec = typename ::MemorySpace::Default::kokkos_space::
2193 execution_space{};
2194 Kokkos::parallel_for(
2195 "::subtract_vector",
2196 Kokkos::RangePolicy<
2197 ::MemorySpace::Default::kokkos_space::execution_space>(
2198 exec, 0, size),
2199 KOKKOS_LAMBDA(size_type i) { data.values(i) -= v_data.values(i); });
2200 exec.fence();
2201 }
2202
2203 static void
2205 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2206 const size_type size,
2207 Number a,
2210 &data)
2211 {
2212 auto exec = typename ::MemorySpace::Default::kokkos_space::
2213 execution_space{};
2214 Kokkos::parallel_for(
2215 "::add_factor",
2216 Kokkos::RangePolicy<
2217 ::MemorySpace::Default::kokkos_space::execution_space>(
2218 exec, 0, size),
2219 KOKKOS_LAMBDA(size_type i) { data.values(i) += a; });
2220 exec.fence();
2221 }
2222
2223 static void
2225 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2226 const size_type size,
2227 const Number a,
2228 const ::MemorySpace::
2229 MemorySpaceData<Number, ::MemorySpace::Default> &v_data,
2230 ::MemorySpace::MemorySpaceData<Number,
2231 ::MemorySpace::Default>
2232 &data)
2233 {
2234 auto exec = typename ::MemorySpace::Default::kokkos_space::
2235 execution_space{};
2236 Kokkos::parallel_for(
2237 "::add_av",
2238 Kokkos::RangePolicy<
2239 ::MemorySpace::Default::kokkos_space::execution_space>(
2240 exec, 0, size),
2242 data.values(i) += a * v_data.values(i);
2243 });
2244 exec.fence();
2245 }
2246
2247 static void
2249 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2250 const size_type size,
2251 const Number a,
2252 const Number b,
2253 const ::MemorySpace::
2254 MemorySpaceData<Number, ::MemorySpace::Default> &v_data,
2255 const ::MemorySpace::
2256 MemorySpaceData<Number, ::MemorySpace::Default> &w_data,
2257 ::MemorySpace::MemorySpaceData<Number,
2258 ::MemorySpace::Default>
2259 &data)
2260 {
2261 auto exec = typename ::MemorySpace::Default::kokkos_space::
2262 execution_space{};
2263 Kokkos::parallel_for(
2264 "::add_avpbw",
2265 Kokkos::RangePolicy<
2266 ::MemorySpace::Default::kokkos_space::execution_space>(
2267 exec, 0, size),
2269 data.values(i) += a * v_data.values(i) + b * w_data.values(i);
2270 });
2271 exec.fence();
2272 }
2273
2274 static void
2276 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2277 const size_type size,
2278 const Number x,
2279 const ::MemorySpace::
2280 MemorySpaceData<Number, ::MemorySpace::Default> &v_data,
2281 ::MemorySpace::MemorySpaceData<Number,
2282 ::MemorySpace::Default>
2283 &data)
2284 {
2285 auto exec = typename ::MemorySpace::Default::kokkos_space::
2286 execution_space{};
2287 Kokkos::parallel_for(
2288 "::sadd_xv",
2289 Kokkos::RangePolicy<
2290 ::MemorySpace::Default::kokkos_space::execution_space>(
2291 exec, 0, size),
2293 data.values(i) = x * data.values(i) + v_data.values(i);
2294 });
2295 exec.fence();
2296 }
2297
2298 static void
2300 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2301 const size_type size,
2302 const Number x,
2303 const Number a,
2304 const ::MemorySpace::
2305 MemorySpaceData<Number, ::MemorySpace::Default> &v_data,
2306 ::MemorySpace::MemorySpaceData<Number,
2307 ::MemorySpace::Default>
2308 &data)
2309 {
2310 auto exec = typename ::MemorySpace::Default::kokkos_space::
2311 execution_space{};
2312 Kokkos::parallel_for(
2313 "::sadd_xav",
2314 Kokkos::RangePolicy<
2315 ::MemorySpace::Default::kokkos_space::execution_space>(
2316 exec, 0, size),
2318 data.values(i) = x * data.values(i) + a * v_data.values(i);
2319 });
2320 exec.fence();
2321 }
2322
2323 static void
2325 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2326 const size_type size,
2327 const Number x,
2328 const Number a,
2329 const Number b,
2330 const ::MemorySpace::
2331 MemorySpaceData<Number, ::MemorySpace::Default> &v_data,
2332 const ::MemorySpace::
2333 MemorySpaceData<Number, ::MemorySpace::Default> &w_data,
2334 ::MemorySpace::MemorySpaceData<Number,
2335 ::MemorySpace::Default>
2336 &data)
2337 {
2338 auto exec = typename ::MemorySpace::Default::kokkos_space::
2339 execution_space{};
2340 Kokkos::parallel_for(
2341 "::sadd_xavbw",
2342 Kokkos::RangePolicy<
2343 ::MemorySpace::Default::kokkos_space::execution_space>(
2344 exec, 0, size),
2346 data.values(i) =
2347 x * data.values(i) + a * v_data.values(i) + b * w_data.values(i);
2348 });
2349 exec.fence();
2350 }
2351
2352 static void
2354 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2355 const size_type size,
2356 const Number factor,
2359 &data)
2360 {
2361 auto exec = typename ::MemorySpace::Default::kokkos_space::
2362 execution_space{};
2363 Kokkos::parallel_for(
2364 "::multiply_factor",
2365 Kokkos::RangePolicy<
2366 ::MemorySpace::Default::kokkos_space::execution_space>(
2367 exec, 0, size),
2368 KOKKOS_LAMBDA(size_type i) { data.values(i) *= factor; });
2369 exec.fence();
2370 }
2371
2372 static void
2374 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2375 const size_type size,
2376 const ::MemorySpace::
2377 MemorySpaceData<Number, ::MemorySpace::Default> &v_data,
2378 ::MemorySpace::MemorySpaceData<Number,
2379 ::MemorySpace::Default>
2380 &data)
2381 {
2382 auto exec = typename ::MemorySpace::Default::kokkos_space::
2383 execution_space{};
2384 Kokkos::parallel_for(
2385 "::scale",
2386 Kokkos::RangePolicy<
2387 ::MemorySpace::Default::kokkos_space::execution_space>(
2388 exec, 0, size),
2389 KOKKOS_LAMBDA(size_type i) { data.values(i) *= v_data.values(i); });
2390 exec.fence();
2391 }
2392
2393 static void
2395 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2396 const size_type size,
2397 const Number a,
2398 const ::MemorySpace::
2399 MemorySpaceData<Number, ::MemorySpace::Default> &v_data,
2400 ::MemorySpace::MemorySpaceData<Number,
2401 ::MemorySpace::Default>
2402 &data)
2403 {
2404 auto exec = typename ::MemorySpace::Default::kokkos_space::
2405 execution_space{};
2406 Kokkos::parallel_for(
2407 "::equ_au",
2408 Kokkos::RangePolicy<
2409 ::MemorySpace::Default::kokkos_space::execution_space>(
2410 exec, 0, size),
2412 data.values(i) = a * v_data.values(i);
2413 });
2414 exec.fence();
2415 }
2416
2417 static void
2419 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2420 const size_type size,
2421 const Number a,
2422 const Number b,
2423 const ::MemorySpace::
2424 MemorySpaceData<Number, ::MemorySpace::Default> &v_data,
2425 const ::MemorySpace::
2426 MemorySpaceData<Number, ::MemorySpace::Default> &w_data,
2427 ::MemorySpace::MemorySpaceData<Number,
2428 ::MemorySpace::Default>
2429 &data)
2430 {
2431 auto exec = typename ::MemorySpace::Default::kokkos_space::
2432 execution_space{};
2433 Kokkos::parallel_for(
2434 "::equ_aubv",
2435 Kokkos::RangePolicy<
2436 ::MemorySpace::Default::kokkos_space::execution_space>(
2437 exec, 0, size),
2439 data.values(i) = a * v_data.values(i) + b * w_data.values(i);
2440 });
2441 exec.fence();
2442 }
2443
2444 static Number
2445 dot(const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2446 const size_type size,
2447 const ::MemorySpace::
2448 MemorySpaceData<Number, ::MemorySpace::Default> &v_data,
2449 ::MemorySpace::MemorySpaceData<Number,
2450 ::MemorySpace::Default>
2451 &data)
2452 {
2453 Number result;
2454
2455 auto exec = typename ::MemorySpace::Default::kokkos_space::
2456 execution_space{};
2457 Kokkos::parallel_reduce(
2458 "::dot",
2459 Kokkos::RangePolicy<
2460 ::MemorySpace::Default::kokkos_space::execution_space>(
2461 exec, 0, size),
2462 KOKKOS_LAMBDA(size_type i, Number & update) {
2463 update += data.values(i) * v_data.values(i);
2464 },
2465 result);
2466
2468 return result;
2469 }
2470
2471 template <typename real_type>
2472 static void
2473 norm_2(const std::shared_ptr<::parallel::internal::TBBPartitioner>
2474 &thread_loop_partitioner,
2475 const size_type size,
2476 real_type &sum,
2477 ::MemorySpace::
2478 MemorySpaceData<Number, ::MemorySpace::Default> &data)
2479 {
2480 sum = dot(thread_loop_partitioner, size, data, data);
2481 }
2482
2483 static Number
2485 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2486 const size_type size,
2487 const ::MemorySpace::
2488 MemorySpaceData<Number, ::MemorySpace::Default> &data)
2489 {
2490 Number result;
2491
2492 auto exec = typename ::MemorySpace::Default::kokkos_space::
2493 execution_space{};
2494 Kokkos::parallel_reduce(
2495 "::mean_value",
2496 Kokkos::RangePolicy<
2497 ::MemorySpace::Default::kokkos_space::execution_space>(
2498 exec, 0, size),
2499 KOKKOS_LAMBDA(size_type i, Number & update) {
2500 update += data.values(i);
2501 },
2502 result);
2503
2505 return result;
2506 }
2507
2508 template <typename real_type>
2509 static void
2511 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2512 const size_type size,
2513 real_type &sum,
2516 &data,
2517 const size_type optional_offset = 0)
2518 {
2519 auto exec = typename ::MemorySpace::Default::kokkos_space::
2520 execution_space{};
2521 Kokkos::parallel_reduce(
2522 "::norm_1",
2523 Kokkos::RangePolicy<
2524 ::MemorySpace::Default::kokkos_space::execution_space>(
2526 KOKKOS_LAMBDA(size_type i, Number & update) {
2527#if KOKKOS_VERSION < 30400
2528 update += std::abs(data.values(i));
2529#elif KOKKOS_VERSION < 30700
2530 update += Kokkos::Experimental::fabs(data.values(i));
2531#else
2532 update += Kokkos::abs(data.values(i));
2533#endif
2534 },
2535 sum);
2536 }
2537
2538 template <typename real_type>
2539 static void
2541 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2542 const size_type size,
2543 real_type &sum,
2544 real_type exp,
2547 &data)
2548 {
2549 auto exec = typename ::MemorySpace::Default::kokkos_space::
2550 execution_space{};
2551 Kokkos::parallel_reduce(
2552 "::norm_p",
2553 Kokkos::RangePolicy<
2554 ::MemorySpace::Default::kokkos_space::execution_space>(
2555 exec, 0, size),
2556 KOKKOS_LAMBDA(size_type i, Number & update) {
2557#if KOKKOS_VERSION < 30400
2558 update += std::pow(fabs(data.values(i)), exp);
2559#elif KOKKOS_VERSION < 30700
2560 update += Kokkos::Experimental::pow(
2561 Kokkos::Experimental::fabs(data.values(i)), exp);
2562#else
2563 update += Kokkos::pow(Kokkos::abs(data.values(i)), exp);
2564#endif
2565 },
2566 sum);
2567 }
2568
2569 static Number
2571 const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2572 const size_type size,
2573 const Number a,
2574 const ::MemorySpace::
2575 MemorySpaceData<Number, ::MemorySpace::Default> &v_data,
2576 const ::MemorySpace::
2577 MemorySpaceData<Number, ::MemorySpace::Default> &w_data,
2578 ::MemorySpace::MemorySpaceData<Number,
2579 ::MemorySpace::Default>
2580 &data)
2581 {
2582 Number res;
2583
2584 auto exec = typename ::MemorySpace::Default::kokkos_space::
2585 execution_space{};
2586 Kokkos::parallel_reduce(
2587 "::add_and_dot",
2588 Kokkos::RangePolicy<
2589 ::MemorySpace::Default::kokkos_space::execution_space>(
2590 exec, 0, size),
2591 KOKKOS_LAMBDA(size_type i, Number & update) {
2592 data.values(i) += a * v_data.values(i);
2593 update +=
2595 w_data.values(i)));
2596 },
2597 res);
2598
2599 return res;
2600 }
2601
2602 template <typename MemorySpace2>
2603 static void
2605 const std::shared_ptr<::parallel::internal::TBBPartitioner>
2606 &thread_loop_partitioner,
2607 const size_type size,
2609 const ::MemorySpace::MemorySpaceData<Number, MemorySpace2>
2610 &v_data,
2613 &data,
2614 std::enable_if_t<
2615 std::is_same_v<MemorySpace2, ::MemorySpace::Default>,
2616 int> = 0)
2617 {
2619 {
2620 copy(thread_loop_partitioner, size, v_data, data);
2621 }
2622 else if (operation == VectorOperation::add)
2623 {
2624 add_vector(thread_loop_partitioner, size, v_data, data);
2625 }
2626 else
2627 {
2629 }
2630 }
2631
2632 template <typename MemorySpace2>
2633 static void
2635 const std::shared_ptr<::parallel::internal::TBBPartitioner>
2636 & /*thread_loop_partitioner*/,
2637 const size_type size,
2639 const ::MemorySpace::MemorySpaceData<Number, MemorySpace2>
2640 &v_data,
2643 &data,
2644 std::enable_if_t<
2645 std::is_same_v<MemorySpace2, ::MemorySpace::Host>,
2646 int> = 0)
2647 {
2649 {
2650 Kokkos::deep_copy(
2651 Kokkos::subview(data.values,
2652 Kokkos::pair<size_type, size_type>(0, size)),
2653 Kokkos::subview(v_data.values,
2654 Kokkos::pair<size_type, size_type>(0, size)));
2655 }
2656 else
2657 {
2659 }
2660 }
2661 };
2662 } // namespace VectorOperations
2663} // namespace internal
2664
2666
2667#endif
static unsigned int n_threads()
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
Definition tensor.h:851
#define DEAL_II_ALWAYS_INLINE
Definition config.h:109
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition config.h:156
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:518
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:519
#define DEAL_II_FALLTHROUGH
Definition config.h:215
Point< 2 > first
Definition grid_out.cc:4629
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIsFinite(number)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
std::vector< index_type > data
Definition mpi.cc:740
std::size_t size
Definition mpi.cc:739
unsigned int minimum_parallel_grain_size
Definition parallel.cc:50
void accumulate_recursive(const Operation &op, const size_type first, const size_type last, ResultType &result)
size_type do_accumulate(const Operation op, const size_type vec_size, const size_type start_index, ResultType *outer_results, std::bool_constant< false >)
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)
const unsigned int vector_accumulation_recursion_threshold
void parallel_for(Iterator x_begin, Iterator x_end, const Functor &functor, const unsigned int grainsize)
Definition parallel.h:122
::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:90
AddAndDot(Number *const X, const Number *const V, const Number *const W, const Number a)
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
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)
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 > &, const size_type size, VectorOperation::values operation, const ::MemorySpace::MemorySpaceData< Number, MemorySpace2 > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data, std::enable_if_t< std::is_same_v< MemorySpace2, ::MemorySpace::Default >, int >=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 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 > &thread_loop_partitioner, const size_type size, VectorOperation::values operation, const ::MemorySpace::MemorySpaceData< Number, MemorySpace2 > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Host > &data, std::enable_if_t< std::is_same_v< MemorySpace2, ::MemorySpace::Host >, int >=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 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, const size_type optional_offset=0)
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 subtract_vector(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &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::Default > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &data)
static Number dot(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &data)
static void add_av(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const Number a, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &data)
static void copy(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &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::Default > &v_data, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &w_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &data)
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::Default > &v_data, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &w_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &data)
static Number mean_value(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &data)
static void norm_p(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, real_type &sum, real_type exp, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &data)
static void add_factor(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, Number a, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &data)
static void add_vector(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &data)
static void norm_1(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, real_type &sum, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &data, const size_type optional_offset=0)
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::Default > &v_data, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &w_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &data)
static void set(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const Number s, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &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::Default > &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::Default > &data, std::enable_if_t< std::is_same_v< MemorySpace2, ::MemorySpace::Host >, int >=0)
static void sadd_xv(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const Number x, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &data)
static void multiply_factor(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const Number factor, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &data)
static void scale(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &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::Default > &v_data, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &w_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &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::Default > &data, std::enable_if_t< std::is_same_v< MemorySpace2, ::MemorySpace::Default >, int >=0)
static void equ_au(const std::shared_ptr<::parallel::internal::TBBPartitioner > &, const size_type size, const Number a, const ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &v_data, ::MemorySpace::MemorySpaceData< Number, ::MemorySpace::Default > &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 real_type abs(const number &x)
Definition numbers.h:575
static constexpr real_type abs_square(const number &x)
Definition numbers.h:566