Tpetra parallel linear algebra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 // @HEADER
41 
42 #ifndef TPETRA_UTIL_HPP
43 #define TPETRA_UTIL_HPP
44 
54 #include "Tpetra_ConfigDefs.hpp"
55 #include "Kokkos_DualView.hpp"
56 #include "Teuchos_Assert.hpp"
57 #include "Teuchos_CommHelpers.hpp"
58 #include "Teuchos_OrdinalTraits.hpp"
59 #include "Teuchos_TypeNameTraits.hpp"
60 #include "Teuchos_Utils.hpp"
61 #include <algorithm>
62 #include <iterator>
63 #include <sstream>
64 
65 #if defined(HAVE_TPETRA_THROW_EFFICIENCY_WARNINGS) || defined(HAVE_TPETRA_PRINT_EFFICIENCY_WARNINGS)
66 #define TPETRA_EFFICIENCY_WARNING(throw_exception_test,Exception,msg) \
100 { \
101  const bool tpetraEfficiencyWarningTest = (throw_exception_test); \
102  if (tpetraEfficiencyWarningTest) { \
103  std::ostringstream errStream; \
104  errStream << Teuchos::typeName(*this) << ":" << std::endl; \
105  errStream << "Efficiency warning: " << #throw_exception_test << std::endl; \
106  errStream << msg; \
107  std::string err = errStream.str(); \
108  if (TPETRA_PRINTS_EFFICIENCY_WARNINGS && tpetraEfficiencyWarningTest) { \
109  std::cerr << err << std::endl; \
110  } \
111  TEUCHOS_TEST_FOR_EXCEPTION(TPETRA_THROWS_EFFICIENCY_WARNINGS && tpetraEfficiencyWarningTest, Exception, err); \
112  } \
113 }
114 #else
115 #define TPETRA_EFFICIENCY_WARNING(throw_exception_test,Exception,msg)
149 #endif
150 
151 // handle an abuse warning, according to HAVE_TPETRA_THROW_ABUSE_WARNINGS and HAVE_TPETRA_PRINT_ABUSE_WARNINGS
152 #if defined(HAVE_TPETRA_THROW_ABUSE_WARNINGS) || defined(HAVE_TPETRA_PRINT_ABUSE_WARNINGS)
153 #define TPETRA_ABUSE_WARNING(throw_exception_test,Exception,msg) \
155 { \
156  std::ostringstream errStream; \
157  errStream << Teuchos::typeName(*this) << msg; \
158  std::string err = errStream.str(); \
159  if (TPETRA_PRINTS_ABUSE_WARNINGS && (throw_exception_test)) { \
160  std::cerr << err << std::endl; \
161  } \
162  TEUCHOS_TEST_FOR_EXCEPTION(TPETRA_THROWS_ABUSE_WARNINGS && (throw_exception_test), Exception, err); \
163 }
164 #else
165 #define TPETRA_ABUSE_WARNING(throw_exception_test,Exception,msg)
167 #endif
168 
198 #define SHARED_TEST_FOR_EXCEPTION(throw_exception_test,Exception,msg,comm) \
199 { \
200  using Teuchos::outArg; \
201  const int lcl_throw_exception = (throw_exception_test) ? Teuchos::rank(comm)+1 : 0; \
202  int gbl_throw; \
203  Teuchos::reduceAll(comm,Teuchos::REDUCE_MAX,lcl_throw_exception,outArg(gbl_throw)); \
204  TEUCHOS_TEST_FOR_EXCEPTION(gbl_throw,Exception, \
205  msg << " Failure on at least one process, including process " << gbl_throw-1 << "." << std::endl); \
206 }
207 
208 #ifdef HAVE_TEUCHOS_DEBUG
209 #define SWITCHED_TEST_FOR_EXCEPTION(throw_exception_test,Exception,msg,comm) \
211 { \
212  SHARED_TEST_FOR_EXCEPTION(throw_exception_test,Exception,msg,comm); \
213 }
214 #else
215 #define SWITCHED_TEST_FOR_EXCEPTION(throw_exception_test,Exception,msg,comm) \
217 { \
218  TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test,Exception,msg); \
219 }
220 #endif
221 
222 namespace Tpetra {
223 
235  template<typename MapType, typename KeyArgType, typename ValueArgType>
236  typename MapType::iterator efficientAddOrUpdate(MapType& m,
237  const KeyArgType & k,
238  const ValueArgType & v)
239  {
240  typename MapType::iterator lb = m.lower_bound(k);
241  if(lb != m.end() && !(m.key_comp()(k, lb->first))) {
242  lb->second = v;
243  return(lb);
244  }
245  else {
246  typedef typename MapType::value_type MVT;
247  return(m.insert(lb, MVT(k, v)));
248  }
249  }
250 
258  namespace SortDetails {
259 
268  template<class IT1>
269  bool isAlreadySorted(const IT1& first, const IT1& last){
270  typedef typename std::iterator_traits<IT1>::difference_type DT;
271  DT myit = Teuchos::OrdinalTraits<DT>::one();
272  const DT sz = last - first;
273  for(;myit < sz; ++myit){
274  if(first[myit] < first[myit-1]){
275  return false;
276  }
277  }
278  return true;
279  }
280 
290  template<class IT>
291  IT getPivot(const IT& first, const IT& last){
292  IT pivot(first+(last-first)/2);
293  if(*first<=*pivot && *(last-1)<=*first) pivot=first;
294  else if(*(last-1)<=*pivot && *first<= *(last-1)) pivot = last-1;
295  return pivot;
296  }
297 
312  template<class IT1, class IT2>
313  IT1 partition2(
314  const IT1& first1,
315  const IT1& last1,
316  const IT2& first2,
317  const IT2& last2,
318  const IT1& pivot)
319  {
320  typename std::iterator_traits<IT1>::value_type piv(*pivot);
321  std::swap(*pivot, *(last1-1));
322  std::swap(first2[(pivot-first1)], *(last2-1));
323  IT1 store1=first1;
324  for(IT1 it=first1; it!=last1-1; ++it){
325  if(*it<=piv){
326  std::swap(*store1, *it);
327  std::swap(first2[(store1-first1)], first2[(it-first1)]);
328  ++store1;
329  }
330  }
331  std::swap(*(last1-1), *store1);
332  std::swap(*(last2-1), first2[store1-first1]);
333  return store1;
334  }
335 
352  template<class IT1, class IT2, class IT3>
353  IT1 partition3(
354  const IT1& first1,
355  const IT1& last1,
356  const IT2& first2,
357  const IT2& last2,
358  const IT3& first3,
359  const IT3& last3,
360  const IT1& pivot)
361  {
362  typename std::iterator_traits<IT1>::value_type piv(*pivot);
363  std::swap(*pivot, *(last1-1));
364  std::swap(first2[(pivot-first1)], *(last2-1));
365  std::swap(first3[(pivot-first1)], *(last3-1));
366  IT1 store1=first1;
367  for(IT1 it=first1; it!=last1-1; ++it){
368  if(*it<=piv){
369  std::swap(*store1, *it);
370  std::swap(first2[(store1-first1)], first2[(it-first1)]);
371  std::swap(first3[(store1-first1)], first3[(it-first1)]);
372  ++store1;
373  }
374  }
375  std::swap(*(last1-1), *store1);
376  std::swap(*(last2-1), first2[store1-first1]);
377  std::swap(*(last3-1), first3[store1-first1]);
378  return store1;
379  }
380 
391  template<class IT1, class IT2>
392  void quicksort2(
393  const IT1& first1,
394  const IT1& last1,
395  const IT2& first2,
396  const IT2& last2)
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 = partition2(first1, last1, first2, last2, pivot);
403  quicksort2(first1, pivot, first2, first2+(pivot-first1));
404  quicksort2(pivot+1, last1, first2+(pivot-first1)+1, last2);
405  }
406  }
407 
420  template<class IT1, class IT2, class IT3>
421  void quicksort3(
422  const IT1& first1,
423  const IT1& last1,
424  const IT2& first2,
425  const IT2& last2,
426  const IT3& first3,
427  const IT3& last3)
428  {
429  typedef typename std::iterator_traits<IT1>::difference_type DT;
430  DT DT1 = Teuchos::OrdinalTraits<DT>::one();
431  if(last1-first1 > DT1){
432  IT1 pivot = getPivot(first1, last1);
433  pivot = partition3(first1, last1, first2, last2, first3, last3, pivot);
434  quicksort3(first1, pivot, first2, first2+(pivot-first1), first3, first3+(pivot-first1));
435  quicksort3(pivot+1, last1, first2+(pivot-first1)+1, last2, first3+(pivot-first1)+1, last3);
436  }
437  }
438 
445  template<class IT1, class IT2, class IT3>
446  void sh_sort3(
447  const IT1& first1,
448  const IT1& last1,
449  const IT2& first2,
450  const IT2& /* last2 */,
451  const IT3& first3,
452  const IT3& /* last3 */)
453  {
454  typedef typename std::iterator_traits<IT1>::difference_type DT;
455  DT n = last1 - first1;
456  DT m = n / 2;
457  DT z = Teuchos::OrdinalTraits<DT>::zero();
458  while (m > z)
459  {
460  DT max = n - m;
461  for (DT j = 0; j < max; j++)
462  {
463  for (DT k = j; k >= 0; k-=m)
464  {
465  if (first1[k+m] >= first1[k])
466  break;
467  std::swap(first1[k+m], first1[k]);
468  std::swap(first2[k+m], first2[k]);
469  std::swap(first3[k+m], first3[k]);
470  }
471  }
472  m = m/2;
473  }
474  }
475 
482  template<class IT1, class IT2>
483  void sh_sort2(
484  const IT1& first1,
485  const IT1& last1,
486  const IT2& first2,
487  const IT2& /* last2 */)
488  {
489  typedef typename std::iterator_traits<IT1>::difference_type DT;
490  DT n = last1 - first1;
491  DT m = n / 2;
492  DT z = Teuchos::OrdinalTraits<DT>::zero();
493  while (m > z)
494  {
495  DT max = n - m;
496  for (DT j = 0; j < max; j++)
497  {
498  for (DT k = j; k >= 0; k-=m)
499  {
500  if (first1[k+m] >= first1[k])
501  break;
502  std::swap(first1[k+m], first1[k]);
503  std::swap(first2[k+m], first2[k]);
504  }
505  }
506  m = m/2;
507  }
508  }
509 
510  } //end namespace SortDetails
511 
512 
531  template<class IT1, class IT2>
532  void sort2(const IT1 &first1, const IT1 &last1, const IT2 &first2) {
533  // Quicksort uses best-case N log N time whether or not the input
534  // data is sorted. However, the common case in Tpetra is that the
535  // input data are sorted, so we first check whether this is the
536  // case before reverting to Quicksort.
537  if(SortDetails::isAlreadySorted(first1, last1)){
538  return;
539  }
540  SortDetails::sh_sort2(first1, last1, first2, first2+(last1-first1));
541 #ifdef HAVE_TPETRA_DEBUG
542  if(!SortDetails::isAlreadySorted(first1, last1)){
543  std::cout << "Trouble: sort() did not sort !!" << std::endl;
544  return;
545  }
546 #endif
547  }
548 
549 
565  template<class IT1, class IT2, class IT3>
566  void sort3(const IT1 &first1, const IT1 &last1, const IT2 &first2,
567  const IT3 &first3)
568  {
569  // Quicksort uses best-case N log N time whether or not the input
570  // data is sorted. However, the common case in Tpetra is that the
571  // input data are sorted, so we first check whether this is the
572  // case before reverting to Quicksort.
573  if(SortDetails::isAlreadySorted(first1, last1)){
574  return;
575  }
576  SortDetails::sh_sort3(first1, last1, first2, first2+(last1-first1), first3,
577  first3+(last1-first1));
578 
579 #ifdef HAVE_TPETRA_DEBUG
580  if(!SortDetails::isAlreadySorted(first1, last1)){
581  std::cout << " Trouble sort did not actually sort... !!!!!!" <<
582  std::endl;
583  return;
584  }
585 #endif
586  }
587 
631  template<class IT1, class IT2>
632  void
633  merge2 (IT1& indResultOut, IT2& valResultOut,
634  IT1 indBeg, IT1 indEnd,
635  IT2 valBeg, IT2 /* valEnd */)
636  {
637  if (indBeg == indEnd) {
638  indResultOut = indBeg; // It's allowed for indResultOut to alias indEnd
639  valResultOut = valBeg; // It's allowed for valResultOut to alias valEnd
640  }
641  else {
642  IT1 indResult = indBeg;
643  IT2 valResult = valBeg;
644  if (indBeg != indEnd) {
645  ++indBeg;
646  ++valBeg;
647  while (indBeg != indEnd) {
648  if (*indResult == *indBeg) { // adjacent column indices equal
649  *valResult += *valBeg; // merge entries by adding their values together
650  } else { // adjacent column indices not equal
651  *(++indResult) = *indBeg; // shift over the index
652  *(++valResult) = *valBeg; // shift over the value
653  }
654  ++indBeg;
655  ++valBeg;
656  }
657  ++indResult; // exclusive end of merged result
658  ++valResult; // exclusive end of merged result
659 
660  // mfh 24 Feb 2014: Setting these is technically correct, but
661  // since the resulting variables aren't used after this point,
662  // it may result in a build error.
663  //
664  // indEnd = indResult;
665  // valEnd = valResult;
666  }
667  indResultOut = indResult;
668  valResultOut = valResult;
669  }
670  }
671 
720  template<class IT1, class IT2, class BinaryFunction>
721  void
722  merge2 (IT1& indResultOut, IT2& valResultOut,
723  IT1 indBeg, IT1 indEnd,
724  IT2 valBeg, IT2 valEnd,
725  BinaryFunction f)
726  {
727  if (indBeg == indEnd) {
728  indResultOut = indBeg; // It's allowed for indResultOut to alias indEnd
729  valResultOut = valBeg; // It's allowed for valResultOut to alias valEnd
730  }
731  else {
732  IT1 indResult = indBeg;
733  IT2 valResult = valBeg;
734  if (indBeg != indEnd) {
735  ++indBeg;
736  ++valBeg;
737  while (indBeg != indEnd) {
738  if (*indResult == *indBeg) { // adjacent column indices equal
739  *valResult = f (*valResult, *valBeg); // merge entries via values
740  } else { // adjacent column indices not equal
741  *(++indResult) = *indBeg; // shift over the index
742  *(++valResult) = *valBeg; // shift over the value
743  }
744  ++indBeg;
745  ++valBeg;
746  }
747  ++indResult; // exclusive end of merged result
748  ++valResult; // exclusive end of merged result
749  indEnd = indResult;
750  valEnd = valResult;
751  }
752  indResultOut = indResult;
753  valResultOut = valResult;
754  }
755  }
756 
757 
786  template<class KeyInputIterType, class ValueInputIterType,
787  class KeyOutputIterType, class ValueOutputIterType,
788  class BinaryFunction>
789  void
790  keyValueMerge (KeyInputIterType keyBeg1, KeyInputIterType keyEnd1,
791  ValueInputIterType valBeg1, ValueInputIterType valEnd1,
792  KeyInputIterType keyBeg2, KeyInputIterType keyEnd2,
793  ValueInputIterType valBeg2, ValueInputIterType valEnd2,
794  KeyOutputIterType keyOut, ValueOutputIterType valOut,
795  BinaryFunction f)
796  {
797  while (keyBeg1 != keyEnd1 && keyBeg2 != keyEnd2) {
798  if (*keyBeg1 < *keyBeg2) {
799  *keyOut++ = *keyBeg1++;
800  *valOut++ = *valBeg1++;
801  } else if (*keyBeg1 > *keyBeg2) {
802  *keyOut++ = *keyBeg2++;
803  *valOut++ = *valBeg2++;
804  } else { // *keyBeg1 == *keyBeg2
805  *keyOut++ = *keyBeg1;
806  *valOut++ = f (*valBeg1, *valBeg2);
807  ++keyBeg1;
808  ++valBeg1;
809  ++keyBeg2;
810  ++valBeg2;
811  }
812  }
813  // At most one of the two sequences will be nonempty.
814  std::copy (keyBeg1, keyEnd1, keyOut);
815  std::copy (valBeg1, valEnd1, valOut);
816  std::copy (keyBeg2, keyEnd2, keyOut);
817  std::copy (valBeg2, valEnd2, valOut);
818  }
819 
820  template<class KeyInputIterType>
821  size_t
822  keyMergeCount (KeyInputIterType keyBeg1, KeyInputIterType keyEnd1,
823  KeyInputIterType keyBeg2, KeyInputIterType keyEnd2)
824  {
825  size_t count = 0;
826  while (keyBeg1 != keyEnd1 && keyBeg2 != keyEnd2) {
827  if (*keyBeg1 < *keyBeg2) {
828  ++keyBeg1;
829  } else if (*keyBeg1 > *keyBeg2) {
830  ++keyBeg2;
831  } else { // *keyBeg1 == *keyBeg2
832  ++keyBeg1;
833  ++keyBeg2;
834  }
835  ++count;
836  }
837  // At most one of the two sequences will be nonempty.
838  count += static_cast<size_t> (keyEnd1 - keyBeg1) +
839  static_cast<size_t> (keyEnd2 - keyBeg2);
840  return count;
841  }
842 
843  namespace Details {
863  bool
864  congruent (const Teuchos::Comm<int>& comm1,
865  const Teuchos::Comm<int>& comm2);
866 
876  template<class DualViewType>
877  Teuchos::ArrayView<typename DualViewType::t_dev::value_type>
878  getArrayViewFromDualView (const DualViewType& x)
879  {
880  static_assert (static_cast<int> (DualViewType::t_dev::rank) == 1,
881  "The input DualView must have rank 1.");
882  TEUCHOS_TEST_FOR_EXCEPTION
883  (x.need_sync_host (), std::logic_error, "The "
884  "input Kokkos::DualView was most recently modified on device, but this "
885  "function needs the host view of the data to be the most recently "
886  "modified.");
887 
888  auto x_host = x.view_host ();
889  typedef typename DualViewType::t_dev::value_type value_type;
890  // mfh 15 Jan 2019: In debug mode, Teuchos::ArrayView's
891  // constructor throws if the pointer is nonnull but the length
892  // is nonpositive.
893  const auto len = x_host.extent (0);
894  return Teuchos::ArrayView<value_type> (len != 0 ? x_host.data () : nullptr,
895  len);
896  }
897 
914  template<class T, class DT>
915  Kokkos::DualView<T*, DT>
916  getDualViewCopyFromArrayView (const Teuchos::ArrayView<const T>& x_av,
917  const char label[],
918  const bool leaveOnHost)
919  {
920  using Kokkos::MemoryUnmanaged;
921  typedef typename DT::memory_space DMS;
922  typedef Kokkos::HostSpace HMS;
923 
924  const size_t len = static_cast<size_t> (x_av.size ());
925  Kokkos::View<const T*, HMS, MemoryUnmanaged> x_in (x_av.getRawPtr (), len);
926  Kokkos::DualView<T*, DT> x_out (label, len);
927  if (leaveOnHost) {
928  x_out.modify_host ();
929  Kokkos::deep_copy (x_out.view_host (), x_in);
930  }
931  else {
932  x_out.template modify<DMS> ();
933  Kokkos::deep_copy (x_out.template view<DMS> (), x_in);
934  }
935  return x_out;
936  }
937 
945  template<class DualViewType>
946  std::string dualViewStatusToString (const DualViewType& dv, const char name[])
947  {
948  const auto host = dv.need_sync_device();
949  const auto dev = dv.need_sync_host();
950 
951  std::ostringstream os;
952  os << name << ": {size: " << dv.extent (0)
953  << ", sync: {host: " << host << ", dev: " << dev << "}";
954  return os.str ();
955  }
956 
957  } // namespace Details
958 } // namespace Tpetra
959 
960 #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:65
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 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.