40 #ifndef TPETRA_UTIL_HPP
41 #define TPETRA_UTIL_HPP
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"
67 #if defined(HAVE_TPETRA_PRINT_EFFICIENCY_WARNINGS)
68 #define TPETRA_EFFICIENCY_WARNING(throw_exception_test,msg) \
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; \
100 std::string err = errStream.str(); \
101 if (TPETRA_PRINTS_EFFICIENCY_WARNINGS && tpetraEfficiencyWarningTest) { \
102 std::cerr << err << std::endl; \
107 #define TPETRA_EFFICIENCY_WARNING(throw_exception_test,msg)
135 #if defined(HAVE_TPETRA_THROW_ABUSE_WARNINGS) || defined(HAVE_TPETRA_PRINT_ABUSE_WARNINGS)
136 #define TPETRA_ABUSE_WARNING(throw_exception_test,Exception,msg) \
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; \
145 TEUCHOS_TEST_FOR_EXCEPTION(TPETRA_THROWS_ABUSE_WARNINGS && (throw_exception_test), Exception, err); \
148 #define TPETRA_ABUSE_WARNING(throw_exception_test,Exception,msg)
181 #define SHARED_TEST_FOR_EXCEPTION(throw_exception_test,Exception,msg,comm) \
183 using Teuchos::outArg; \
184 const int lcl_throw_exception = (throw_exception_test) ? Teuchos::rank(comm)+1 : 0; \
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); \
204 template<
typename MapType,
typename KeyArgType,
typename ValueArgType>
206 const KeyArgType & k,
207 const ValueArgType & v)
209 typename MapType::iterator lb = m.lower_bound(k);
210 if(lb != m.end() && !(m.key_comp()(k, lb->first))) {
215 typedef typename MapType::value_type MVT;
216 return(m.insert(lb, MVT(k, v)));
227 namespace SortDetails {
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]){
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;
281 template<
class IT1,
class IT2>
289 typename std::iterator_traits<IT1>::value_type piv(*pivot);
290 std::swap(*pivot, *(last1-1));
291 std::swap(first2[(pivot-first1)], *(last2-1));
293 for(IT1 it=first1; it!=last1-1; ++it){
295 std::swap(*store1, *it);
296 std::swap(first2[(store1-first1)], first2[(it-first1)]);
300 std::swap(*(last1-1), *store1);
301 std::swap(*(last2-1), first2[store1-first1]);
321 template<
class IT1,
class IT2,
class IT3>
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));
336 for(IT1 it=first1; it!=last1-1; ++it){
338 std::swap(*store1, *it);
339 std::swap(first2[(store1-first1)], first2[(it-first1)]);
340 std::swap(first3[(store1-first1)], first3[(it-first1)]);
344 std::swap(*(last1-1), *store1);
345 std::swap(*(last2-1), first2[store1-first1]);
346 std::swap(*(last3-1), first3[store1-first1]);
360 template<
class IT1,
class IT2>
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);
389 template<
class IT1,
class IT2,
class IT3>
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);
414 template<
class IT1,
class IT2,
class IT3>
423 typedef typename std::iterator_traits<IT1>::difference_type DT;
424 DT n = last1 - first1;
426 DT z = Teuchos::OrdinalTraits<DT>::zero();
430 for (DT j = 0; j < max; j++)
432 for (DT k = j; k >= 0; k-=m)
434 if (first1[k+m] >= first1[k])
436 std::swap(first1[k+m], first1[k]);
437 std::swap(first2[k+m], first2[k]);
438 std::swap(first3[k+m], first3[k]);
451 template<
class IT1,
class IT2>
458 typedef typename std::iterator_traits<IT1>::difference_type DT;
459 DT n = last1 - first1;
461 DT z = Teuchos::OrdinalTraits<DT>::zero();
465 for (DT j = 0; j < max; j++)
467 for (DT k = j; k >= 0; k-=m)
469 if (first1[k+m] >= first1[k])
471 std::swap(first1[k+m], first1[k]);
472 std::swap(first2[k+m], first2[k]);
483 template<
typename IT1>
484 std::vector<typename std::iterator_traits<IT1>::difference_type> sort_indexes(IT1 first, IT1 last) {
486 typedef typename std::iterator_traits<IT1>::difference_type DT;
488 DT length = last - first;
489 std::vector<DT> idx(length);
490 std::iota(idx.begin(), idx.end(), 0);
492 std::stable_sort(idx.begin(), idx.end(),
493 [&first](
size_t i1,
size_t i2) {
return first[i1] < first[i2];});
501 template<
typename IT1,
typename IT2>
508 typedef typename std::iterator_traits<IT1>::difference_type DT;
509 typedef typename std::iterator_traits<IT1>::value_type T;
513 Kokkos::View<T*, Kokkos::HostSpace> tmp(Kokkos::view_alloc(Kokkos::WithoutInitializing,
"tmp"), n);
515 Kokkos::parallel_for(
"apply_permutation_1",
516 Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>(0, n),
518 tmp(j) = first[indices[j]];
520 Kokkos::parallel_for(
"apply_permutation_2",
521 Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>(0, n),
533 template<
class IT1,
class IT2>
534 void std_sort2(
const IT1& first1,
539 auto indices = sort_indexes(first1, last1);
540 apply_permutation(first1, last1, indices);
541 apply_permutation(first2, last2, indices);
550 template<
class IT1,
class IT2,
class IT3>
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);
586 template<
class IT1,
class IT2>
587 void sort2(
const IT1 &first1,
const IT1 &last1,
const IT2 &first2,
const bool stableSort=
false) {
592 if(SortDetails::isAlreadySorted(first1, last1)){
596 SortDetails::std_sort2(first1, last1, first2, first2+(last1-first1));
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;
624 template<
class View1,
class View2>
625 void sort2(View1 &view1,
const size_t &size, View2 &view2) {
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);
632 sort2(view1_rcp.begin(),view1_rcp.end(),view2_rcp.begin());
644 void sort(View &view,
const size_t &size) {
648 Teuchos::ArrayRCP<typename View::non_const_value_type> view_rcp = Kokkos::Compat::persistingView(view, 0, size);
650 std::sort(view_rcp.begin(),view_rcp.end());
665 Teuchos::ArrayRCP<typename View::non_const_value_type> view_rcp = Kokkos::Compat::persistingView(view, 0, size);
667 std::sort(view_rcp.rbegin(),view_rcp.rend());
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)
696 if(SortDetails::isAlreadySorted(first1, last1)){
700 SortDetails::std_sort3(first1, last1, first2, first2+(last1-first1), first3,
701 first3+(last1-first1));
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... !!!!!!" <<
757 template<
class IT1,
class IT2>
759 merge2 (IT1& indResultOut, IT2& valResultOut,
760 IT1 indBeg, IT1 indEnd,
763 if (indBeg == indEnd) {
764 indResultOut = indBeg;
765 valResultOut = valBeg;
768 IT1 indResult = indBeg;
769 IT2 valResult = valBeg;
770 if (indBeg != indEnd) {
773 while (indBeg != indEnd) {
774 if (*indResult == *indBeg) {
775 *valResult += *valBeg;
777 *(++indResult) = *indBeg;
778 *(++valResult) = *valBeg;
793 indResultOut = indResult;
794 valResultOut = valResult;
846 template<
class IT1,
class IT2,
class BinaryFunction>
848 merge2 (IT1& indResultOut, IT2& valResultOut,
849 IT1 indBeg, IT1 indEnd,
850 IT2 valBeg, IT2 valEnd,
853 if (indBeg == indEnd) {
854 indResultOut = indBeg;
855 valResultOut = valBeg;
858 IT1 indResult = indBeg;
859 IT2 valResult = valBeg;
860 if (indBeg != indEnd) {
863 while (indBeg != indEnd) {
864 if (*indResult == *indBeg) {
865 *valResult = f (*valResult, *valBeg);
867 *(++indResult) = *indBeg;
868 *(++valResult) = *valBeg;
878 indResultOut = indResult;
879 valResultOut = valResult;
912 template<
class KeyInputIterType,
class ValueInputIterType,
913 class KeyOutputIterType,
class ValueOutputIterType,
914 class BinaryFunction>
917 ValueInputIterType valBeg1, ValueInputIterType valEnd1,
918 KeyInputIterType keyBeg2, KeyInputIterType keyEnd2,
919 ValueInputIterType valBeg2, ValueInputIterType valEnd2,
920 KeyOutputIterType keyOut, ValueOutputIterType valOut,
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++;
931 *keyOut++ = *keyBeg1;
932 *valOut++ = f (*valBeg1, *valBeg2);
940 std::copy (keyBeg1, keyEnd1, keyOut);
941 std::copy (valBeg1, valEnd1, valOut);
942 std::copy (keyBeg2, keyEnd2, keyOut);
943 std::copy (valBeg2, valEnd2, valOut);
946 template<
class KeyInputIterType>
948 keyMergeCount (KeyInputIterType keyBeg1, KeyInputIterType keyEnd1,
949 KeyInputIterType keyBeg2, KeyInputIterType keyEnd2)
952 while (keyBeg1 != keyEnd1 && keyBeg2 != keyEnd2) {
953 if (*keyBeg1 < *keyBeg2) {
955 }
else if (*keyBeg1 > *keyBeg2) {
964 count +=
static_cast<size_t> (keyEnd1 - keyBeg1) +
965 static_cast<size_t> (keyEnd2 - keyBeg2);
990 congruent (
const Teuchos::Comm<int>& comm1,
991 const Teuchos::Comm<int>& comm2);
1002 template<
class DualViewType>
1003 Teuchos::ArrayView<typename DualViewType::t_dev::value_type>
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 "
1014 auto x_host = x.view_host ();
1015 typedef typename DualViewType::t_dev::value_type value_type;
1019 const auto len = x_host.extent (0);
1020 return Teuchos::ArrayView<value_type> (len != 0 ? x_host.data () :
nullptr,
1040 template<
class T,
class DT>
1041 Kokkos::DualView<T*, DT>
1044 const bool leaveOnHost)
1046 using Kokkos::MemoryUnmanaged;
1047 typedef typename DT::memory_space DMS;
1048 typedef typename DT::execution_space execution_space;
1049 typedef Kokkos::HostSpace HMS;
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);
1055 x_out.modify_host ();
1060 x_out.template modify<DMS> ();
1074 template<
class DualViewType>
1077 const auto host = dv.need_sync_device();
1078 const auto dev = dv.need_sync_host();
1080 std::ostringstream os;
1081 os << name <<
": {size: " << dv.extent (0)
1082 <<
", sync: {host: " << host <<
", dev: " << dev <<
"}";
1090 template<
class ArrayType>
1095 const size_t maxNumToPrint)
1097 out << name <<
": [";
1099 const size_t numEnt(x.size());
1100 if (maxNumToPrint == 0) {
1106 const size_t numToPrint = numEnt > maxNumToPrint ?
1107 maxNumToPrint : numEnt;
1109 for ( ; k < numToPrint; ++k) {
1111 if (k +
size_t(1) < numToPrint) {
1125 std::unique_ptr<std::string>
1127 const char prefix[]);
1136 std::unique_ptr<std::string>
1138 const char functionName[]);
1146 std::unique_ptr<std::string>
1148 const char className[],
1149 const char methodName[]);
1154 #endif // TPETRA_UTIL_HPP
bool congruent(const Teuchos::Comm< int > &comm1, const Teuchos::Comm< int > &comm2)
Whether the two communicators are congruent.
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.