Tpetra parallel linear algebra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Tpetra_Util.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Tpetra: Templated Linear Algebra Services Package
5 // Copyright (2008) Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // ************************************************************************
38 // @HEADER
39 
40 #ifndef TPETRA_UTIL_HPP
41 #define TPETRA_UTIL_HPP
42 
52 #include "Tpetra_ConfigDefs.hpp"
53 #include "Kokkos_DualView.hpp"
54 #include "KokkosCompat_View.hpp"
55 #include "Teuchos_Assert.hpp"
56 #include "Teuchos_CommHelpers.hpp"
57 #include "Teuchos_OrdinalTraits.hpp"
58 #include "Teuchos_TypeNameTraits.hpp"
59 #include "Teuchos_Utils.hpp"
60 #include <numeric>
61 #include <algorithm>
62 #include <iterator>
63 #include <memory>
64 #include <ostream>
65 #include <sstream>
66 
67 #if defined(HAVE_TPETRA_PRINT_EFFICIENCY_WARNINGS)
68 #define TPETRA_EFFICIENCY_WARNING(throw_exception_test,msg) \
93 { \
94  const bool tpetraEfficiencyWarningTest = (throw_exception_test); \
95  if (tpetraEfficiencyWarningTest) { \
96  std::ostringstream errStream; \
97  errStream << Teuchos::typeName(*this) << ":" << std::endl; \
98  errStream << "Efficiency warning: " << #throw_exception_test << std::endl; \
99  errStream << msg; \
100  std::string err = errStream.str(); \
101  if (TPETRA_PRINTS_EFFICIENCY_WARNINGS && tpetraEfficiencyWarningTest) { \
102  std::cerr << err << std::endl; \
103  } \
104  } \
105 }
106 #else
107 #define TPETRA_EFFICIENCY_WARNING(throw_exception_test,msg)
132 #endif
133 
134 // handle an abuse warning, according to HAVE_TPETRA_THROW_ABUSE_WARNINGS and HAVE_TPETRA_PRINT_ABUSE_WARNINGS
135 #if defined(HAVE_TPETRA_THROW_ABUSE_WARNINGS) || defined(HAVE_TPETRA_PRINT_ABUSE_WARNINGS)
136 #define TPETRA_ABUSE_WARNING(throw_exception_test,Exception,msg) \
138 { \
139  std::ostringstream errStream; \
140  errStream << Teuchos::typeName(*this) << msg; \
141  std::string err = errStream.str(); \
142  if (TPETRA_PRINTS_ABUSE_WARNINGS && (throw_exception_test)) { \
143  std::cerr << err << std::endl; \
144  } \
145  TEUCHOS_TEST_FOR_EXCEPTION(TPETRA_THROWS_ABUSE_WARNINGS && (throw_exception_test), Exception, err); \
146 }
147 #else
148 #define TPETRA_ABUSE_WARNING(throw_exception_test,Exception,msg)
150 #endif
151 
181 #define SHARED_TEST_FOR_EXCEPTION(throw_exception_test,Exception,msg,comm) \
182 { \
183  using Teuchos::outArg; \
184  const int lcl_throw_exception = (throw_exception_test) ? Teuchos::rank(comm)+1 : 0; \
185  int gbl_throw; \
186  Teuchos::reduceAll(comm,Teuchos::REDUCE_MAX,lcl_throw_exception,outArg(gbl_throw)); \
187  TEUCHOS_TEST_FOR_EXCEPTION(gbl_throw,Exception, \
188  msg << " Failure on at least one process, including process " << gbl_throw-1 << "." << std::endl); \
189 }
190 
191 namespace Tpetra {
192 
204  template<typename MapType, typename KeyArgType, typename ValueArgType>
205  typename MapType::iterator efficientAddOrUpdate(MapType& m,
206  const KeyArgType & k,
207  const ValueArgType & v)
208  {
209  typename MapType::iterator lb = m.lower_bound(k);
210  if(lb != m.end() && !(m.key_comp()(k, lb->first))) {
211  lb->second = v;
212  return(lb);
213  }
214  else {
215  typedef typename MapType::value_type MVT;
216  return(m.insert(lb, MVT(k, v)));
217  }
218  }
219 
227  namespace SortDetails {
228 
237  template<class IT1>
238  bool isAlreadySorted(const IT1& first, const IT1& last){
239  typedef typename std::iterator_traits<IT1>::difference_type DT;
240  DT myit = Teuchos::OrdinalTraits<DT>::one();
241  const DT sz = last - first;
242  for(;myit < sz; ++myit){
243  if(first[myit] < first[myit-1]){
244  return false;
245  }
246  }
247  return true;
248  }
249 
259  template<class IT>
260  IT getPivot(const IT& first, const IT& last){
261  IT pivot(first+(last-first)/2);
262  if(*first<=*pivot && *(last-1)<=*first) pivot=first;
263  else if(*(last-1)<=*pivot && *first<= *(last-1)) pivot = last-1;
264  return pivot;
265  }
266 
281  template<class IT1, class IT2>
282  IT1 partition2(
283  const IT1& first1,
284  const IT1& last1,
285  const IT2& first2,
286  const IT2& last2,
287  const IT1& pivot)
288  {
289  typename std::iterator_traits<IT1>::value_type piv(*pivot);
290  std::swap(*pivot, *(last1-1));
291  std::swap(first2[(pivot-first1)], *(last2-1));
292  IT1 store1=first1;
293  for(IT1 it=first1; it!=last1-1; ++it){
294  if(*it<=piv){
295  std::swap(*store1, *it);
296  std::swap(first2[(store1-first1)], first2[(it-first1)]);
297  ++store1;
298  }
299  }
300  std::swap(*(last1-1), *store1);
301  std::swap(*(last2-1), first2[store1-first1]);
302  return store1;
303  }
304 
321  template<class IT1, class IT2, class IT3>
322  IT1 partition3(
323  const IT1& first1,
324  const IT1& last1,
325  const IT2& first2,
326  const IT2& last2,
327  const IT3& first3,
328  const IT3& last3,
329  const IT1& pivot)
330  {
331  typename std::iterator_traits<IT1>::value_type piv(*pivot);
332  std::swap(*pivot, *(last1-1));
333  std::swap(first2[(pivot-first1)], *(last2-1));
334  std::swap(first3[(pivot-first1)], *(last3-1));
335  IT1 store1=first1;
336  for(IT1 it=first1; it!=last1-1; ++it){
337  if(*it<=piv){
338  std::swap(*store1, *it);
339  std::swap(first2[(store1-first1)], first2[(it-first1)]);
340  std::swap(first3[(store1-first1)], first3[(it-first1)]);
341  ++store1;
342  }
343  }
344  std::swap(*(last1-1), *store1);
345  std::swap(*(last2-1), first2[store1-first1]);
346  std::swap(*(last3-1), first3[store1-first1]);
347  return store1;
348  }
349 
360  template<class IT1, class IT2>
361  void quicksort2(
362  const IT1& first1,
363  const IT1& last1,
364  const IT2& first2,
365  const IT2& last2)
366  {
367  typedef typename std::iterator_traits<IT1>::difference_type DT;
368  DT DT1 = Teuchos::OrdinalTraits<DT>::one();
369  if(last1-first1 > DT1){
370  IT1 pivot = getPivot(first1, last1);
371  pivot = partition2(first1, last1, first2, last2, pivot);
372  quicksort2(first1, pivot, first2, first2+(pivot-first1));
373  quicksort2(pivot+1, last1, first2+(pivot-first1)+1, last2);
374  }
375  }
376 
389  template<class IT1, class IT2, class IT3>
390  void quicksort3(
391  const IT1& first1,
392  const IT1& last1,
393  const IT2& first2,
394  const IT2& last2,
395  const IT3& first3,
396  const IT3& last3)
397  {
398  typedef typename std::iterator_traits<IT1>::difference_type DT;
399  DT DT1 = Teuchos::OrdinalTraits<DT>::one();
400  if(last1-first1 > DT1){
401  IT1 pivot = getPivot(first1, last1);
402  pivot = partition3(first1, last1, first2, last2, first3, last3, pivot);
403  quicksort3(first1, pivot, first2, first2+(pivot-first1), first3, first3+(pivot-first1));
404  quicksort3(pivot+1, last1, first2+(pivot-first1)+1, last2, first3+(pivot-first1)+1, last3);
405  }
406  }
407 
414  template<class IT1, class IT2, class IT3>
415  void sh_sort3(
416  const IT1& first1,
417  const IT1& last1,
418  const IT2& first2,
419  const IT2& /* last2 */,
420  const IT3& first3,
421  const IT3& /* last3 */)
422  {
423  typedef typename std::iterator_traits<IT1>::difference_type DT;
424  DT n = last1 - first1;
425  DT m = n / 2;
426  DT z = Teuchos::OrdinalTraits<DT>::zero();
427  while (m > z)
428  {
429  DT max = n - m;
430  for (DT j = 0; j < max; j++)
431  {
432  for (DT k = j; k >= 0; k-=m)
433  {
434  if (first1[k+m] >= first1[k])
435  break;
436  std::swap(first1[k+m], first1[k]);
437  std::swap(first2[k+m], first2[k]);
438  std::swap(first3[k+m], first3[k]);
439  }
440  }
441  m = m/2;
442  }
443  }
444 
451  template<class IT1, class IT2>
452  void sh_sort2(
453  const IT1& first1,
454  const IT1& last1,
455  const IT2& first2,
456  const IT2& /* last2 */)
457  {
458  typedef typename std::iterator_traits<IT1>::difference_type DT;
459  DT n = last1 - first1;
460  DT m = n / 2;
461  DT z = Teuchos::OrdinalTraits<DT>::zero();
462  while (m > z)
463  {
464  DT max = n - m;
465  for (DT j = 0; j < max; j++)
466  {
467  for (DT k = j; k >= 0; k-=m)
468  {
469  if (first1[k+m] >= first1[k])
470  break;
471  std::swap(first1[k+m], first1[k]);
472  std::swap(first2[k+m], first2[k]);
473  }
474  }
475  m = m/2;
476  }
477  }
478 
483  template<typename IT1>
484  std::vector<typename std::iterator_traits<IT1>::difference_type> sort_indexes(IT1 first, IT1 last) {
485 
486  typedef typename std::iterator_traits<IT1>::difference_type DT;
487 
488  DT length = last - first;
489  std::vector<DT> idx(length);
490  std::iota(idx.begin(), idx.end(), 0);
491 
492  std::stable_sort(idx.begin(), idx.end(),
493  [&first](size_t i1, size_t i2) {return first[i1] < first[i2];});
494 
495  return idx;
496  }
497 
501  template<typename IT1, typename IT2>
502  void
503  apply_permutation(
504  IT1 first,
505  IT1 last,
506  IT2 indices)
507  {
508  typedef typename std::iterator_traits<IT1>::difference_type DT;
509  typedef typename std::iterator_traits<IT1>::value_type T;
510 
511  DT n = last - first;
512 
513  Kokkos::View<T*, Kokkos::HostSpace> tmp(Kokkos::view_alloc(Kokkos::WithoutInitializing, "tmp"), n);
514 
515  Kokkos::parallel_for("apply_permutation_1",
516  Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>(0, n),
517  [=](const int j) {
518  tmp(j) = first[indices[j]];
519  });
520  Kokkos::parallel_for("apply_permutation_2",
521  Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>(0, n),
522  [=](const int j) {
523  first[j] = tmp(j);
524  });
525  }
526 
533  template<class IT1, class IT2>
534  void std_sort2(const IT1& first1,
535  const IT1& last1,
536  const IT2& first2,
537  const IT2& last2)
538  {
539  auto indices = sort_indexes(first1, last1);
540  apply_permutation(first1, last1, indices);
541  apply_permutation(first2, last2, indices);
542  }
543 
550  template<class IT1, class IT2, class IT3>
551  void std_sort3(
552  const IT1& first1,
553  const IT1& last1,
554  const IT2& first2,
555  const IT2& last2,
556  const IT3& first3,
557  const IT3& last3)
558  {
559  auto indices = sort_indexes(first1, last1);
560  apply_permutation(first1, last1, indices);
561  apply_permutation(first2, last2, indices);
562  apply_permutation(first3, last3, indices);
563  }
564 
565  } //end namespace SortDetails
566 
567 
586  template<class IT1, class IT2>
587  void sort2(const IT1 &first1, const IT1 &last1, const IT2 &first2, const bool stableSort=false) {
588  // Quicksort uses best-case N log N time whether or not the input
589  // data is sorted. However, the common case in Tpetra is that the
590  // input data are sorted, so we first check whether this is the
591  // case before reverting to Quicksort.
592  if(SortDetails::isAlreadySorted(first1, last1)){
593  return;
594  }
595  if(stableSort)
596  SortDetails::std_sort2(first1, last1, first2, first2+(last1-first1));
597  else
598  SortDetails::sh_sort2(first1, last1, first2, first2+(last1-first1));
599 #ifdef HAVE_TPETRA_DEBUG
600  if(!SortDetails::isAlreadySorted(first1, last1)){
601  std::cout << "Trouble: sort() did not sort !!" << std::endl;
602  return;
603  }
604 #endif
605  }
606 
607 
624  template<class View1, class View2>
625  void sort2(View1 &view1, const size_t &size, View2 &view2) {
626  // NOTE: This assumes the view is host-accessible.
627 
628  // Wrap the views as rcps (this happens to preserve the reference counting, but that doesn't really matter here)
629  Teuchos::ArrayRCP<typename View1::non_const_value_type> view1_rcp = Kokkos::Compat::persistingView(view1, 0, size);
630  Teuchos::ArrayRCP<typename View2::non_const_value_type> view2_rcp = Kokkos::Compat::persistingView(view2, 0, size);
631 
632  sort2(view1_rcp.begin(),view1_rcp.end(),view2_rcp.begin());
633  }
634 
643  template<class View>
644  void sort(View &view, const size_t &size) {
645  // NOTE: This assumes the view is host-accessible.
646 
647  // Wrap the view as rcps (this happens to preserve the reference counting, but that doesn't really matter here)
648  Teuchos::ArrayRCP<typename View::non_const_value_type> view_rcp = Kokkos::Compat::persistingView(view, 0, size);
649 
650  std::sort(view_rcp.begin(),view_rcp.end());
651  }
652 
661  template<class View>
662  void reverse_sort(View &view, const size_t &size) {
663  // NOTE: This assumes the view is host-accessible.
664  // Wrap the view as rcps (this happens to preserve the reference counting, but that doesn't really matter here)
665  Teuchos::ArrayRCP<typename View::non_const_value_type> view_rcp = Kokkos::Compat::persistingView(view, 0, size);
666 
667  std::sort(view_rcp.rbegin(),view_rcp.rend());
668  }
669 
670 
671 
672 
688  template<class IT1, class IT2, class IT3>
689  void sort3(const IT1 &first1, const IT1 &last1, const IT2 &first2,
690  const IT3 &first3, const bool stableSort=false)
691  {
692  // Quicksort uses best-case N log N time whether or not the input
693  // data is sorted. However, the common case in Tpetra is that the
694  // input data are sorted, so we first check whether this is the
695  // case before reverting to Quicksort.
696  if(SortDetails::isAlreadySorted(first1, last1)){
697  return;
698  }
699  if(stableSort)
700  SortDetails::std_sort3(first1, last1, first2, first2+(last1-first1), first3,
701  first3+(last1-first1));
702  else
703  SortDetails::sh_sort3(first1, last1, first2, first2+(last1-first1), first3,
704  first3+(last1-first1));
705 #ifdef HAVE_TPETRA_DEBUG
706  if(!SortDetails::isAlreadySorted(first1, last1)){
707  std::cout << " Trouble sort did not actually sort... !!!!!!" <<
708  std::endl;
709  return;
710  }
711 #endif
712  }
713 
757  template<class IT1, class IT2>
758  void
759  merge2 (IT1& indResultOut, IT2& valResultOut,
760  IT1 indBeg, IT1 indEnd,
761  IT2 valBeg, IT2 /* valEnd */)
762  {
763  if (indBeg == indEnd) {
764  indResultOut = indBeg; // It's allowed for indResultOut to alias indEnd
765  valResultOut = valBeg; // It's allowed for valResultOut to alias valEnd
766  }
767  else {
768  IT1 indResult = indBeg;
769  IT2 valResult = valBeg;
770  if (indBeg != indEnd) {
771  ++indBeg;
772  ++valBeg;
773  while (indBeg != indEnd) {
774  if (*indResult == *indBeg) { // adjacent column indices equal
775  *valResult += *valBeg; // merge entries by adding their values together
776  } else { // adjacent column indices not equal
777  *(++indResult) = *indBeg; // shift over the index
778  *(++valResult) = *valBeg; // shift over the value
779  }
780  ++indBeg;
781  ++valBeg;
782  }
783  ++indResult; // exclusive end of merged result
784  ++valResult; // exclusive end of merged result
785 
786  // mfh 24 Feb 2014: Setting these is technically correct, but
787  // since the resulting variables aren't used after this point,
788  // it may result in a build error.
789  //
790  // indEnd = indResult;
791  // valEnd = valResult;
792  }
793  indResultOut = indResult;
794  valResultOut = valResult;
795  }
796  }
797 
846  template<class IT1, class IT2, class BinaryFunction>
847  void
848  merge2 (IT1& indResultOut, IT2& valResultOut,
849  IT1 indBeg, IT1 indEnd,
850  IT2 valBeg, IT2 valEnd,
851  BinaryFunction f)
852  {
853  if (indBeg == indEnd) {
854  indResultOut = indBeg; // It's allowed for indResultOut to alias indEnd
855  valResultOut = valBeg; // It's allowed for valResultOut to alias valEnd
856  }
857  else {
858  IT1 indResult = indBeg;
859  IT2 valResult = valBeg;
860  if (indBeg != indEnd) {
861  ++indBeg;
862  ++valBeg;
863  while (indBeg != indEnd) {
864  if (*indResult == *indBeg) { // adjacent column indices equal
865  *valResult = f (*valResult, *valBeg); // merge entries via values
866  } else { // adjacent column indices not equal
867  *(++indResult) = *indBeg; // shift over the index
868  *(++valResult) = *valBeg; // shift over the value
869  }
870  ++indBeg;
871  ++valBeg;
872  }
873  ++indResult; // exclusive end of merged result
874  ++valResult; // exclusive end of merged result
875  indEnd = indResult;
876  valEnd = valResult;
877  }
878  indResultOut = indResult;
879  valResultOut = valResult;
880  }
881  }
882 
883 
912  template<class KeyInputIterType, class ValueInputIterType,
913  class KeyOutputIterType, class ValueOutputIterType,
914  class BinaryFunction>
915  void
916  keyValueMerge (KeyInputIterType keyBeg1, KeyInputIterType keyEnd1,
917  ValueInputIterType valBeg1, ValueInputIterType valEnd1,
918  KeyInputIterType keyBeg2, KeyInputIterType keyEnd2,
919  ValueInputIterType valBeg2, ValueInputIterType valEnd2,
920  KeyOutputIterType keyOut, ValueOutputIterType valOut,
921  BinaryFunction f)
922  {
923  while (keyBeg1 != keyEnd1 && keyBeg2 != keyEnd2) {
924  if (*keyBeg1 < *keyBeg2) {
925  *keyOut++ = *keyBeg1++;
926  *valOut++ = *valBeg1++;
927  } else if (*keyBeg1 > *keyBeg2) {
928  *keyOut++ = *keyBeg2++;
929  *valOut++ = *valBeg2++;
930  } else { // *keyBeg1 == *keyBeg2
931  *keyOut++ = *keyBeg1;
932  *valOut++ = f (*valBeg1, *valBeg2);
933  ++keyBeg1;
934  ++valBeg1;
935  ++keyBeg2;
936  ++valBeg2;
937  }
938  }
939  // At most one of the two sequences will be nonempty.
940  std::copy (keyBeg1, keyEnd1, keyOut);
941  std::copy (valBeg1, valEnd1, valOut);
942  std::copy (keyBeg2, keyEnd2, keyOut);
943  std::copy (valBeg2, valEnd2, valOut);
944  }
945 
946  template<class KeyInputIterType>
947  size_t
948  keyMergeCount (KeyInputIterType keyBeg1, KeyInputIterType keyEnd1,
949  KeyInputIterType keyBeg2, KeyInputIterType keyEnd2)
950  {
951  size_t count = 0;
952  while (keyBeg1 != keyEnd1 && keyBeg2 != keyEnd2) {
953  if (*keyBeg1 < *keyBeg2) {
954  ++keyBeg1;
955  } else if (*keyBeg1 > *keyBeg2) {
956  ++keyBeg2;
957  } else { // *keyBeg1 == *keyBeg2
958  ++keyBeg1;
959  ++keyBeg2;
960  }
961  ++count;
962  }
963  // At most one of the two sequences will be nonempty.
964  count += static_cast<size_t> (keyEnd1 - keyBeg1) +
965  static_cast<size_t> (keyEnd2 - keyBeg2);
966  return count;
967  }
968 
969  namespace Details {
989  bool
990  congruent (const Teuchos::Comm<int>& comm1,
991  const Teuchos::Comm<int>& comm2);
992 
1002  template<class DualViewType>
1003  Teuchos::ArrayView<typename DualViewType::t_dev::value_type>
1004  getArrayViewFromDualView (const DualViewType& x)
1005  {
1006  static_assert (static_cast<int> (DualViewType::t_dev::rank) == 1,
1007  "The input DualView must have rank 1.");
1008  TEUCHOS_TEST_FOR_EXCEPTION
1009  (x.need_sync_host (), std::logic_error, "The "
1010  "input Kokkos::DualView was most recently modified on device, but this "
1011  "function needs the host view of the data to be the most recently "
1012  "modified.");
1013 
1014  auto x_host = x.view_host ();
1015  typedef typename DualViewType::t_dev::value_type value_type;
1016  // mfh 15 Jan 2019: In debug mode, Teuchos::ArrayView's
1017  // constructor throws if the pointer is nonnull but the length
1018  // is nonpositive.
1019  const auto len = x_host.extent (0);
1020  return Teuchos::ArrayView<value_type> (len != 0 ? x_host.data () : nullptr,
1021  len);
1022  }
1023 
1040  template<class T, class DT>
1041  Kokkos::DualView<T*, DT>
1042  getDualViewCopyFromArrayView (const Teuchos::ArrayView<const T>& x_av,
1043  const char label[],
1044  const bool leaveOnHost)
1045  {
1046  using Kokkos::MemoryUnmanaged;
1047  typedef typename DT::memory_space DMS;
1048  typedef typename DT::execution_space execution_space;
1049  typedef Kokkos::HostSpace HMS;
1050 
1051  const size_t len = static_cast<size_t> (x_av.size ());
1052  Kokkos::View<const T*, HMS, MemoryUnmanaged> x_in (x_av.getRawPtr (), len);
1053  Kokkos::DualView<T*, DT> x_out (label, len);
1054  if (leaveOnHost) {
1055  x_out.modify_host ();
1056  // DEEP_COPY REVIEW - NOT TESTED FOR CUDA BUILD
1057  Kokkos::deep_copy (x_out.view_host (), x_in);
1058  }
1059  else {
1060  x_out.template modify<DMS> ();
1061  // DEEP_COPY REVIEW - HOST-TO-DEVICE
1062  Kokkos::deep_copy (execution_space(), x_out.template view<DMS> (), x_in);
1063  }
1064  return x_out;
1065  }
1066 
1074  template<class DualViewType>
1075  std::string dualViewStatusToString (const DualViewType& dv, const char name[])
1076  {
1077  const auto host = dv.need_sync_device();
1078  const auto dev = dv.need_sync_host();
1079 
1080  std::ostringstream os;
1081  os << name << ": {size: " << dv.extent (0)
1082  << ", sync: {host: " << host << ", dev: " << dev << "}";
1083  return os.str ();
1084  }
1085 
1090  template<class ArrayType>
1091  void
1092  verbosePrintArray(std::ostream& out,
1093  const ArrayType& x,
1094  const char name[],
1095  const size_t maxNumToPrint)
1096  {
1097  out << name << ": [";
1098 
1099  const size_t numEnt(x.size());
1100  if (maxNumToPrint == 0) {
1101  if (numEnt != 0) {
1102  out << "...";
1103  }
1104  }
1105  else {
1106  const size_t numToPrint = numEnt > maxNumToPrint ?
1107  maxNumToPrint : numEnt;
1108  size_t k = 0;
1109  for ( ; k < numToPrint; ++k) {
1110  out << x[k];
1111  if (k + size_t(1) < numToPrint) {
1112  out << ", ";
1113  }
1114  }
1115  if (k < numEnt) {
1116  out << ", ...";
1117  }
1118  }
1119  out << "]";
1120  }
1121 
1125  std::unique_ptr<std::string>
1126  createPrefix(const int myRank,
1127  const char prefix[]);
1128 
1136  std::unique_ptr<std::string>
1137  createPrefix(const Teuchos::Comm<int>* comm,
1138  const char functionName[]);
1139 
1146  std::unique_ptr<std::string>
1147  createPrefix(const Teuchos::Comm<int>*,
1148  const char className[],
1149  const char methodName[]);
1150 
1151  } // namespace Details
1152 } // namespace Tpetra
1153 
1154 #endif // TPETRA_UTIL_HPP
bool congruent(const Teuchos::Comm< int > &comm1, const Teuchos::Comm< int > &comm2)
Whether the two communicators are congruent.
Definition: Tpetra_Util.cpp:64
void merge2(IT1 &indResultOut, IT2 &valResultOut, IT1 indBeg, IT1 indEnd, IT2 valBeg, IT2)
Merge values in place, additively, with the same index.
void sort3(const IT1 &first1, const IT1 &last1, const IT2 &first2, const IT3 &first3, const bool stableSort=false)
Sort the first array, and apply the same permutation to the second and third arrays.
void verbosePrintArray(std::ostream &out, const ArrayType &x, const char name[], const size_t maxNumToPrint)
Print min(x.size(), maxNumToPrint) entries of x.
void sort(View &view, const size_t &size)
Convenience wrapper for std::sort for host-accessible views.
void deep_copy(MultiVector< DS, DL, DG, DN > &dst, const MultiVector< SS, SL, SG, SN > &src)
Copy the contents of the MultiVector src into dst.
Kokkos::DualView< T *, DT > getDualViewCopyFromArrayView(const Teuchos::ArrayView< const T > &x_av, const char label[], const bool leaveOnHost)
Get a 1-D Kokkos::DualView which is a deep copy of the input Teuchos::ArrayView (which views host mem...
void keyValueMerge(KeyInputIterType keyBeg1, KeyInputIterType keyEnd1, ValueInputIterType valBeg1, ValueInputIterType valEnd1, KeyInputIterType keyBeg2, KeyInputIterType keyEnd2, ValueInputIterType valBeg2, ValueInputIterType valEnd2, KeyOutputIterType keyOut, ValueOutputIterType valOut, BinaryFunction f)
Merge two sorted (by keys) sequences of unique (key,value) pairs by combining pairs with equal keys...
void sort2(const IT1 &first1, const IT1 &last1, const IT2 &first2, const bool stableSort=false)
Sort the first array, and apply the resulting permutation to the second array.
MapType::iterator efficientAddOrUpdate(MapType &m, const KeyArgType &k, const ValueArgType &v)
Efficiently insert or replace an entry in an std::map.
std::string dualViewStatusToString(const DualViewType &dv, const char name[])
Return the status of the given Kokkos::DualView, as a human-readable string.
Teuchos::ArrayView< typename DualViewType::t_dev::value_type > getArrayViewFromDualView(const DualViewType &x)
Get a Teuchos::ArrayView which views the host Kokkos::View of the input 1-D Kokkos::DualView.
void reverse_sort(View &view, const size_t &size)
Convenience wrapper for a reversed std::sort for host-accessible views.
std::unique_ptr< std::string > createPrefix(const int myRank, const char prefix[])
Create string prefix for each line of verbose output.