Reference documentation for deal.II version GIT 5dd0dfed0b 2022-12-08 21:45:02+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\}}\)
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 
25 #include <deal.II/base/parallel.h>
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 
38 namespace internal
39 {
40  namespace VectorOperations
41  {
43 
44  template <typename T>
45  bool
46  is_non_negative(const T &t)
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)
108  : functor(functor)
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 =
117  std::min(static_cast<size_type>(4 * MultithreadInfo::n_threads()),
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;
142  const size_type end;
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>
201  struct Vector_set
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
211  operator()(const size_type begin, const size_type end) const
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>
235  struct Vector_copy
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
246  operator()(const size_type begin, const size_type end) const
247  {
249 
250 #ifdef DEAL_II_HAVE_CXX17
251  if constexpr (std::is_trivially_copyable<Number>() &&
252  std::is_same<Number, OtherNumber>::value)
253 #else
254  if (std::is_trivially_copyable<Number>() &&
255  std::is_same<Number, OtherNumber>::value)
256 #endif
257  std::memcpy(dst + begin, src + begin, (end - begin) * sizeof(Number));
258  else
259  {
261  for (size_type i = begin; i < end; ++i)
262  dst[i] = src[i];
263  }
264  }
265 
266  const OtherNumber *const src;
267  Number *const dst;
268  };
269 
270  template <typename Number>
272  {
273  Vectorization_multiply_factor(Number *const val, const Number factor)
274  : val(val)
275  , stored_factor(factor)
276  {}
277 
278  void
279  operator()(const size_type begin, const size_type end) const
280  {
281  // create a local copy of the variable to help the compiler with the
282  // aliasing analysis
283  const Number factor = stored_factor;
284 
286  {
288  for (size_type i = begin; i < end; ++i)
289  val[i] *= factor;
290  }
291  else
292  {
293  for (size_type i = begin; i < end; ++i)
294  val[i] *= factor;
295  }
296  }
297 
298  Number *const val;
299  const Number stored_factor;
300  };
301 
302  template <typename Number>
304  {
305  Vectorization_add_av(Number *const val,
306  const Number *const v_val,
307  const Number factor)
308  : val(val)
309  , v_val(v_val)
310  , stored_factor(factor)
311  {}
312 
313  void
314  operator()(const size_type begin, const size_type end) const
315  {
316  // create a local copy of the variable to help the compiler with the
317  // aliasing analysis
318  const Number factor = stored_factor;
320  {
322  for (size_type i = begin; i < end; ++i)
323  val[i] += factor * v_val[i];
324  }
325  else
326  {
327  for (size_type i = begin; i < end; ++i)
328  val[i] += factor * v_val[i];
329  }
330  }
331 
332  Number *const val;
333  const Number *const v_val;
334  const Number stored_factor;
335  };
336 
337  template <typename Number>
339  {
341  const Number *const v_val,
342  const Number a,
343  const Number x)
344  : val(val)
345  , v_val(v_val)
346  , stored_a(a)
347  , stored_x(x)
348  {}
349 
350  void
351  operator()(const size_type begin, const size_type end) const
352  {
353  // create a local copy of the variable to help the compiler with the
354  // aliasing analysis
355  const Number x = stored_x, a = stored_a;
356 
358  {
360  for (size_type i = begin; i < end; ++i)
361  val[i] = x * val[i] + a * v_val[i];
362  }
363  else
364  {
365  for (size_type i = begin; i < end; ++i)
366  val[i] = x * val[i] + a * v_val[i];
367  }
368  }
369 
370  Number *const val;
371  const Number *const v_val;
372  const Number stored_a;
373  const Number stored_x;
374  };
375 
376  template <typename Number>
378  {
379  Vectorization_subtract_v(Number *val, const Number *const v_val)
380  : val(val)
381  , v_val(v_val)
382  {}
383 
384  void
385  operator()(const size_type begin, const size_type end) const
386  {
388  {
390  for (size_type i = begin; i < end; ++i)
391  val[i] -= v_val[i];
392  }
393  else
394  {
395  for (size_type i = begin; i < end; ++i)
396  val[i] -= v_val[i];
397  }
398  }
399 
400  Number *const val;
401  const Number *const v_val;
402  };
403 
404  template <typename Number>
406  {
407  Vectorization_add_factor(Number *const val, const Number factor)
408  : val(val)
409  , stored_factor(factor)
410  {}
411 
412  void
413  operator()(const size_type begin, const size_type end) const
414  {
415  const Number factor = stored_factor;
416 
418  {
420  for (size_type i = begin; i < end; ++i)
421  val[i] += factor;
422  }
423  else
424  {
425  for (size_type i = begin; i < end; ++i)
426  val[i] += factor;
427  }
428  }
429 
430  Number *const val;
431  const Number stored_factor;
432  };
433 
434  template <typename Number>
436  {
437  Vectorization_add_v(Number *const val, const Number *const v_val)
438  : val(val)
439  , v_val(v_val)
440  {}
441 
442  void
443  operator()(const size_type begin, const size_type end) const
444  {
446  {
448  for (size_type i = begin; i < end; ++i)
449  val[i] += v_val[i];
450  }
451  else
452  {
453  for (size_type i = begin; i < end; ++i)
454  val[i] += v_val[i];
455  }
456  }
457 
458  Number *const val;
459  const Number *const v_val;
460  };
461 
462  template <typename Number>
464  {
466  const Number *const v_val,
467  const Number *const w_val,
468  const Number a,
469  const Number b)
470  : val(val)
471  , v_val(v_val)
472  , w_val(w_val)
473  , stored_a(a)
474  , stored_b(b)
475  {}
476 
477  void
478  operator()(const size_type begin, const size_type end) const
479  {
480  const Number a = stored_a, b = stored_b;
481 
483  {
485  for (size_type i = begin; i < end; ++i)
486  val[i] = val[i] + a * v_val[i] + b * w_val[i];
487  }
488  else
489  {
490  for (size_type i = begin; i < end; ++i)
491  val[i] = val[i] + a * v_val[i] + b * w_val[i];
492  }
493  }
494 
495  Number *const val;
496  const Number *const v_val;
497  const Number *const w_val;
498  const Number stored_a;
499  const Number stored_b;
500  };
501 
502  template <typename Number>
504  {
505  Vectorization_sadd_xv(Number *const val,
506  const Number *const v_val,
507  const Number x)
508  : val(val)
509  , v_val(v_val)
510  , stored_x(x)
511  {}
512 
513  void
514  operator()(const size_type begin, const size_type end) const
515  {
516  const Number x = stored_x;
517 
519  {
521  for (size_type i = begin; i < end; ++i)
522  val[i] = x * val[i] + v_val[i];
523  }
524  else
525  {
526  for (size_type i = begin; i < end; ++i)
527  val[i] = x * val[i] + v_val[i];
528  }
529  }
530 
531  Number *const val;
532  const Number *const v_val;
533  const Number stored_x;
534  };
535 
536  template <typename Number>
538  {
540  const Number *v_val,
541  const Number *w_val,
542  Number x,
543  Number a,
544  Number b)
545  : val(val)
546  , v_val(v_val)
547  , w_val(w_val)
548  , stored_x(x)
549  , stored_a(a)
550  , stored_b(b)
551  {}
552 
553  void
554  operator()(const size_type begin, const size_type end) const
555  {
556  const Number x = stored_x, a = stored_a, b = stored_b;
557 
559  {
561  for (size_type i = begin; i < end; ++i)
562  val[i] = x * val[i] + a * v_val[i] + b * w_val[i];
563  }
564  else
565  {
566  for (size_type i = begin; i < end; ++i)
567  val[i] = x * val[i] + a * v_val[i] + b * w_val[i];
568  }
569  }
570 
571  Number *const val;
572  const Number *const v_val;
573  const Number *const w_val;
574  const Number stored_x;
575  const Number stored_a;
576  const Number stored_b;
577  };
578 
579  template <typename Number>
581  {
582  Vectorization_scale(Number *const val, const Number *const v_val)
583  : val(val)
584  , v_val(v_val)
585  {}
586 
587  void
588  operator()(const size_type begin, const size_type end) const
589  {
591  {
593  for (size_type i = begin; i < end; ++i)
594  val[i] *= v_val[i];
595  }
596  else
597  {
598  for (size_type i = begin; i < end; ++i)
599  val[i] *= v_val[i];
600  }
601  }
602 
603  Number *const val;
604  const Number *const v_val;
605  };
606 
607  template <typename Number>
609  {
610  Vectorization_equ_au(Number *const val,
611  const Number *const u_val,
612  const Number a)
613  : val(val)
614  , u_val(u_val)
615  , stored_a(a)
616  {}
617 
618  void
619  operator()(const size_type begin, const size_type end) const
620  {
621  const Number a = stored_a;
622 
624  {
626  for (size_type i = begin; i < end; ++i)
627  val[i] = a * u_val[i];
628  }
629  else
630  {
631  for (size_type i = begin; i < end; ++i)
632  val[i] = a * u_val[i];
633  }
634  }
635 
636  Number *const val;
637  const Number *const u_val;
638  const Number stored_a;
639  };
640 
641  template <typename Number>
643  {
645  const Number *const u_val,
646  const Number *const v_val,
647  const Number a,
648  const Number b)
649  : val(val)
650  , u_val(u_val)
651  , v_val(v_val)
652  , stored_a(a)
653  , stored_b(b)
654  {}
655 
656  void
657  operator()(const size_type begin, const size_type end) const
658  {
659  const Number a = stored_a, b = stored_b;
660 
662  {
664  for (size_type i = begin; i < end; ++i)
665  val[i] = a * u_val[i] + b * v_val[i];
666  }
667  else
668  {
669  for (size_type i = begin; i < end; ++i)
670  val[i] = a * u_val[i] + b * v_val[i];
671  }
672  }
673 
674  Number *const val;
675  const Number *const u_val;
676  const Number *const v_val;
677  const Number stored_a;
678  const Number stored_b;
679  };
680 
681  template <typename Number>
683  {
685  const Number *u_val,
686  const Number *v_val,
687  const Number *w_val,
688  const Number a,
689  const Number b,
690  const Number c)
691  : val(val)
692  , u_val(u_val)
693  , v_val(v_val)
694  , w_val(w_val)
695  , stored_a(a)
696  , stored_b(b)
697  , stored_c(c)
698  {}
699 
700  void
701  operator()(const size_type begin, const size_type end) const
702  {
703  const Number a = stored_a, b = stored_b, c = stored_c;
704 
706  {
708  for (size_type i = begin; i < end; ++i)
709  val[i] = a * u_val[i] + b * v_val[i] + c * w_val[i];
710  }
711  else
712  {
713  for (size_type i = begin; i < end; ++i)
714  val[i] = a * u_val[i] + b * v_val[i] + c * w_val[i];
715  }
716  }
717 
718  Number *const val;
719  const Number *const u_val;
720  const Number *const v_val;
721  const Number *const w_val;
722  const Number stored_a;
723  const Number stored_b;
724  const Number stored_c;
725  };
726 
727  template <typename Number>
729  {
730  Vectorization_ratio(Number *val, const Number *a_val, const Number *b_val)
731  : val(val)
732  , a_val(a_val)
733  , b_val(b_val)
734  {}
735 
736  void
737  operator()(const size_type begin, const size_type end) const
738  {
740  {
742  for (size_type i = begin; i < end; ++i)
743  val[i] = a_val[i] / b_val[i];
744  }
745  else
746  {
747  for (size_type i = begin; i < end; ++i)
748  val[i] = a_val[i] / b_val[i];
749  }
750  }
751 
752  Number *const val;
753  const Number *const a_val;
754  const Number *const b_val;
755  };
756 
757 
758 
759  // All sums over all the vector entries (l2-norm, inner product, etc.) are
760  // performed with the same code, using a templated operation defined
761  // here. There are always two versions defined, a standard one that covers
762  // most cases and a vectorized one which is only for equal types and float
763  // and double.
764  template <typename Number, typename Number2>
765  struct Dot
766  {
767  static constexpr bool vectorizes = std::is_same<Number, Number2>::value &&
769 
770  Dot(const Number *const X, const Number2 *const Y)
771  : X(X)
772  , Y(Y)
773  {}
774 
775  Number
776  operator()(const size_type i) const
777  {
778  return X[i] * Number(numbers::NumberTraits<Number2>::conjugate(Y[i]));
779  }
780 
782  do_vectorized(const size_type i) const
783  {
785  x.load(X + i);
786  y.load(Y + i);
787 
788  // the following operation in VectorizedArray does an element-wise
789  // scalar product without taking into account complex values and
790  // the need to take the complex-conjugate of one argument. this
791  // may be a bug, but because all VectorizedArray classes only
792  // work on real scalars, it doesn't really matter very much.
793  // in any case, assert that we really don't get here for
794  // complex-valued objects
795  static_assert(numbers::NumberTraits<Number>::is_complex == false,
796  "This operation is not correctly implemented for "
797  "complex-valued objects.");
798  return x * y;
799  }
800 
801  const Number *const X;
802  const Number2 *const Y;
803  };
804 
805  template <typename Number, typename RealType>
806  struct Norm2
807  {
808  static const bool vectorizes = VectorizedArray<Number>::size() > 1;
809 
810  Norm2(const Number *const X)
811  : X(X)
812  {}
813 
814  RealType
815  operator()(const size_type i) const
816  {
818  }
819 
821  do_vectorized(const size_type i) const
822  {
824  x.load(X + i);
825  return x * x;
826  }
827 
828  const Number *const X;
829  };
830 
831  template <typename Number, typename RealType>
832  struct Norm1
833  {
834  static const bool vectorizes = VectorizedArray<Number>::size() > 1;
835 
836  Norm1(const Number *X)
837  : X(X)
838  {}
839 
840  RealType
841  operator()(const size_type i) const
842  {
844  }
845 
847  do_vectorized(const size_type i) const
848  {
850  x.load(X + i);
851  return std::abs(x);
852  }
853 
854  const Number *X;
855  };
856 
857  template <typename Number, typename RealType>
858  struct NormP
859  {
860  static const bool vectorizes = VectorizedArray<Number>::size() > 1;
861 
862  NormP(const Number *X, RealType p)
863  : X(X)
864  , p(p)
865  {}
866 
867  RealType
868  operator()(const size_type i) const
869  {
870  return std::pow(numbers::NumberTraits<Number>::abs(X[i]), p);
871  }
872 
874  do_vectorized(const size_type i) const
875  {
877  x.load(X + i);
878  return std::pow(std::abs(x), p);
879  }
880 
881  const Number * X;
882  const RealType p;
883  };
884 
885  template <typename Number>
886  struct MeanValue
887  {
888  static const bool vectorizes = VectorizedArray<Number>::size() > 1;
889 
890  MeanValue(const Number *X)
891  : X(X)
892  {}
893 
894  Number
895  operator()(const size_type i) const
896  {
897  return X[i];
898  }
899 
901  do_vectorized(const size_type i) const
902  {
904  x.load(X + i);
905  return x;
906  }
907 
908  const Number *X;
909  };
910 
911  template <typename Number>
912  struct AddAndDot
913  {
914  static const bool vectorizes = VectorizedArray<Number>::size() > 1;
915 
916  AddAndDot(Number *const X,
917  const Number *const V,
918  const Number *const W,
919  const Number a)
920  : X(X)
921  , V(V)
922  , W(W)
923  , a(a)
924  {}
925 
926  Number
927  operator()(const size_type i) const
928  {
929  X[i] += a * V[i];
930  return X[i] * Number(numbers::NumberTraits<Number>::conjugate(W[i]));
931  }
932 
934  do_vectorized(const size_type i) const
935  {
937  x.load(X + i);
938  v.load(V + i);
939  x += a * v;
940  x.store(X + i);
941  // may only load from W after storing in X because the pointers might
942  // point to the same memory
943  w.load(W + i);
944 
945  // the following operation in VectorizedArray does an element-wise
946  // scalar product without taking into account complex values and
947  // the need to take the complex-conjugate of one argument. this
948  // may be a bug, but because all VectorizedArray classes only
949  // work on real scalars, it doesn't really matter very much.
950  // in any case, assert that we really don't get here for
951  // complex-valued objects
952  static_assert(numbers::NumberTraits<Number>::is_complex == false,
953  "This operation is not correctly implemented for "
954  "complex-valued objects.");
955  return x * w;
956  }
957 
958  Number *const X;
959  const Number *const V;
960  const Number *const W;
961  const Number a;
962  };
963 
964 
965 
966  // this is the main working loop for all vector sums using the templated
967  // operation above. it accumulates the sums using a block-wise summation
968  // algorithm with post-update. this blocked algorithm has been proposed in
969  // a similar form by Castaldo, Whaley and Chronopoulos (SIAM
970  // J. Sci. Comput. 31, 1156-1174, 2008) and we use the smallest possible
971  // block size, 2. Sometimes it is referred to as pairwise summation. The
972  // worst case error made by this algorithm is on the order O(eps *
973  // log2(vec_size)), whereas a naive summation is O(eps * vec_size). Even
974  // though the Kahan summation is even more accurate with an error O(eps)
975  // by carrying along remainders not captured by the main sum, that involves
976  // additional costs which are not worthwhile. See the Wikipedia article on
977  // the Kahan summation algorithm.
978 
979  // The algorithm implemented here has the additional benefit that it is
980  // easily parallelized without changing the order of how the elements are
981  // added (floating point addition is not associative). For the same vector
982  // size and minimum_parallel_grainsize, the blocks are always the
983  // same and added pairwise.
984 
985  // The depth of recursion is controlled by the 'magic' parameter
986  // vector_accumulation_recursion_threshold: If the length is below
987  // vector_accumulation_recursion_threshold * 32 (32 is the part of code we
988  // unroll), a straight loop instead of recursion will be used. At the
989  // innermost level, eight values are added consecutively in order to better
990  // balance multiplications and additions.
991 
992  // Loops are unrolled as follows: the range [first,last) is broken into
993  // @p n_chunks each of size 32 plus the @p remainder.
994  // accumulate_regular() does the work on 32*n_chunks elements employing SIMD
995  // if possible and stores the result of the operation for each chunk in @p outer_results.
996 
997  // The code returns the result as the last argument in order to make
998  // spawning tasks simpler and use automatic template deduction.
999 
1000 
1007 
1008  template <typename Operation, typename ResultType>
1009  void
1010  accumulate_recursive(const Operation &op,
1011  const size_type first,
1012  const size_type last,
1013  ResultType & result)
1014  {
1015  if (first == last)
1016  {
1017  result = ResultType();
1018  return;
1019  }
1020 
1021  const size_type vec_size = last - first;
1022  if (vec_size <= vector_accumulation_recursion_threshold * 32)
1023  {
1024  // The vector is short enough so we perform the summation. We store
1025  // the number of chunks (each 32 indices) for the given vector
1026  // length; all results are stored in outer_results[0,n_chunks). We
1027  // keep twice the number around to be able to do the pairwise
1028  // summation with a single for loop (see the loop over j below)
1029  ResultType outer_results[vector_accumulation_recursion_threshold * 2];
1030 
1031  // Select between the regular version and vectorized version based
1032  // on the number types we are given. To choose the vectorized
1033  // version often enough, we need to have all tasks but the last one
1034  // to be divisible by the vectorization length
1035  size_type n_chunks = do_accumulate(
1036  op,
1037  vec_size,
1038  first,
1039  outer_results,
1040  std::integral_constant<bool, Operation::vectorizes>());
1041 
1042  AssertIndexRange(n_chunks,
1044 
1045  // now sum the results from the chunks stored in
1046  // outer_results[0,n_chunks) recursively
1047  unsigned int j = 0;
1048  constexpr unsigned int n_lanes = VectorizedArray<ResultType>::size();
1049  for (; j + 2 * n_lanes - 1 < n_chunks;
1050  j += 2 * n_lanes, n_chunks += n_lanes)
1051  {
1053  a.load(outer_results + j);
1054  b.load(outer_results + j + n_lanes);
1055  a += b;
1056  a.store(outer_results + n_chunks);
1057  }
1058 
1059  // In the vectorized case, we know the loop bounds and can do things
1060  // more efficiently
1061  if (Operation::vectorizes)
1062  {
1063  AssertDimension(j + n_lanes, n_chunks);
1064  AssertIndexRange(n_chunks,
1066  ResultType *result_ptr = outer_results + j;
1067  if (n_lanes >= 16)
1068  for (unsigned int i = 0; i < 8; ++i)
1069  result_ptr[i] = result_ptr[i] + result_ptr[i + 8];
1070  if (n_lanes >= 8)
1071  for (unsigned int i = 0; i < 4; ++i)
1072  result_ptr[i] = result_ptr[i] + result_ptr[i + 4];
1073  if (n_lanes >= 4)
1074  for (unsigned int i = 0; i < 2; ++i)
1075  result_ptr[i] = result_ptr[i] + result_ptr[i + 2];
1076  result = result_ptr[0] + result_ptr[1];
1077  }
1078  else
1079  {
1080  // Without vectorization, we do not know the exact bounds, so we
1081  // need to continue the variable-length pairwise summation loop
1082  // from above
1083  for (; j + 1 < n_chunks; j += 2, ++n_chunks)
1084  outer_results[n_chunks] =
1085  outer_results[j] + outer_results[j + 1];
1086 
1087  AssertIndexRange(n_chunks,
1089  Assert(n_chunks > 0, ExcInternalError());
1090  result = outer_results[n_chunks - 1];
1091  }
1092  }
1093  else
1094  {
1095  // split vector into four pieces and work on the pieces
1096  // recursively. Make pieces (except last) divisible by one fourth the
1097  // recursion threshold.
1098  const size_type new_size =
1099  (vec_size / (vector_accumulation_recursion_threshold * 32)) *
1101  Assert(first + 3 * new_size < last, ExcInternalError());
1102  ResultType r0, r1, r2, r3;
1103  accumulate_recursive(op, first, first + new_size, r0);
1104  accumulate_recursive(op, first + new_size, first + 2 * new_size, r1);
1106  first + 2 * new_size,
1107  first + 3 * new_size,
1108  r2);
1109  accumulate_recursive(op, first + 3 * new_size, last, r3);
1110  result = (r0 + r1) + (r2 + r3);
1111  }
1112  }
1113 
1114 
1115  // this is the inner working routine for the accumulation loops below. We
1116  // pulled this part out of the regular accumulate routine because we might
1117  // do this thing vectorized (see specialized function below; this is the
1118  // un-vectorized version). As opposed to the vector add functions above,
1119  // we here pass the functor 'op' by value, because we cannot create a copy
1120  // of the scalar inline, and instead make sure that the numbers get local
1121  // (and thus definitely not aliased) for the compiler
1122  template <typename Operation, typename ResultType>
1123  size_type
1124  do_accumulate(const Operation op,
1125  const size_type vec_size,
1126  const size_type start_index,
1127  ResultType * outer_results,
1128  std::integral_constant<bool, false>)
1129  {
1130  // Create local copy to indicate no aliasing to the compiler
1131  size_type index = start_index;
1132 
1133  // choose each chunk to have a width of 32, thereby the index
1134  // is incremented by 4*8 for each @p i.
1135  size_type n_chunks = vec_size / 32;
1136  for (size_type i = 0; i < n_chunks; ++i)
1137  {
1138  ResultType r = {};
1139  for (unsigned int k = 0; k < 2; ++k)
1140  {
1141  ResultType r0 = op(index);
1142  ResultType r1 = op(index + 1);
1143  ResultType r2 = op(index + 2);
1144  ResultType r3 = op(index + 3);
1145  index += 4;
1146  for (size_type j = 1; j < 4; ++j, index += 4)
1147  {
1148  r0 += op(index);
1149  r1 += op(index + 1);
1150  r2 += op(index + 2);
1151  r3 += op(index + 3);
1152  }
1153  r += (r0 + r1) + (r2 + r3);
1154  }
1155  outer_results[i] = r;
1156  }
1157 
1158  if (n_chunks * 32 < vec_size)
1159  {
1160  const size_type remainder = vec_size - n_chunks * 32;
1161  const size_type inner_chunks = remainder / 8;
1162  const size_type remainder_inner = remainder % 8;
1163  ResultType r0 = ResultType(), r1 = ResultType(), r2 = ResultType();
1164  switch (inner_chunks)
1165  {
1166  case 3:
1167  r2 = op(index++);
1168  for (size_type j = 1; j < 8; ++j)
1169  r2 += op(index++);
1171  case 2:
1172  r1 = op(index++);
1173  for (size_type j = 1; j < 8; ++j)
1174  r1 += op(index++);
1175  r1 += r2;
1177  case 1:
1178  r2 = op(index++);
1179  for (size_type j = 1; j < 8; ++j)
1180  r2 += op(index++);
1182  default:
1183  for (size_type j = 0; j < remainder_inner; ++j)
1184  r0 += op(index++);
1185  outer_results[n_chunks++] = (r0 + r2) + r1;
1186  break;
1187  }
1188  }
1189 
1190  // make sure we worked through all indices
1191  AssertDimension(index, start_index + vec_size);
1192 
1193  return n_chunks;
1194  }
1195 
1196 
1197 
1198  // this is the inner working routine for the accumulation loops
1199  // below. This is the specialized case where we can vectorize. We request
1200  // the 'do_vectorized' routine of the operation instead of the regular one
1201  // which does several operations at once. As above, pass in the functor by
1202  // value to create a local copy of the scalar factors in the function (if
1203  // there are any).
1204  template <typename Operation, typename Number>
1205  size_type
1206  do_accumulate(const Operation op,
1207  const size_type vec_size,
1208  const size_type start_index,
1209  Number * outer_results,
1210  std::integral_constant<bool, true>)
1211  {
1212  // Create local copy to indicate no aliasing to the compiler
1213  size_type index = start_index;
1214 
1215  // we start from @p index and workout @p n_chunks each of size 32.
1216  // in order employ SIMD and work on @p nvecs at a time, we split this
1217  // loop yet again:
1218  // First we work on (n_chunks/nvecs) chunks, where each chunk processes
1219  // nvecs*(4*8) elements.
1220 
1221  constexpr size_type n_lanes = VectorizedArray<Number>::size();
1222  const size_type regular_chunks = vec_size / (32 * n_lanes);
1223  for (size_type i = 0; i < regular_chunks; ++i)
1224  {
1225  VectorizedArray<Number> r = {};
1226  for (unsigned int k = 0; k < 2; ++k)
1227  {
1228  VectorizedArray<Number> r0 = op.do_vectorized(index);
1229  VectorizedArray<Number> r1 = op.do_vectorized(index + n_lanes);
1231  op.do_vectorized(index + 2 * n_lanes);
1233  op.do_vectorized(index + 3 * n_lanes);
1234  index += n_lanes * 4;
1235  for (size_type j = 1; j < 4; ++j, index += n_lanes * 4)
1236  {
1237  r0 += op.do_vectorized(index);
1238  r1 += op.do_vectorized(index + n_lanes);
1239  r2 += op.do_vectorized(index + 2 * n_lanes);
1240  r3 += op.do_vectorized(index + 3 * n_lanes);
1241  }
1242  r += (r0 + r1) + (r2 + r3);
1243  }
1244  r.store(&outer_results[i * n_lanes]);
1245  }
1246 
1247  // If we are treating a case where the vector length is not divisible by
1248  // the vectorization length, need a cleanup loop
1249  // The remaining chunks are processed one by one starting from
1250  // regular_chunks * n_lanes; We do as much as possible with 2 SIMD
1251  // operations within each chunk. Here we assume that n_lanes < 32/2 = 16
1252  // as well as 16 % n_lanes == 0.
1253  static_assert(n_lanes <= 16 && 16 % n_lanes == 0,
1254  "VectorizedArray::size() must be 1, 2, 4, 8, or 16");
1255  size_type n_chunks = regular_chunks * n_lanes;
1256  const size_type start_irregular = regular_chunks * n_lanes * 32;
1257  if (start_irregular < vec_size)
1258  {
1260  r1 = VectorizedArray<Number>();
1261  const size_type remainder = vec_size - start_irregular;
1262  const size_type loop_length = remainder / (2 * n_lanes);
1263  for (size_type j = 0; j < loop_length; ++j, index += 2 * n_lanes)
1264  {
1265  r0 += op.do_vectorized(index);
1266  r1 += op.do_vectorized(index + n_lanes);
1267  }
1268  Number scalar_part = Number();
1269  size_type last = remainder % (2 * n_lanes);
1270  if (last > 0)
1271  {
1272  if (last >= n_lanes)
1273  {
1274  r0 += op.do_vectorized(index);
1275  index += n_lanes;
1276  last -= n_lanes;
1277  }
1278  for (unsigned int i = 0; i < last; ++i)
1279  scalar_part += op(index++);
1280  }
1281 
1282  r0 += r1;
1283  r0.store(&outer_results[n_chunks]);
1284  outer_results[n_chunks] += scalar_part;
1285 
1286  // update n_chunks to denote range of entries to sum up in
1287  // outer_results[].
1288  n_chunks += n_lanes;
1289  }
1290 
1291  // make sure we worked through all indices
1292  AssertDimension(index, start_index + vec_size);
1293 
1294  return n_chunks;
1295  }
1296 
1297 
1298 
1299 #ifdef DEAL_II_WITH_TBB
1328  template <typename Operation, typename ResultType>
1330  {
1331  static const unsigned int threshold_array_allocate = 512;
1332 
1333  TBBReduceFunctor(const Operation &op,
1334  const size_type start,
1335  const size_type end)
1336  : op(op)
1337  , start(start)
1338  , end(end)
1339  {
1340  const size_type vec_size = end - start;
1341  // set chunk size for sub-tasks
1342  const unsigned int gs =
1344  n_chunks =
1345  std::min(static_cast<size_type>(4 * MultithreadInfo::n_threads()),
1346  vec_size / gs);
1347  chunk_size = vec_size / n_chunks;
1348 
1349  // round to next multiple of 512 (or leave it at the minimum grain size
1350  // if that happens to be smaller). this is advantageous because our
1351  // algorithm favors lengths of a power of 2 due to pairwise summation ->
1352  // at most one 'oddly' sized chunk
1353  if (chunk_size > 512)
1354  chunk_size = ((chunk_size + 511) / 512) * 512;
1355  n_chunks = (vec_size + chunk_size - 1) / chunk_size;
1356  AssertIndexRange((n_chunks - 1) * chunk_size, vec_size);
1357  AssertIndexRange(vec_size, n_chunks * chunk_size + 1);
1358 
1360  {
1361  // make sure we allocate an even number of elements,
1362  // access to the new last element is needed in do_sum()
1363  large_array.resize(2 * ((n_chunks + 1) / 2));
1364  array_ptr = large_array.data();
1365  }
1366  else
1367  array_ptr = &small_array[0];
1368  }
1369 
1374  void
1375  operator()(const tbb::blocked_range<size_type> &range) const
1376  {
1377  for (size_type i = range.begin(); i < range.end(); ++i)
1379  start + i * chunk_size,
1380  std::min(start + (i + 1) * chunk_size, end),
1381  array_ptr[i]);
1382  }
1383 
1384  ResultType
1385  do_sum() const
1386  {
1387  while (n_chunks > 1)
1388  {
1389  if (n_chunks % 2 == 1)
1390  array_ptr[n_chunks++] = ResultType();
1391  for (size_type i = 0; i < n_chunks; i += 2)
1392  array_ptr[i / 2] = array_ptr[i] + array_ptr[i + 1];
1393  n_chunks /= 2;
1394  }
1395  return array_ptr[0];
1396  }
1397 
1398  const Operation &op;
1401 
1402  mutable unsigned int n_chunks;
1403  unsigned int chunk_size;
1405  std::vector<ResultType> large_array;
1406  // this variable either points to small_array or large_array depending on
1407  // the number of threads we want to feed
1408  mutable ResultType *array_ptr;
1409  };
1410 #endif
1411 
1412 
1413 
1418  template <typename Operation, typename ResultType>
1419 #ifndef DEBUG
1421 #endif
1422  inline void
1424  const Operation &op,
1425  const size_type start,
1426  const size_type end,
1427  ResultType & result,
1428  const std::shared_ptr<::parallel::internal::TBBPartitioner>
1429  &partitioner)
1430  {
1431 #ifdef DEAL_II_WITH_TBB
1432  const size_type vec_size = end - start;
1433  // only go to the parallel function in case there are at least 4 parallel
1434  // items, otherwise the overhead is too large
1435  if (vec_size >=
1438  {
1439  Assert(partitioner.get() != nullptr,
1441  "Unexpected initialization of Vector that does "
1442  "not set the TBB partitioner to a usable state."));
1443  std::shared_ptr<tbb::affinity_partitioner> tbb_partitioner =
1444  partitioner->acquire_one_partitioner();
1445 
1446  TBBReduceFunctor<Operation, ResultType> generic_functor(op,
1447  start,
1448  end);
1449  // We use a minimum grain size of 1 here since the grains at this
1450  // stage of dividing the work refer to the number of vector chunks
1451  // that are processed by (possibly different) threads in the
1452  // parallelized for loop (i.e., they do not refer to individual
1453  // vector entries). The number of chunks here is calculated inside
1454  // TBBForFunctor. See also GitHub issue #2496 for further discussion
1455  // of this strategy.
1457  static_cast<size_type>(0),
1458  static_cast<size_type>(generic_functor.n_chunks),
1459  generic_functor,
1460  1,
1461  tbb_partitioner);
1462  partitioner->release_one_partitioner(tbb_partitioner);
1463  result = generic_functor.do_sum();
1464  }
1465  else
1466  accumulate_recursive(op, start, end, result);
1467 #else
1468  accumulate_recursive(op, start, end, result);
1469  (void)partitioner;
1470 #endif
1471  }
1472 
1473 
1474  template <typename Number, typename Number2, typename MemorySpace>
1475  struct functions
1476  {
1477  static void
1479  const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1480  /*thread_loop_partitioner*/,
1481  const size_type /*size*/,
1482  const ::MemorySpace::MemorySpaceData<Number2, MemorySpace>
1483  & /*v_data*/,
1485  {
1486  static_assert(
1487  std::is_same<MemorySpace, ::MemorySpace::CUDA>::value &&
1488  std::is_same<Number, Number2>::value,
1489  "For the CUDA MemorySpace Number and Number2 should be the same type");
1490  }
1491 
1492  static void
1494  const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1495  /*thread_loop_partitioner*/,
1496  const size_type /*size*/,
1497  const Number /*s*/,
1499  {}
1500 
1501  static void
1503  const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1504  /*thread_loop_partitioner*/,
1505  const size_type /*size*/,
1506  const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1507  & /*v_data*/,
1509  {}
1510 
1511  static void
1513  const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1514  /*thread_loop_partitioner*/,
1515  const size_type /*size*/,
1516  const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1517  & /*v_data*/,
1519  {}
1520 
1521  static void
1523  const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1524  /*thread_loop_partitioner*/,
1525  const size_type /*size*/,
1526  Number /*a*/,
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 ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1537  & /*v_data*/,
1539  {}
1540 
1541  static void
1543  const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1544  /*thread_loop_partitioner*/,
1545  const size_type /*size*/,
1546  const Number /*a*/,
1547  const Number /*b*/,
1548  const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1549  & /*v_data*/,
1550  const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1551  & /*w_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 ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1562  & /*v_data*/,
1564  {}
1565 
1566  static void
1568  const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1569  /*thread_loop_partitioner*/,
1570  const size_type /*size*/,
1571  const Number /*x*/,
1572  const Number /*a*/,
1573  const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1574  & /*v_data*/,
1576  {}
1577 
1578  static void
1580  const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1581  /*thread_loop_partitioner*/,
1582  const size_type /*size*/,
1583  const Number /*x*/,
1584  const Number /*a*/,
1585  const Number /*b*/,
1586  const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1587  & /*v_data*/,
1588  const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1589  & /*w_data*/,
1591  {}
1592 
1593  static void
1595  const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1596  /*thread_loop_partitioner*/,
1597  const size_type /*size*/,
1598  const Number /*factor*/,
1600  {}
1601 
1602  static void
1604  const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1605  /*thread_loop_partitioner*/,
1606  const size_type /*size*/,
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 ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1619  & /*v_data*/,
1621  {}
1622 
1623  static void
1625  const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1626  /*thread_loop_partitioner*/,
1627  const size_type /*size*/,
1628  const Number /*a*/,
1629  const Number /*b*/,
1630  const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1631  & /*v_data*/,
1632  const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1633  & /*w_data*/,
1635  {}
1636 
1637  static Number
1639  const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1640  /*thread_loop_partitioner*/,
1641  const size_type /*size*/,
1642  const ::MemorySpace::MemorySpaceData<Number2, MemorySpace>
1643  & /*v_data*/,
1645  {
1646  return Number();
1647  }
1648 
1649  template <typename real_type>
1650  static void
1652  const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1653  /*thread_loop_partitioner*/,
1654  const size_type /*size*/,
1655  real_type & /*sum*/,
1656  const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1657  & /*v_data*/,
1659  {}
1660 
1661  static Number
1663  const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1664  /*thread_loop_partitioner*/,
1665  const size_type /*size*/,
1666  const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1667  & /*data*/)
1668  {
1669  return Number();
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  Number * /*values*/,
1680  Number * /*values_dev*/)
1681  {}
1682 
1683  template <typename real_type>
1684  static void
1686  const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1687  /*thread_loop_partitioner*/,
1688  const size_type /*size*/,
1689  real_type & /*sum*/,
1690  real_type /*p*/,
1692  {}
1693 
1694  static Number
1696  const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1697  /*thread_loop_partitioner*/,
1698  const size_type /*size*/,
1699  const Number /*a*/,
1700  const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1701  & /*v_data*/,
1702  const ::MemorySpace::MemorySpaceData<Number, MemorySpace>
1703  & /*w_data*/,
1705  {
1706  return Number();
1707  }
1708 
1709  template <typename MemorySpace2>
1710  static void
1712  const std::shared_ptr<::parallel::internal::TBBPartitioner> &
1713  /*thread_loop_partitioner*/,
1714  const size_type /*size*/,
1715  VectorOperation::values /*operation*/,
1716  const ::MemorySpace::MemorySpaceData<Number, MemorySpace2>
1717  & /*v_data*/,
1719  {}
1720  };
1721 
1722 
1723 
1724  template <typename Number, typename Number2>
1725  struct functions<Number, Number2, ::MemorySpace::Host>
1726  {
1727  static void
1728  copy(const std::shared_ptr<::parallel::internal::TBBPartitioner>
1729  & thread_loop_partitioner,
1730  const size_type size,
1731  const ::MemorySpace::
1732  MemorySpaceData<Number2, ::MemorySpace::Host> &v_data,
1735  &data)
1736  {
1737  Vector_copy<Number, Number2> copier(v_data.values.get(),
1738  data.values.get());
1739  parallel_for(copier, 0, size, thread_loop_partitioner);
1740  }
1741 
1742  static void
1743  set(const std::shared_ptr<::parallel::internal::TBBPartitioner>
1744  & thread_loop_partitioner,
1745  const size_type size,
1746  const Number s,
1749  &data)
1750  {
1751  Vector_set<Number> setter(s, data.values.get());
1752  parallel_for(setter, 0, size, thread_loop_partitioner);
1753  }
1754 
1755  static void
1757  const std::shared_ptr<::parallel::internal::TBBPartitioner>
1758  & thread_loop_partitioner,
1759  const size_type size,
1760  const ::MemorySpace::
1761  MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1764  &data)
1765  {
1766  Vectorization_add_v<Number> vector_add(data.values.get(),
1767  v_data.values.get());
1768  parallel_for(vector_add, 0, size, thread_loop_partitioner);
1769  }
1770 
1771  static void
1773  const std::shared_ptr<::parallel::internal::TBBPartitioner>
1774  & thread_loop_partitioner,
1775  const size_type size,
1776  const ::MemorySpace::
1777  MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1780  &data)
1781  {
1782  Vectorization_subtract_v<Number> vector_subtract(data.values.get(),
1783  v_data.values.get());
1784  parallel_for(vector_subtract, 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  Number a,
1795  &data)
1796  {
1797  Vectorization_add_factor<Number> vector_add(data.values.get(), a);
1798  parallel_for(vector_add, 0, size, thread_loop_partitioner);
1799  }
1800 
1801  static void
1802  add_av(const std::shared_ptr<::parallel::internal::TBBPartitioner>
1803  & thread_loop_partitioner,
1804  const size_type size,
1805  const Number a,
1806  const ::MemorySpace::
1807  MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1810  &data)
1811  {
1812  Vectorization_add_av<Number> vector_add(data.values.get(),
1813  v_data.values.get(),
1814  a);
1815  parallel_for(vector_add, 0, size, thread_loop_partitioner);
1816  }
1817 
1818  static void
1820  const std::shared_ptr<::parallel::internal::TBBPartitioner>
1821  & thread_loop_partitioner,
1822  const size_type size,
1823  const Number a,
1824  const Number b,
1825  const ::MemorySpace::
1826  MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1827  const ::MemorySpace::
1828  MemorySpaceData<Number, ::MemorySpace::Host> &w_data,
1831  &data)
1832  {
1834  data.values.get(), v_data.values.get(), w_data.values.get(), a, b);
1835  parallel_for(vector_add, 0, size, thread_loop_partitioner);
1836  }
1837 
1838  static void
1840  const std::shared_ptr<::parallel::internal::TBBPartitioner>
1841  & thread_loop_partitioner,
1842  const size_type size,
1843  const Number x,
1844  const ::MemorySpace::
1845  MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1848  &data)
1849  {
1850  Vectorization_sadd_xv<Number> vector_sadd(data.values.get(),
1851  v_data.values.get(),
1852  x);
1853  parallel_for(vector_sadd, 0, size, thread_loop_partitioner);
1854  }
1855 
1856  static void
1858  const std::shared_ptr<::parallel::internal::TBBPartitioner>
1859  & thread_loop_partitioner,
1860  const size_type size,
1861  const Number x,
1862  const Number a,
1863  const ::MemorySpace::
1864  MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1867  &data)
1868  {
1869  Vectorization_sadd_xav<Number> vector_sadd(data.values.get(),
1870  v_data.values.get(),
1871  a,
1872  x);
1873  parallel_for(vector_sadd, 0, size, thread_loop_partitioner);
1874  }
1875 
1876  static void
1878  const std::shared_ptr<::parallel::internal::TBBPartitioner>
1879  & thread_loop_partitioner,
1880  const size_type size,
1881  const Number x,
1882  const Number a,
1883  const Number b,
1884  const ::MemorySpace::
1885  MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1886  const ::MemorySpace::
1887  MemorySpaceData<Number, ::MemorySpace::Host> &w_data,
1890  &data)
1891  {
1893  data.values.get(), v_data.values.get(), w_data.values.get(), x, a, b);
1894  parallel_for(vector_sadd, 0, size, thread_loop_partitioner);
1895  }
1896 
1897  static void
1899  const std::shared_ptr<::parallel::internal::TBBPartitioner>
1900  & thread_loop_partitioner,
1901  const size_type size,
1902  const Number factor,
1905  &data)
1906  {
1907  Vectorization_multiply_factor<Number> vector_multiply(data.values.get(),
1908  factor);
1909  parallel_for(vector_multiply, 0, size, thread_loop_partitioner);
1910  }
1911 
1912  static void
1913  scale(const std::shared_ptr<::parallel::internal::TBBPartitioner>
1914  & thread_loop_partitioner,
1915  const size_type size,
1916  const ::MemorySpace::
1917  MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1920  &data)
1921  {
1922  Vectorization_scale<Number> vector_scale(data.values.get(),
1923  v_data.values.get());
1924  parallel_for(vector_scale, 0, size, thread_loop_partitioner);
1925  }
1926 
1927  static void
1928  equ_au(const std::shared_ptr<::parallel::internal::TBBPartitioner>
1929  & thread_loop_partitioner,
1930  const size_type size,
1931  const Number a,
1932  const ::MemorySpace::
1933  MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1936  &data)
1937  {
1938  Vectorization_equ_au<Number> vector_equ(data.values.get(),
1939  v_data.values.get(),
1940  a);
1941  parallel_for(vector_equ, 0, size, thread_loop_partitioner);
1942  }
1943 
1944  static void
1946  const std::shared_ptr<::parallel::internal::TBBPartitioner>
1947  & thread_loop_partitioner,
1948  const size_type size,
1949  const Number a,
1950  const Number b,
1951  const ::MemorySpace::
1952  MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
1953  const ::MemorySpace::
1954  MemorySpaceData<Number, ::MemorySpace::Host> &w_data,
1957  &data)
1958  {
1959  Vectorization_equ_aubv<Number> vector_equ(
1960  data.values.get(), v_data.values.get(), w_data.values.get(), a, b);
1961  parallel_for(vector_equ, 0, size, thread_loop_partitioner);
1962  }
1963 
1964  static Number
1965  dot(const std::shared_ptr<::parallel::internal::TBBPartitioner>
1966  & thread_loop_partitioner,
1967  const size_type size,
1968  const ::MemorySpace::
1969  MemorySpaceData<Number2, ::MemorySpace::Host> &v_data,
1972  &data)
1973  {
1974  Number sum;
1976  data.values.get(), v_data.values.get());
1978  dot, 0, size, sum, thread_loop_partitioner);
1980 
1981  return sum;
1982  }
1983 
1984  template <typename real_type>
1985  static void
1986  norm_2(const std::shared_ptr<::parallel::internal::TBBPartitioner>
1987  & thread_loop_partitioner,
1988  const size_type size,
1989  real_type & sum,
1992  &data)
1993  {
1994  Norm2<Number, real_type> norm2(data.values.get());
1995  parallel_reduce(norm2, 0, size, sum, thread_loop_partitioner);
1996  }
1997 
1998  static Number
2000  const std::shared_ptr<::parallel::internal::TBBPartitioner>
2001  & thread_loop_partitioner,
2002  const size_type size,
2003  const ::MemorySpace::
2004  MemorySpaceData<Number, ::MemorySpace::Host> &data)
2005  {
2006  Number sum;
2007  MeanValue<Number> mean(data.values.get());
2008  parallel_reduce(mean, 0, size, sum, thread_loop_partitioner);
2009 
2010  return sum;
2011  }
2012 
2013  template <typename real_type>
2014  static void
2015  norm_1(const std::shared_ptr<::parallel::internal::TBBPartitioner>
2016  & thread_loop_partitioner,
2017  const size_type size,
2018  real_type & sum,
2021  &data)
2022  {
2023  Norm1<Number, real_type> norm1(data.values.get());
2024  parallel_reduce(norm1, 0, size, sum, thread_loop_partitioner);
2025  }
2026 
2027  template <typename real_type>
2028  static void
2029  norm_p(const std::shared_ptr<::parallel::internal::TBBPartitioner>
2030  & thread_loop_partitioner,
2031  const size_type size,
2032  real_type & sum,
2033  const real_type p,
2036  &data)
2037  {
2038  NormP<Number, real_type> normp(data.values.get(), p);
2039  parallel_reduce(normp, 0, size, sum, thread_loop_partitioner);
2040  }
2041 
2042  static Number
2044  const std::shared_ptr<::parallel::internal::TBBPartitioner>
2045  & thread_loop_partitioner,
2046  const size_type size,
2047  const Number a,
2048  const ::MemorySpace::
2049  MemorySpaceData<Number, ::MemorySpace::Host> &v_data,
2050  const ::MemorySpace::
2051  MemorySpaceData<Number, ::MemorySpace::Host> &w_data,
2054  &data)
2055  {
2056  Number sum;
2057  AddAndDot<Number> adder(data.values.get(),
2058  v_data.values.get(),
2059  w_data.values.get(),
2060  a);
2061  parallel_reduce(adder, 0, size, sum, thread_loop_partitioner);
2062 
2063  return sum;
2064  }
2065 
2066  template <typename MemorySpace2>
2067  static void
2069  const std::shared_ptr<::parallel::internal::TBBPartitioner>
2070  & thread_loop_partitioner,
2071  const size_type size,
2072  VectorOperation::values operation,
2073  const ::MemorySpace::MemorySpaceData<Number, MemorySpace2>
2074  &v_data,
2077  &data,
2078  std::enable_if_t<
2079  std::is_same<MemorySpace2, ::MemorySpace::Host>::value,
2080  int> = 0)
2081  {
2082  if (operation == VectorOperation::insert)
2083  {
2084  copy(thread_loop_partitioner, size, v_data, data);
2085  }
2086  else if (operation == VectorOperation::add)
2087  {
2088  add_vector(thread_loop_partitioner, size, v_data, data);
2089  }
2090  else
2091  {
2092  AssertThrow(false, ExcNotImplemented());
2093  }
2094  }
2095 
2096 #ifdef DEAL_II_COMPILER_CUDA_AWARE
2097  template <typename MemorySpace2>
2098  static void
2100  const std::shared_ptr<::parallel::internal::TBBPartitioner>
2101  & /*thread_loop_partitioner*/,
2102  const size_type size,
2103  VectorOperation::values operation,
2104  const ::MemorySpace::MemorySpaceData<Number, MemorySpace2>
2105  &v_data,
2108  &data,
2109  std::enable_if_t<
2110  std::is_same<MemorySpace2, ::MemorySpace::CUDA>::value,
2111  int> = 0)
2112  {
2113  if (operation == VectorOperation::insert)
2114  {
2115  cudaError_t cuda_error_code = cudaMemcpy(data.values.get(),
2116  v_data.values_dev.get(),
2117  size * sizeof(Number),
2118  cudaMemcpyDeviceToHost);
2119  AssertCuda(cuda_error_code);
2120  }
2121  else
2122  {
2123  AssertThrow(false, ExcNotImplemented());
2124  }
2125  }
2126 #endif
2127  };
2128 
2129 
2130 
2131 #ifdef DEAL_II_COMPILER_CUDA_AWARE
2132  template <typename Number>
2133  struct functions<Number, Number, ::MemorySpace::CUDA>
2134  {
2135  static const int block_size =
2137  static const int chunk_size =
2139 
2140  static void
2142  const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2143  const size_type size,
2144  const ::MemorySpace::
2145  MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2148  &data)
2149  {
2150  cudaError_t cuda_error_code = cudaMemcpy(data.values_dev.get(),
2151  v_data.values_dev.get(),
2152  size * sizeof(Number),
2153  cudaMemcpyDeviceToDevice);
2154  AssertCuda(cuda_error_code);
2155  }
2156 
2157  static void
2158  set(const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2159  const size_type size,
2160  const Number s,
2163  &data)
2164  {
2165  const int n_blocks = 1 + size / (chunk_size * block_size);
2166  ::LinearAlgebra::CUDAWrappers::kernel::set<Number>
2167  <<<n_blocks, block_size>>>(data.values_dev.get(), s, size);
2168  AssertCudaKernel();
2169  }
2170 
2171  static void
2173  const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2174  const size_type size,
2175  const ::MemorySpace::
2176  MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2179  &data)
2180  {
2181  const int n_blocks = 1 + size / (chunk_size * block_size);
2182  ::LinearAlgebra::CUDAWrappers::kernel::add_aV<Number>
2183  <<<n_blocks, block_size>>>(data.values_dev.get(),
2184  1.,
2185  v_data.values_dev.get(),
2186  size);
2187  AssertCudaKernel();
2188  }
2189 
2190  static void
2192  const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2193  const size_type size,
2194  const ::MemorySpace::
2195  MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2198  &data)
2199  {
2200  const int n_blocks = 1 + size / (chunk_size * block_size);
2201  ::LinearAlgebra::CUDAWrappers::kernel::add_aV<Number>
2202  <<<n_blocks, block_size>>>(data.values_dev.get(),
2203  -1.,
2204  v_data.values_dev.get(),
2205  size);
2206  AssertCudaKernel();
2207  }
2208 
2209  static void
2211  const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2212  const size_type size,
2213  Number a,
2216  &data)
2217  {
2218  const int n_blocks = 1 + size / (chunk_size * block_size);
2219  ::LinearAlgebra::CUDAWrappers::kernel::vec_add<Number>
2220  <<<n_blocks, block_size>>>(data.values_dev.get(), a, size);
2221  AssertCudaKernel();
2222  }
2223 
2224  static void
2226  const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2227  const size_type size,
2228  const Number a,
2229  const ::MemorySpace::
2230  MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2233  &data)
2234  {
2235  const int n_blocks = 1 + size / (chunk_size * block_size);
2236  ::LinearAlgebra::CUDAWrappers::kernel::add_aV<Number>
2237  <<<n_blocks, block_size>>>(data.values_dev.get(),
2238  a,
2239  v_data.values_dev.get(),
2240  size);
2241  AssertCudaKernel();
2242  }
2243 
2244  static void
2246  const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2247  const size_type size,
2248  const Number a,
2249  const Number b,
2250  const ::MemorySpace::
2251  MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2252  const ::MemorySpace::
2253  MemorySpaceData<Number, ::MemorySpace::CUDA> &w_data,
2256  &data)
2257  {
2258  const int n_blocks = 1 + size / (chunk_size * block_size);
2259  ::LinearAlgebra::CUDAWrappers::kernel::add_aVbW<Number>
2260  <<<dim3(n_blocks, 1), dim3(block_size)>>>(data.values_dev.get(),
2261  a,
2262  v_data.values_dev.get(),
2263  b,
2264  w_data.values_dev.get(),
2265  size);
2266  AssertCudaKernel();
2267  }
2268 
2269  static void
2271  const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2272  const size_type size,
2273  const Number x,
2274  const ::MemorySpace::
2275  MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2278  &data)
2279  {
2280  const int n_blocks = 1 + size / (chunk_size * block_size);
2281  ::LinearAlgebra::CUDAWrappers::kernel::sadd<Number>
2282  <<<dim3(n_blocks, 1), dim3(block_size)>>>(
2283  x, data.values_dev.get(), 1., v_data.values_dev.get(), size);
2284  AssertCudaKernel();
2285  }
2286 
2287  static void
2289  const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2290  const size_type size,
2291  const Number x,
2292  const Number a,
2293  const ::MemorySpace::
2294  MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2297  &data)
2298  {
2299  const int n_blocks = 1 + size / (chunk_size * block_size);
2300  ::LinearAlgebra::CUDAWrappers::kernel::sadd<Number>
2301  <<<dim3(n_blocks, 1), dim3(block_size)>>>(
2302  x, data.values_dev.get(), a, v_data.values_dev.get(), size);
2303  AssertCudaKernel();
2304  }
2305 
2306  static void
2308  const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2309  const size_type size,
2310  const Number x,
2311  const Number a,
2312  const Number b,
2313  const ::MemorySpace::
2314  MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2315  const ::MemorySpace::
2316  MemorySpaceData<Number, ::MemorySpace::CUDA> &w_data,
2319  &data)
2320  {
2321  const int n_blocks = 1 + size / (chunk_size * block_size);
2322  ::LinearAlgebra::CUDAWrappers::kernel::sadd<Number>
2323  <<<dim3(n_blocks, 1), dim3(block_size)>>>(x,
2324  data.values_dev.get(),
2325  a,
2326  v_data.values_dev.get(),
2327  b,
2328  w_data.values_dev.get(),
2329  size);
2330  AssertCudaKernel();
2331  }
2332 
2333  static void
2335  const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2336  const size_type size,
2337  const Number factor,
2340  &data)
2341  {
2342  const int n_blocks = 1 + size / (chunk_size * block_size);
2343  ::LinearAlgebra::CUDAWrappers::kernel::vec_scale<Number>
2344  <<<n_blocks, block_size>>>(data.values_dev.get(), factor, size);
2345  AssertCudaKernel();
2346  }
2347 
2348  static void
2350  const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2351  const size_type size,
2352  const ::MemorySpace::
2353  MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2356  &data)
2357  {
2358  const int n_blocks = 1 + size / (chunk_size * block_size);
2359  ::LinearAlgebra::CUDAWrappers::kernel::scale<Number>
2360  <<<dim3(n_blocks, 1), dim3(block_size)>>>(data.values_dev.get(),
2361  v_data.values_dev.get(),
2362  size);
2363  AssertCudaKernel();
2364  }
2365 
2366  static void
2368  const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2369  const size_type size,
2370  const Number a,
2371  const ::MemorySpace::
2372  MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2375  &data)
2376  {
2377  const int n_blocks = 1 + size / (chunk_size * block_size);
2378  ::LinearAlgebra::CUDAWrappers::kernel::equ<Number>
2379  <<<dim3(n_blocks, 1), dim3(block_size)>>>(data.values_dev.get(),
2380  a,
2381  v_data.values_dev.get(),
2382  size);
2383  AssertCudaKernel();
2384  }
2385 
2386  static void
2388  const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2389  const size_type size,
2390  const Number a,
2391  const Number b,
2392  const ::MemorySpace::
2393  MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2394  const ::MemorySpace::
2395  MemorySpaceData<Number, ::MemorySpace::CUDA> &w_data,
2398  &data)
2399  {
2400  const int n_blocks = 1 + size / (chunk_size * block_size);
2401  ::LinearAlgebra::CUDAWrappers::kernel::equ<Number>
2402  <<<dim3(n_blocks, 1), dim3(block_size)>>>(data.values_dev.get(),
2403  a,
2404  v_data.values_dev.get(),
2405  b,
2406  w_data.values_dev.get(),
2407  size);
2408  AssertCudaKernel();
2409  }
2410 
2411  static Number
2412  dot(const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2413  const size_type size,
2414  const ::MemorySpace::
2415  MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2418  &data)
2419  {
2420  Number * result_device;
2421  cudaError_t error_code = cudaMalloc(&result_device, sizeof(Number));
2422  AssertCuda(error_code);
2423  error_code = cudaMemset(result_device, 0, sizeof(Number));
2424  AssertCuda(error_code);
2425 
2426  const int n_blocks = 1 + size / (chunk_size * block_size);
2428  Number,
2430  <<<dim3(n_blocks, 1), dim3(block_size)>>>(result_device,
2431  data.values_dev.get(),
2432  v_data.values_dev.get(),
2433  static_cast<unsigned int>(
2434  size));
2435  AssertCudaKernel();
2436 
2437  // Copy the result back to the host
2438  Number result;
2439  error_code = cudaMemcpy(&result,
2440  result_device,
2441  sizeof(Number),
2442  cudaMemcpyDeviceToHost);
2443  AssertCuda(error_code);
2444  // Free the memory on the device
2445  error_code = cudaFree(result_device);
2446  AssertCuda(error_code);
2447 
2448  AssertIsFinite(result);
2449 
2450  return result;
2451  }
2452 
2453  template <typename real_type>
2454  static void
2455  norm_2(const std::shared_ptr<::parallel::internal::TBBPartitioner>
2456  & thread_loop_partitioner,
2457  const size_type size,
2458  real_type & sum,
2461  &data)
2462  {
2463  sum = dot(thread_loop_partitioner, size, data, data);
2464  }
2465 
2466  static Number
2468  const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2469  const size_type size,
2470  const ::MemorySpace::
2471  MemorySpaceData<Number, ::MemorySpace::CUDA> &data)
2472  {
2473  Number * result_device;
2474  cudaError_t error_code = cudaMalloc(&result_device, sizeof(Number));
2475  AssertCuda(error_code);
2476  error_code = cudaMemset(result_device, 0, sizeof(Number));
2477 
2478  const int n_blocks = 1 + size / (chunk_size * block_size);
2480  Number,
2482  <<<dim3(n_blocks, 1), dim3(block_size)>>>(result_device,
2483  data.values_dev.get(),
2484  size);
2485 
2486  // Copy the result back to the host
2487  Number result;
2488  error_code = cudaMemcpy(&result,
2489  result_device,
2490  sizeof(Number),
2491  cudaMemcpyDeviceToHost);
2492  AssertCuda(error_code);
2493  // Free the memory on the device
2494  error_code = cudaFree(result_device);
2495  AssertCuda(error_code);
2496 
2497  return result;
2498  }
2499 
2500  template <typename real_type>
2501  static void
2503  const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2504  const size_type size,
2505  real_type & sum,
2508  &data)
2509  {
2510  Number * result_device;
2511  cudaError_t error_code = cudaMalloc(&result_device, sizeof(Number));
2512  AssertCuda(error_code);
2513  error_code = cudaMemset(result_device, 0, sizeof(Number));
2514 
2515  const int n_blocks = 1 + size / (chunk_size * block_size);
2517  Number,
2519  <<<dim3(n_blocks, 1), dim3(block_size)>>>(result_device,
2520  data.values_dev.get(),
2521  size);
2522 
2523  // Copy the result back to the host
2524  error_code = cudaMemcpy(&sum,
2525  result_device,
2526  sizeof(Number),
2527  cudaMemcpyDeviceToHost);
2528  AssertCuda(error_code);
2529  // Free the memory on the device
2530  error_code = cudaFree(result_device);
2531  AssertCuda(error_code);
2532  }
2533 
2534  template <typename real_type>
2535  static void
2537  const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2538  const size_type,
2539  real_type &,
2540  real_type,
2542  ::MemorySpace::CUDA> &)
2543  {
2544  Assert(false, ExcNotImplemented());
2545  }
2546 
2547  static Number
2549  const std::shared_ptr<::parallel::internal::TBBPartitioner> &,
2550  const size_type size,
2551  const Number a,
2552  const ::MemorySpace::
2553  MemorySpaceData<Number, ::MemorySpace::CUDA> &v_data,
2554  const ::MemorySpace::
2555  MemorySpaceData<Number, ::MemorySpace::CUDA> &w_data,
2558  &data)
2559  {
2560  Number * res_d;
2561  cudaError_t error_code = cudaMalloc(&res_d, sizeof(Number));
2562  AssertCuda(error_code);
2563  error_code = cudaMemset(res_d, 0, sizeof(Number));
2564  AssertCuda(error_code);
2565 
2566  const int n_blocks = 1 + size / (chunk_size * block_size);
2567  ::LinearAlgebra::CUDAWrappers::kernel::add_and_dot<Number>
2568  <<<dim3(n_blocks, 1), dim3(block_size)>>>(res_d,
2569  data.values_dev.get(),
2570  v_data.values_dev.get(),
2571  w_data.values_dev.get(),
2572  a,
2573  size);
2574 
2575  Number res;
2576  error_code =
2577  cudaMemcpy(&res, res_d, sizeof(Number), cudaMemcpyDeviceToHost);
2578  AssertCuda(error_code);
2579  error_code = cudaFree(res_d);
2580 
2581  return res;
2582  }
2583 
2584  template <typename MemorySpace2>
2585  static void
2587  const std::shared_ptr<::parallel::internal::TBBPartitioner>
2588  & thread_loop_partitioner,
2589  const size_type size,
2590  VectorOperation::values operation,
2591  const ::MemorySpace::MemorySpaceData<Number, MemorySpace2>
2592  &v_data,
2595  &data,
2596  std::enable_if_t<
2597  std::is_same<MemorySpace2, ::MemorySpace::CUDA>::value,
2598  int> = 0)
2599  {
2600  if (operation == VectorOperation::insert)
2601  {
2602  copy(thread_loop_partitioner, size, v_data, data);
2603  }
2604  else if (operation == VectorOperation::add)
2605  {
2606  add_vector(thread_loop_partitioner, size, v_data, data);
2607  }
2608  else
2609  {
2610  AssertThrow(false, ExcNotImplemented());
2611  }
2612  }
2613 
2614  template <typename MemorySpace2>
2615  static void
2617  const std::shared_ptr<::parallel::internal::TBBPartitioner>
2618  & /*thread_loop_partitioner*/,
2619  const size_type size,
2620  VectorOperation::values operation,
2621  const ::MemorySpace::MemorySpaceData<Number, MemorySpace2>
2622  &v_data,
2625  &data,
2626  std::enable_if_t<
2627  std::is_same<MemorySpace2, ::MemorySpace::Host>::value,
2628  int> = 0)
2629  {
2630  if (operation == VectorOperation::insert)
2631  {
2632  cudaError_t cuda_error_code = cudaMemcpy(data.values_dev.get(),
2633  v_data.values.get(),
2634  size * sizeof(Number),
2635  cudaMemcpyHostToDevice);
2636  AssertCuda(cuda_error_code);
2637  }
2638  else
2639  {
2640  AssertThrow(false, ExcNotImplemented());
2641  }
2642  }
2643  };
2644 #endif
2645  } // namespace VectorOperations
2646 } // namespace internal
2647 
2649 
2650 #endif
static unsigned int n_threads()
void store(OtherNumber *ptr) const
void load(const OtherNumber *ptr)
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:102
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition: config.h:142
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
#define DEAL_II_FALLTHROUGH
Definition: config.h:176
Point< 2 > first
Definition: grid_out.cc:4605
__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)
static ::ExceptionBase & ExcInternalError()
#define AssertCudaKernel()
Definition: exceptions.h:1899
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcNotImplemented()
#define AssertIsFinite(number)
Definition: exceptions.h:1786
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1695
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
#define AssertCuda(error_code)
Definition: exceptions.h:1844
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1611
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_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition: operators.h:50
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 * begin(VectorType &V)
VectorType::value_type * end(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
size_type do_accumulate(const Operation op, const size_type vec_size, const size_type start_index, ResultType *outer_results, std::integral_constant< bool, false >)
unsigned int global_dof_index
Definition: types.h:81
AddAndDot(Number *const X, const Number *const V, const Number *const W, const Number a)
VectorizedArray< Number > do_vectorized(const size_type i) const
Number operator()(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, std::enable_if_t< std::is_same< MemorySpace2, ::MemorySpace::Host >::value, 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 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, std::enable_if_t< std::is_same< MemorySpace2, ::MemorySpace::CUDA >::value, 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 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 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 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, std::enable_if_t< std::is_same< MemorySpace2, ::MemorySpace::Host >::value, int >=0)
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 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 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, std::enable_if_t< std::is_same< MemorySpace2, ::MemorySpace::CUDA >::value, int >=0)
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 real_type abs(const number &x)
Definition: numbers.h:624
static constexpr std::enable_if_t< std::is_same< Dummy, number >::value &&is_cuda_compatible< Dummy >::value, real_type > abs_square(const number &x)