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