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