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"
66 #if defined(HAVE_TPETRA_PRINT_EFFICIENCY_WARNINGS)
67 #define TPETRA_EFFICIENCY_WARNING(throw_exception_test,msg) \
93 const bool tpetraEfficiencyWarningTest = (throw_exception_test); \
94 if (tpetraEfficiencyWarningTest) { \
95 std::ostringstream errStream; \
96 errStream << Teuchos::typeName(*this) << ":" << std::endl; \
97 errStream << "Efficiency warning: " << #throw_exception_test << std::endl; \
99 std::string err = errStream.str(); \
100 if (TPETRA_PRINTS_EFFICIENCY_WARNINGS && tpetraEfficiencyWarningTest) { \
101 std::cerr << err << std::endl; \
106 #define TPETRA_EFFICIENCY_WARNING(throw_exception_test,msg)
134 #if defined(HAVE_TPETRA_THROW_ABUSE_WARNINGS) || defined(HAVE_TPETRA_PRINT_ABUSE_WARNINGS)
135 #define TPETRA_ABUSE_WARNING(throw_exception_test,Exception,msg) \
138 std::ostringstream errStream; \
139 errStream << Teuchos::typeName(*this) << msg; \
140 std::string err = errStream.str(); \
141 if (TPETRA_PRINTS_ABUSE_WARNINGS && (throw_exception_test)) { \
142 std::cerr << err << std::endl; \
144 TEUCHOS_TEST_FOR_EXCEPTION(TPETRA_THROWS_ABUSE_WARNINGS && (throw_exception_test), Exception, err); \
147 #define TPETRA_ABUSE_WARNING(throw_exception_test,Exception,msg)
180 #define SHARED_TEST_FOR_EXCEPTION(throw_exception_test,Exception,msg,comm) \
182 using Teuchos::outArg; \
183 const int lcl_throw_exception = (throw_exception_test) ? Teuchos::rank(comm)+1 : 0; \
185 Teuchos::reduceAll(comm,Teuchos::REDUCE_MAX,lcl_throw_exception,outArg(gbl_throw)); \
186 TEUCHOS_TEST_FOR_EXCEPTION(gbl_throw,Exception, \
187 msg << " Failure on at least one process, including process " << gbl_throw-1 << "." << std::endl); \
203 template<
typename MapType,
typename KeyArgType,
typename ValueArgType>
205 const KeyArgType & k,
206 const ValueArgType & v)
208 typename MapType::iterator lb = m.lower_bound(k);
209 if(lb != m.end() && !(m.key_comp()(k, lb->first))) {
214 typedef typename MapType::value_type MVT;
215 return(m.insert(lb, MVT(k, v)));
226 namespace SortDetails {
237 bool isAlreadySorted(
const IT1& first,
const IT1& last){
238 typedef typename std::iterator_traits<IT1>::difference_type DT;
239 DT myit = Teuchos::OrdinalTraits<DT>::one();
240 const DT sz = last - first;
241 for(;myit < sz; ++myit){
242 if(first[myit] < first[myit-1]){
259 IT getPivot(
const IT& first,
const IT& last){
260 IT pivot(first+(last-first)/2);
261 if(*first<=*pivot && *(last-1)<=*first) pivot=first;
262 else if(*(last-1)<=*pivot && *first<= *(last-1)) pivot = last-1;
280 template<
class IT1,
class IT2>
288 typename std::iterator_traits<IT1>::value_type piv(*pivot);
289 std::swap(*pivot, *(last1-1));
290 std::swap(first2[(pivot-first1)], *(last2-1));
292 for(IT1 it=first1; it!=last1-1; ++it){
294 std::swap(*store1, *it);
295 std::swap(first2[(store1-first1)], first2[(it-first1)]);
299 std::swap(*(last1-1), *store1);
300 std::swap(*(last2-1), first2[store1-first1]);
320 template<
class IT1,
class IT2,
class IT3>
330 typename std::iterator_traits<IT1>::value_type piv(*pivot);
331 std::swap(*pivot, *(last1-1));
332 std::swap(first2[(pivot-first1)], *(last2-1));
333 std::swap(first3[(pivot-first1)], *(last3-1));
335 for(IT1 it=first1; it!=last1-1; ++it){
337 std::swap(*store1, *it);
338 std::swap(first2[(store1-first1)], first2[(it-first1)]);
339 std::swap(first3[(store1-first1)], first3[(it-first1)]);
343 std::swap(*(last1-1), *store1);
344 std::swap(*(last2-1), first2[store1-first1]);
345 std::swap(*(last3-1), first3[store1-first1]);
359 template<
class IT1,
class IT2>
366 typedef typename std::iterator_traits<IT1>::difference_type DT;
367 DT DT1 = Teuchos::OrdinalTraits<DT>::one();
368 if(last1-first1 > DT1){
369 IT1 pivot = getPivot(first1, last1);
370 pivot = partition2(first1, last1, first2, last2, pivot);
371 quicksort2(first1, pivot, first2, first2+(pivot-first1));
372 quicksort2(pivot+1, last1, first2+(pivot-first1)+1, last2);
388 template<
class IT1,
class IT2,
class IT3>
397 typedef typename std::iterator_traits<IT1>::difference_type DT;
398 DT DT1 = Teuchos::OrdinalTraits<DT>::one();
399 if(last1-first1 > DT1){
400 IT1 pivot = getPivot(first1, last1);
401 pivot = partition3(first1, last1, first2, last2, first3, last3, pivot);
402 quicksort3(first1, pivot, first2, first2+(pivot-first1), first3, first3+(pivot-first1));
403 quicksort3(pivot+1, last1, first2+(pivot-first1)+1, last2, first3+(pivot-first1)+1, last3);
413 template<
class IT1,
class IT2,
class IT3>
422 typedef typename std::iterator_traits<IT1>::difference_type DT;
423 DT n = last1 - first1;
425 DT z = Teuchos::OrdinalTraits<DT>::zero();
429 for (DT j = 0; j < max; j++)
431 for (DT k = j; k >= 0; k-=m)
433 if (first1[k+m] >= first1[k])
435 std::swap(first1[k+m], first1[k]);
436 std::swap(first2[k+m], first2[k]);
437 std::swap(first3[k+m], first3[k]);
450 template<
class IT1,
class IT2>
457 typedef typename std::iterator_traits<IT1>::difference_type DT;
458 DT n = last1 - first1;
460 DT z = Teuchos::OrdinalTraits<DT>::zero();
464 for (DT j = 0; j < max; j++)
466 for (DT k = j; k >= 0; k-=m)
468 if (first1[k+m] >= first1[k])
470 std::swap(first1[k+m], first1[k]);
471 std::swap(first2[k+m], first2[k]);
499 template<
class IT1,
class IT2>
500 void sort2(
const IT1 &first1,
const IT1 &last1,
const IT2 &first2) {
505 if(SortDetails::isAlreadySorted(first1, last1)){
508 SortDetails::sh_sort2(first1, last1, first2, first2+(last1-first1));
509 #ifdef HAVE_TPETRA_DEBUG
510 if(!SortDetails::isAlreadySorted(first1, last1)){
511 std::cout <<
"Trouble: sort() did not sort !!" << std::endl;
534 template<
class View1,
class View2>
535 void sort2(View1 &view1,
const size_t &size, View2 &view2) {
539 Teuchos::ArrayRCP<typename View1::non_const_value_type> view1_rcp = Kokkos::Compat::persistingView(view1, 0, size);
540 Teuchos::ArrayRCP<typename View2::non_const_value_type> view2_rcp = Kokkos::Compat::persistingView(view2, 0, size);
542 sort2(view1_rcp.begin(),view1_rcp.end(),view2_rcp.begin());
554 void sort(View &view,
const size_t &size) {
558 Teuchos::ArrayRCP<typename View::non_const_value_type> view_rcp = Kokkos::Compat::persistingView(view, 0, size);
560 std::sort(view_rcp.begin(),view_rcp.end());
575 Teuchos::ArrayRCP<typename View::non_const_value_type> view_rcp = Kokkos::Compat::persistingView(view, 0, size);
577 std::sort(view_rcp.rbegin(),view_rcp.rend());
598 template<
class IT1,
class IT2,
class IT3>
599 void sort3(
const IT1 &first1,
const IT1 &last1,
const IT2 &first2,
606 if(SortDetails::isAlreadySorted(first1, last1)){
609 SortDetails::sh_sort3(first1, last1, first2, first2+(last1-first1), first3,
610 first3+(last1-first1));
612 #ifdef HAVE_TPETRA_DEBUG
613 if(!SortDetails::isAlreadySorted(first1, last1)){
614 std::cout <<
" Trouble sort did not actually sort... !!!!!!" <<
664 template<
class IT1,
class IT2>
666 merge2 (IT1& indResultOut, IT2& valResultOut,
667 IT1 indBeg, IT1 indEnd,
670 if (indBeg == indEnd) {
671 indResultOut = indBeg;
672 valResultOut = valBeg;
675 IT1 indResult = indBeg;
676 IT2 valResult = valBeg;
677 if (indBeg != indEnd) {
680 while (indBeg != indEnd) {
681 if (*indResult == *indBeg) {
682 *valResult += *valBeg;
684 *(++indResult) = *indBeg;
685 *(++valResult) = *valBeg;
700 indResultOut = indResult;
701 valResultOut = valResult;
753 template<
class IT1,
class IT2,
class BinaryFunction>
755 merge2 (IT1& indResultOut, IT2& valResultOut,
756 IT1 indBeg, IT1 indEnd,
757 IT2 valBeg, IT2 valEnd,
760 if (indBeg == indEnd) {
761 indResultOut = indBeg;
762 valResultOut = valBeg;
765 IT1 indResult = indBeg;
766 IT2 valResult = valBeg;
767 if (indBeg != indEnd) {
770 while (indBeg != indEnd) {
771 if (*indResult == *indBeg) {
772 *valResult = f (*valResult, *valBeg);
774 *(++indResult) = *indBeg;
775 *(++valResult) = *valBeg;
785 indResultOut = indResult;
786 valResultOut = valResult;
819 template<
class KeyInputIterType,
class ValueInputIterType,
820 class KeyOutputIterType,
class ValueOutputIterType,
821 class BinaryFunction>
824 ValueInputIterType valBeg1, ValueInputIterType valEnd1,
825 KeyInputIterType keyBeg2, KeyInputIterType keyEnd2,
826 ValueInputIterType valBeg2, ValueInputIterType valEnd2,
827 KeyOutputIterType keyOut, ValueOutputIterType valOut,
830 while (keyBeg1 != keyEnd1 && keyBeg2 != keyEnd2) {
831 if (*keyBeg1 < *keyBeg2) {
832 *keyOut++ = *keyBeg1++;
833 *valOut++ = *valBeg1++;
834 }
else if (*keyBeg1 > *keyBeg2) {
835 *keyOut++ = *keyBeg2++;
836 *valOut++ = *valBeg2++;
838 *keyOut++ = *keyBeg1;
839 *valOut++ = f (*valBeg1, *valBeg2);
847 std::copy (keyBeg1, keyEnd1, keyOut);
848 std::copy (valBeg1, valEnd1, valOut);
849 std::copy (keyBeg2, keyEnd2, keyOut);
850 std::copy (valBeg2, valEnd2, valOut);
853 template<
class KeyInputIterType>
855 keyMergeCount (KeyInputIterType keyBeg1, KeyInputIterType keyEnd1,
856 KeyInputIterType keyBeg2, KeyInputIterType keyEnd2)
859 while (keyBeg1 != keyEnd1 && keyBeg2 != keyEnd2) {
860 if (*keyBeg1 < *keyBeg2) {
862 }
else if (*keyBeg1 > *keyBeg2) {
871 count +=
static_cast<size_t> (keyEnd1 - keyBeg1) +
872 static_cast<size_t> (keyEnd2 - keyBeg2);
897 congruent (
const Teuchos::Comm<int>& comm1,
898 const Teuchos::Comm<int>& comm2);
909 template<
class DualViewType>
910 Teuchos::ArrayView<typename DualViewType::t_dev::value_type>
913 static_assert (static_cast<int> (DualViewType::t_dev::rank) == 1,
914 "The input DualView must have rank 1.");
915 TEUCHOS_TEST_FOR_EXCEPTION
916 (x.need_sync_host (), std::logic_error,
"The "
917 "input Kokkos::DualView was most recently modified on device, but this "
918 "function needs the host view of the data to be the most recently "
921 auto x_host = x.view_host ();
922 typedef typename DualViewType::t_dev::value_type value_type;
926 const auto len = x_host.extent (0);
927 return Teuchos::ArrayView<value_type> (len != 0 ? x_host.data () :
nullptr,
947 template<
class T,
class DT>
948 Kokkos::DualView<T*, DT>
951 const bool leaveOnHost)
953 using Kokkos::MemoryUnmanaged;
954 typedef typename DT::memory_space DMS;
955 typedef typename DT::execution_space execution_space;
956 typedef Kokkos::HostSpace HMS;
958 const size_t len =
static_cast<size_t> (x_av.size ());
959 Kokkos::View<const T*, HMS, MemoryUnmanaged> x_in (x_av.getRawPtr (), len);
960 Kokkos::DualView<T*, DT> x_out (label, len);
962 x_out.modify_host ();
967 x_out.template modify<DMS> ();
981 template<
class DualViewType>
984 const auto host = dv.need_sync_device();
985 const auto dev = dv.need_sync_host();
987 std::ostringstream os;
988 os << name <<
": {size: " << dv.extent (0)
989 <<
", sync: {host: " << host <<
", dev: " << dev <<
"}";
997 template<
class ArrayType>
1002 const size_t maxNumToPrint)
1004 out << name <<
": [";
1006 const size_t numEnt(x.size());
1007 if (maxNumToPrint == 0) {
1013 const size_t numToPrint = numEnt > maxNumToPrint ?
1014 maxNumToPrint : numEnt;
1016 for ( ; k < numToPrint; ++k) {
1018 if (k +
size_t(1) < numToPrint) {
1032 std::unique_ptr<std::string>
1034 const char prefix[]);
1043 std::unique_ptr<std::string>
1045 const char functionName[]);
1053 std::unique_ptr<std::string>
1055 const char className[],
1056 const char methodName[]);
1061 #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)
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...
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.
void sort2(const IT1 &first1, const IT1 &last1, const IT2 &first2)
Sort the first array, and apply the resulting permutation to the second array.
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.