Tpetra_Util.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//          Tpetra: Templated Linear Algebra Services Package
//
// Copyright 2008 NTESS and the Tpetra contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
#ifndef TPETRA_UTIL_HPP
#define TPETRA_UTIL_HPP

 
#include "Tpetra_ConfigDefs.hpp"
#include "Kokkos_DualView.hpp"
#include "KokkosCompat_View.hpp"
#include "Teuchos_Assert.hpp"
#include "Teuchos_CommHelpers.hpp"
#include "Teuchos_OrdinalTraits.hpp"
#include "Teuchos_TypeNameTraits.hpp"
#include "Teuchos_Utils.hpp"
#include <numeric>
#include <algorithm>
#include <iterator>
#include <memory>
#include <ostream>
#include <sstream>
 
#if defined(HAVE_TPETRA_PRINT_EFFICIENCY_WARNINGS)
#define TPETRA_EFFICIENCY_WARNING(throw_exception_test, msg)                     \
  {                                                                              \
    const bool tpetraEfficiencyWarningTest = (throw_exception_test);             \
    if (tpetraEfficiencyWarningTest) {                                           \
      std::ostringstream errStream;                                              \
      errStream << Teuchos::typeName(*this) << ":" << std::endl;                 \
      errStream << "Efficiency warning: " << #throw_exception_test << std::endl; \
      errStream << msg;                                                          \
      std::string err = errStream.str();                                         \
      if (TPETRA_PRINTS_EFFICIENCY_WARNINGS && tpetraEfficiencyWarningTest) {    \
        std::cerr << err << std::endl;                                           \
      }                                                                          \
    }                                                                            \
  }
#else
#define TPETRA_EFFICIENCY_WARNING(throw_exception_test, msg)
#endif
 
// handle an abuse warning, according to HAVE_TPETRA_THROW_ABUSE_WARNINGS and HAVE_TPETRA_PRINT_ABUSE_WARNINGS
#if defined(HAVE_TPETRA_THROW_ABUSE_WARNINGS) || defined(HAVE_TPETRA_PRINT_ABUSE_WARNINGS)
#define TPETRA_ABUSE_WARNING(throw_exception_test, Exception, msg)                                      \
  {                                                                                                     \
    std::ostringstream errStream;                                                                       \
    errStream << Teuchos::typeName(*this) << msg;                                                       \
    std::string err = errStream.str();                                                                  \
    if (TPETRA_PRINTS_ABUSE_WARNINGS && (throw_exception_test)) {                                       \
      std::cerr << err << std::endl;                                                                    \
    }                                                                                                   \
    TEUCHOS_TEST_FOR_EXCEPTION(TPETRA_THROWS_ABUSE_WARNINGS && (throw_exception_test), Exception, err); \
  }
#else
#define TPETRA_ABUSE_WARNING(throw_exception_test, Exception, msg)
#endif

#define SHARED_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg, comm)                                                       \
  {                                                                                                                                 \
    using Teuchos::outArg;                                                                                                          \
    const int lcl_throw_exception = (throw_exception_test) ? Teuchos::rank(comm) + 1 : 0;                                           \
    int gbl_throw;                                                                                                                  \
    Teuchos::reduceAll(comm, Teuchos::REDUCE_MAX, lcl_throw_exception, outArg(gbl_throw));                                          \
    TEUCHOS_TEST_FOR_EXCEPTION(gbl_throw, Exception,                                                                                \
                               msg << " Failure on at least one process, including process " << gbl_throw - 1 << "." << std::endl); \
  }
 
namespace Tpetra {

template <typename MapType, typename KeyArgType, typename ValueArgType>
typename MapType::iterator efficientAddOrUpdate(MapType& m,
                                                const KeyArgType& k,
                                                const ValueArgType& v) {
  typename MapType::iterator lb = m.lower_bound(k);
  if (lb != m.end() && !(m.key_comp()(k, lb->first))) {
    lb->second = v;
    return (lb);
  } else {
    typedef typename MapType::value_type MVT;
    return (m.insert(lb, MVT(k, v)));
  }
}

namespace SortDetails {

template <class IT1>
bool isAlreadySorted(const IT1& first, const IT1& last) {
  typedef typename std::iterator_traits<IT1>::difference_type DT;
  DT myit     = Teuchos::OrdinalTraits<DT>::one();
  const DT sz = last - first;
  for (; myit < sz; ++myit) {
    if (first[myit] < first[myit - 1]) {
      return false;
    }
  }
  return true;
}

template <class IT>
IT getPivot(const IT& first, const IT& last) {
  IT pivot(first + (last - first) / 2);
  if (*first <= *pivot && *(last - 1) <= *first)
    pivot = first;
  else if (*(last - 1) <= *pivot && *first <= *(last - 1))
    pivot = last - 1;
  return pivot;
}

template <class IT1, class IT2>
IT1 partition2(
    const IT1& first1,
    const IT1& last1,
    const IT2& first2,
    const IT2& last2,
    const IT1& pivot) {
  typename std::iterator_traits<IT1>::value_type piv(*pivot);
  std::swap(*pivot, *(last1 - 1));
  std::swap(first2[(pivot - first1)], *(last2 - 1));
  IT1 store1 = first1;
  for (IT1 it = first1; it != last1 - 1; ++it) {
    if (*it <= piv) {
      std::swap(*store1, *it);
      std::swap(first2[(store1 - first1)], first2[(it - first1)]);
      ++store1;
    }
  }
  std::swap(*(last1 - 1), *store1);
  std::swap(*(last2 - 1), first2[store1 - first1]);
  return store1;
}

template <class IT1, class IT2, class IT3>
IT1 partition3(
    const IT1& first1,
    const IT1& last1,
    const IT2& first2,
    const IT2& last2,
    const IT3& first3,
    const IT3& last3,
    const IT1& pivot) {
  typename std::iterator_traits<IT1>::value_type piv(*pivot);
  std::swap(*pivot, *(last1 - 1));
  std::swap(first2[(pivot - first1)], *(last2 - 1));
  std::swap(first3[(pivot - first1)], *(last3 - 1));
  IT1 store1 = first1;
  for (IT1 it = first1; it != last1 - 1; ++it) {
    if (*it <= piv) {
      std::swap(*store1, *it);
      std::swap(first2[(store1 - first1)], first2[(it - first1)]);
      std::swap(first3[(store1 - first1)], first3[(it - first1)]);
      ++store1;
    }
  }
  std::swap(*(last1 - 1), *store1);
  std::swap(*(last2 - 1), first2[store1 - first1]);
  std::swap(*(last3 - 1), first3[store1 - first1]);
  return store1;
}

template <class IT1, class IT2>
void quicksort2(
    const IT1& first1,
    const IT1& last1,
    const IT2& first2,
    const IT2& last2) {
  typedef typename std::iterator_traits<IT1>::difference_type DT;
  DT DT1 = Teuchos::OrdinalTraits<DT>::one();
  if (last1 - first1 > DT1) {
    IT1 pivot = getPivot(first1, last1);
    pivot     = partition2(first1, last1, first2, last2, pivot);
    quicksort2(first1, pivot, first2, first2 + (pivot - first1));
    quicksort2(pivot + 1, last1, first2 + (pivot - first1) + 1, last2);
  }
}

template <class IT1, class IT2, class IT3>
void quicksort3(
    const IT1& first1,
    const IT1& last1,
    const IT2& first2,
    const IT2& last2,
    const IT3& first3,
    const IT3& last3) {
  typedef typename std::iterator_traits<IT1>::difference_type DT;
  DT DT1 = Teuchos::OrdinalTraits<DT>::one();
  if (last1 - first1 > DT1) {
    IT1 pivot = getPivot(first1, last1);
    pivot     = partition3(first1, last1, first2, last2, first3, last3, pivot);
    quicksort3(first1, pivot, first2, first2 + (pivot - first1), first3, first3 + (pivot - first1));
    quicksort3(pivot + 1, last1, first2 + (pivot - first1) + 1, last2, first3 + (pivot - first1) + 1, last3);
  }
}

template <class IT1, class IT2, class IT3>
void sh_sort3(
    const IT1& first1,
    const IT1& last1,
    const IT2& first2,
    const IT2& /* last2 */,
    const IT3& first3,
    const IT3& /* last3 */) {
  typedef typename std::iterator_traits<IT1>::difference_type DT;
  DT n = last1 - first1;
  DT m = n / 2;
  DT z = Teuchos::OrdinalTraits<DT>::zero();
  while (m > z) {
    DT max = n - m;
    for (DT j = 0; j < max; j++) {
      for (DT k = j; k >= 0; k -= m) {
        if (first1[k + m] >= first1[k])
          break;
        std::swap(first1[k + m], first1[k]);
        std::swap(first2[k + m], first2[k]);
        std::swap(first3[k + m], first3[k]);
      }
    }
    m = m / 2;
  }
}

template <class IT1, class IT2>
void sh_sort2(
    const IT1& first1,
    const IT1& last1,
    const IT2& first2,
    const IT2& /* last2 */) {
  typedef typename std::iterator_traits<IT1>::difference_type DT;
  DT n = last1 - first1;
  DT m = n / 2;
  DT z = Teuchos::OrdinalTraits<DT>::zero();
  while (m > z) {
    DT max = n - m;
    for (DT j = 0; j < max; j++) {
      for (DT k = j; k >= 0; k -= m) {
        if (first1[k + m] >= first1[k])
          break;
        std::swap(first1[k + m], first1[k]);
        std::swap(first2[k + m], first2[k]);
      }
    }
    m = m / 2;
  }
}

template <typename IT1>
std::vector<typename std::iterator_traits<IT1>::difference_type> sort_indexes(IT1 first, IT1 last) {
  typedef typename std::iterator_traits<IT1>::difference_type DT;
 
  DT length = last - first;
  std::vector<DT> idx(length);
  std::iota(idx.begin(), idx.end(), 0);
 
  std::stable_sort(idx.begin(), idx.end(),
                   [&first](size_t i1, size_t i2) { return first[i1] < first[i2]; });
 
  return idx;
}

template <typename IT1, typename IT2>
void apply_permutation(
    IT1 first,
    IT1 last,
    IT2 indices) {
  typedef typename std::iterator_traits<IT1>::difference_type DT;
  typedef typename std::iterator_traits<IT1>::value_type T;
 
  DT n = last - first;
 
  Kokkos::View<T*, Kokkos::HostSpace> tmp(Kokkos::view_alloc(Kokkos::WithoutInitializing, "tmp"), n);
 
  Kokkos::parallel_for("apply_permutation_1",
                       Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>(0, n),
                       [=](const int j) {
                         tmp(j) = first[indices[j]];
                       });
  Kokkos::parallel_for("apply_permutation_2",
                       Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>(0, n),
                       [=](const int j) {
                         first[j] = tmp(j);
                       });
}

template <class IT1, class IT2>
void std_sort2(const IT1& first1,
               const IT1& last1,
               const IT2& first2,
               const IT2& last2) {
  auto indices = sort_indexes(first1, last1);
  apply_permutation(first1, last1, indices);
  apply_permutation(first2, last2, indices);
}

template <class IT1, class IT2, class IT3>
void std_sort3(
    const IT1& first1,
    const IT1& last1,
    const IT2& first2,
    const IT2& last2,
    const IT3& first3,
    const IT3& last3) {
  auto indices = sort_indexes(first1, last1);
  apply_permutation(first1, last1, indices);
  apply_permutation(first2, last2, indices);
  apply_permutation(first3, last3, indices);
}
 
}  // end namespace SortDetails

template <class IT1, class IT2>
void sort2(const IT1& first1, const IT1& last1, const IT2& first2, const bool stableSort = false) {
  // Quicksort uses best-case N log N time whether or not the input
  // data is sorted.  However, the common case in Tpetra is that the
  // input data are sorted, so we first check whether this is the
  // case before reverting to Quicksort.
  if (SortDetails::isAlreadySorted(first1, last1)) {
    return;
  }
  if (stableSort)
    SortDetails::std_sort2(first1, last1, first2, first2 + (last1 - first1));
  else
    SortDetails::sh_sort2(first1, last1, first2, first2 + (last1 - first1));
#ifdef HAVE_TPETRA_DEBUG
  if (!SortDetails::isAlreadySorted(first1, last1)) {
    std::cout << "Trouble: sort() did not sort !!" << std::endl;
    return;
  }
#endif
}

template <class View1, class View2>
void sort2(View1& view1, const size_t& size, View2& view2) {
  // NOTE: This assumes the view is host-accessible.
 
  // Wrap the views as rcps (this happens to preserve the reference counting, but that doesn't really matter here)
  Teuchos::ArrayRCP<typename View1::non_const_value_type> view1_rcp = Kokkos::Compat::persistingView(view1, 0, size);
  Teuchos::ArrayRCP<typename View2::non_const_value_type> view2_rcp = Kokkos::Compat::persistingView(view2, 0, size);
 
  sort2(view1_rcp.begin(), view1_rcp.end(), view2_rcp.begin());
}

template <class View>
void sort(View& view, const size_t& size) {
  // NOTE: This assumes the view is host-accessible.
 
  // Wrap the view as rcps (this happens to preserve the reference counting, but that doesn't really matter here)
  Teuchos::ArrayRCP<typename View::non_const_value_type> view_rcp = Kokkos::Compat::persistingView(view, 0, size);
 
  std::sort(view_rcp.begin(), view_rcp.end());
}

template <class View>
void reverse_sort(View& view, const size_t& size) {
  // NOTE: This assumes the view is host-accessible.
  // Wrap the view as rcps (this happens to preserve the reference counting, but that doesn't really matter here)
  Teuchos::ArrayRCP<typename View::non_const_value_type> view_rcp = Kokkos::Compat::persistingView(view, 0, size);
 
  std::sort(view_rcp.rbegin(), view_rcp.rend());
}

template <class IT1, class IT2, class IT3>
void sort3(const IT1& first1, const IT1& last1, const IT2& first2,
           const IT3& first3, const bool stableSort = false) {
  // Quicksort uses best-case N log N time whether or not the input
  // data is sorted.  However, the common case in Tpetra is that the
  // input data are sorted, so we first check whether this is the
  // case before reverting to Quicksort.
  if (SortDetails::isAlreadySorted(first1, last1)) {
    return;
  }
  if (stableSort)
    SortDetails::std_sort3(first1, last1, first2, first2 + (last1 - first1), first3,
                           first3 + (last1 - first1));
  else
    SortDetails::sh_sort3(first1, last1, first2, first2 + (last1 - first1), first3,
                          first3 + (last1 - first1));
#ifdef HAVE_TPETRA_DEBUG
  if (!SortDetails::isAlreadySorted(first1, last1)) {
    std::cout << " Trouble sort did not actually sort... !!!!!!" << std::endl;
    return;
  }
#endif
}

template <class IT1, class IT2>
void merge2(IT1& indResultOut, IT2& valResultOut,
            IT1 indBeg, IT1 indEnd,
            IT2 valBeg, IT2 /* valEnd */) {
  if (indBeg == indEnd) {
    indResultOut = indBeg;  // It's allowed for indResultOut to alias indEnd
    valResultOut = valBeg;  // It's allowed for valResultOut to alias valEnd
  } else {
    IT1 indResult = indBeg;
    IT2 valResult = valBeg;
    if (indBeg != indEnd) {
      ++indBeg;
      ++valBeg;
      while (indBeg != indEnd) {
        if (*indResult == *indBeg) {  // adjacent column indices equal
          *valResult += *valBeg;      // merge entries by adding their values together
        } else {                      // adjacent column indices not equal
          *(++indResult) = *indBeg;   // shift over the index
          *(++valResult) = *valBeg;   // shift over the value
        }
        ++indBeg;
        ++valBeg;
      }
      ++indResult;  // exclusive end of merged result
      ++valResult;  // exclusive end of merged result
 
      // mfh 24 Feb 2014: Setting these is technically correct, but
      // since the resulting variables aren't used after this point,
      // it may result in a build error.
      //
      // indEnd = indResult;
      // valEnd = valResult;
    }
    indResultOut = indResult;
    valResultOut = valResult;
  }
}

template <class IT1, class IT2, class BinaryFunction>
void merge2(IT1& indResultOut, IT2& valResultOut,
            IT1 indBeg, IT1 indEnd,
            IT2 valBeg, IT2 valEnd,
            BinaryFunction f) {
  if (indBeg == indEnd) {
    indResultOut = indBeg;  // It's allowed for indResultOut to alias indEnd
    valResultOut = valBeg;  // It's allowed for valResultOut to alias valEnd
  } else {
    IT1 indResult = indBeg;
    IT2 valResult = valBeg;
    if (indBeg != indEnd) {
      ++indBeg;
      ++valBeg;
      while (indBeg != indEnd) {
        if (*indResult == *indBeg) {            // adjacent column indices equal
          *valResult = f(*valResult, *valBeg);  // merge entries via values
        } else {                                // adjacent column indices not equal
          *(++indResult) = *indBeg;             // shift over the index
          *(++valResult) = *valBeg;             // shift over the value
        }
        ++indBeg;
        ++valBeg;
      }
      ++indResult;  // exclusive end of merged result
      ++valResult;  // exclusive end of merged result
      indEnd = indResult;
      valEnd = valResult;
    }
    indResultOut = indResult;
    valResultOut = valResult;
  }
}

template <class KeyInputIterType, class ValueInputIterType,
          class KeyOutputIterType, class ValueOutputIterType,
          class BinaryFunction>
void keyValueMerge(KeyInputIterType keyBeg1, KeyInputIterType keyEnd1,
                   ValueInputIterType valBeg1, ValueInputIterType valEnd1,
                   KeyInputIterType keyBeg2, KeyInputIterType keyEnd2,
                   ValueInputIterType valBeg2, ValueInputIterType valEnd2,
                   KeyOutputIterType keyOut, ValueOutputIterType valOut,
                   BinaryFunction f) {
  while (keyBeg1 != keyEnd1 && keyBeg2 != keyEnd2) {
    if (*keyBeg1 < *keyBeg2) {
      *keyOut++ = *keyBeg1++;
      *valOut++ = *valBeg1++;
    } else if (*keyBeg1 > *keyBeg2) {
      *keyOut++ = *keyBeg2++;
      *valOut++ = *valBeg2++;
    } else {  // *keyBeg1 == *keyBeg2
      *keyOut++ = *keyBeg1;
      *valOut++ = f(*valBeg1, *valBeg2);
      ++keyBeg1;
      ++valBeg1;
      ++keyBeg2;
      ++valBeg2;
    }
  }
  // At most one of the two sequences will be nonempty.
  std::copy(keyBeg1, keyEnd1, keyOut);
  std::copy(valBeg1, valEnd1, valOut);
  std::copy(keyBeg2, keyEnd2, keyOut);
  std::copy(valBeg2, valEnd2, valOut);
}
 
template <class KeyInputIterType>
size_t
keyMergeCount(KeyInputIterType keyBeg1, KeyInputIterType keyEnd1,
              KeyInputIterType keyBeg2, KeyInputIterType keyEnd2) {
  size_t count = 0;
  while (keyBeg1 != keyEnd1 && keyBeg2 != keyEnd2) {
    if (*keyBeg1 < *keyBeg2) {
      ++keyBeg1;
    } else if (*keyBeg1 > *keyBeg2) {
      ++keyBeg2;
    } else {  // *keyBeg1 == *keyBeg2
      ++keyBeg1;
      ++keyBeg2;
    }
    ++count;
  }
  // At most one of the two sequences will be nonempty.
  count += static_cast<size_t>(keyEnd1 - keyBeg1) +
           static_cast<size_t>(keyEnd2 - keyBeg2);
  return count;
}
 
namespace Details {
bool congruent(const Teuchos::Comm<int>& comm1,
               const Teuchos::Comm<int>& comm2);

template <class DualViewType>
Teuchos::ArrayView<typename DualViewType::t_dev::value_type>
getArrayViewFromDualView(const DualViewType& x) {
  static_assert(static_cast<int>(DualViewType::t_dev::rank) == 1,
                "The input DualView must have rank 1.");
  TEUCHOS_TEST_FOR_EXCEPTION(x.need_sync_host(), std::logic_error,
                             "The "
                             "input Kokkos::DualView was most recently modified on device, but this "
                             "function needs the host view of the data to be the most recently "
                             "modified.");
 
  auto x_host = x.view_host();
  typedef typename DualViewType::t_dev::value_type value_type;
  // mfh 15 Jan 2019: In debug mode, Teuchos::ArrayView's
  // constructor throws if the pointer is nonnull but the length
  // is nonpositive.
  const auto len = x_host.extent(0);
  return Teuchos::ArrayView<value_type>(len != 0 ? x_host.data() : nullptr,
                                        len);
}

template <class T, class DT>
Kokkos::DualView<T*, DT>
getDualViewCopyFromArrayView(const Teuchos::ArrayView<const T>& x_av,
                             const char label[],
                             const bool leaveOnHost) {
  using Kokkos::MemoryUnmanaged;
  typedef typename DT::memory_space DMS;
  typedef typename DT::execution_space execution_space;
  typedef Kokkos::HostSpace HMS;
 
  const size_t len = static_cast<size_t>(x_av.size());
  Kokkos::View<const T*, HMS, MemoryUnmanaged> x_in(x_av.getRawPtr(), len);
  Kokkos::DualView<T*, DT> x_out(label, len);
  if (leaveOnHost) {
    x_out.modify_host();
    // DEEP_COPY REVIEW - NOT TESTED FOR CUDA BUILD
    Kokkos::deep_copy(x_out.view_host(), x_in);
  } else {
    x_out.template modify<DMS>();
    // DEEP_COPY REVIEW - HOST-TO-DEVICE
    Kokkos::deep_copy(execution_space(), x_out.template view<DMS>(), x_in);
  }
  return x_out;
}

template <class DualViewType>
std::string dualViewStatusToString(const DualViewType& dv, const char name[]) {
  const auto host = dv.need_sync_device();
  const auto dev  = dv.need_sync_host();
 
  std::ostringstream os;
  os << name << ": {size: " << dv.extent(0)
     << ", sync: {host: " << host << ", dev: " << dev << "}";
  return os.str();
}

template <class ArrayType>
void verbosePrintArray(std::ostream& out,
                       const ArrayType& x,
                       const char name[],
                       const size_t maxNumToPrint) {
  out << name << ": [";
 
  const size_t numEnt(x.size());
  if (maxNumToPrint == 0) {
    if (numEnt != 0) {
      out << "...";
    }
  } else {
    const size_t numToPrint = numEnt > maxNumToPrint ? maxNumToPrint : numEnt;
    size_t k                = 0;
    for (; k < numToPrint; ++k) {
      out << x[k];
      if (k + size_t(1) < numToPrint) {
        out << ", ";
      }
    }
    if (k < numEnt) {
      out << ", ...";
    }
  }
  out << "]";
}

std::unique_ptr<std::string>
createPrefix(const int myRank,
             const char prefix[]);

std::unique_ptr<std::string>
createPrefix(const Teuchos::Comm<int>* comm,
             const char functionName[]);

std::unique_ptr<std::string>
createPrefix(const Teuchos::Comm<int>*,
             const char className[],
             const char methodName[]);
 
}  // namespace Details
}  // namespace Tpetra
 
#endif  // TPETRA_UTIL_HPP