Tpetra parallel linear algebra  Version of the Day
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 
71 #include "Tpetra_ConfigDefs.hpp" // for map, vector, string, and iostream
72 #include <iterator>
73 #include <algorithm>
74 #include <Teuchos_Utils.hpp>
75 #include <Teuchos_Assert.hpp>
76 #include <sstream>
77 
78 #if defined(HAVE_TPETRA_THROW_EFFICIENCY_WARNINGS) || defined(HAVE_TPETRA_PRINT_EFFICIENCY_WARNINGS)
79 #define TPETRA_EFFICIENCY_WARNING(throw_exception_test,Exception,msg) \
113 { \
114  const bool tpetraEfficiencyWarningTest = (throw_exception_test); \
115  if (tpetraEfficiencyWarningTest) { \
116  std::ostringstream errStream; \
117  errStream << Teuchos::typeName(*this) << ":" << std::endl; \
118  errStream << "Efficiency warning: " << #throw_exception_test << std::endl; \
119  errStream << msg; \
120  std::string err = errStream.str(); \
121  if (TPETRA_PRINTS_EFFICIENCY_WARNINGS && tpetraEfficiencyWarningTest) { \
122  std::cerr << err << std::endl; \
123  } \
124  TEUCHOS_TEST_FOR_EXCEPTION(TPETRA_THROWS_EFFICIENCY_WARNINGS && tpetraEfficiencyWarningTest, Exception, err); \
125  } \
126 }
127 #else
128 #define TPETRA_EFFICIENCY_WARNING(throw_exception_test,Exception,msg)
162 #endif
163 
164 // handle an abuse warning, according to HAVE_TPETRA_THROW_ABUSE_WARNINGS and HAVE_TPETRA_PRINT_ABUSE_WARNINGS
165 #if defined(HAVE_TPETRA_THROW_ABUSE_WARNINGS) || defined(HAVE_TPETRA_PRINT_ABUSE_WARNINGS)
166 #define TPETRA_ABUSE_WARNING(throw_exception_test,Exception,msg) \
168 { \
169  std::ostringstream errStream; \
170  errStream << Teuchos::typeName(*this) << msg; \
171  std::string err = errStream.str(); \
172  if (TPETRA_PRINTS_ABUSE_WARNINGS && (throw_exception_test)) { \
173  std::cerr << err << std::endl; \
174  } \
175  TEUCHOS_TEST_FOR_EXCEPTION(TPETRA_THROWS_ABUSE_WARNINGS && (throw_exception_test), Exception, err); \
176 }
177 #else
178 #define TPETRA_ABUSE_WARNING(throw_exception_test,Exception,msg)
180 #endif
181 
211 #define SHARED_TEST_FOR_EXCEPTION(throw_exception_test,Exception,msg,comm) \
212 { \
213  using Teuchos::outArg; \
214  const int lcl_throw_exception = (throw_exception_test) ? Teuchos::rank(comm)+1 : 0; \
215  int gbl_throw; \
216  Teuchos::reduceAll(comm,Teuchos::REDUCE_MAX,lcl_throw_exception,outArg(gbl_throw)); \
217  TEUCHOS_TEST_FOR_EXCEPTION(gbl_throw,Exception, \
218  msg << " Failure on at least one process, including process " << gbl_throw-1 << "." << std::endl); \
219 }
220 
221 #ifdef HAVE_TEUCHOS_DEBUG
222 #define SWITCHED_TEST_FOR_EXCEPTION(throw_exception_test,Exception,msg,comm) \
224 { \
225  SHARED_TEST_FOR_EXCEPTION(throw_exception_test,Exception,msg,comm); \
226 }
227 #else
228 #define SWITCHED_TEST_FOR_EXCEPTION(throw_exception_test,Exception,msg,comm) \
230 { \
231  TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test,Exception,msg); \
232 }
233 #endif
234 
235 namespace Tpetra {
236 
248  template<typename MapType, typename KeyArgType, typename ValueArgType>
249  typename MapType::iterator efficientAddOrUpdate(MapType& m,
250  const KeyArgType & k,
251  const ValueArgType & v)
252  {
253  typename MapType::iterator lb = m.lower_bound(k);
254  if(lb != m.end() && !(m.key_comp()(k, lb->first))) {
255  lb->second = v;
256  return(lb);
257  }
258  else {
259  typedef typename MapType::value_type MVT;
260  return(m.insert(lb, MVT(k, v)));
261  }
262  }
263 
271  namespace SortDetails {
272 
281  template<class IT1>
282  bool isAlreadySorted(const IT1& first, const IT1& last){
283  typedef typename std::iterator_traits<IT1>::difference_type DT;
284  DT myit =OrdinalTraits<DT>::one();
285  const DT sz = last - first;
286  for(;myit < sz; ++myit){
287  if(first[myit] < first[myit-1]){
288  return false;
289  }
290  }
291  return true;
292  }
293 
303  template<class IT>
304  IT getPivot(const IT& first, const IT& last){
305  IT pivot(first+(last-first)/2);
306  if(*first<=*pivot && *(last-1)<=*first) pivot=first;
307  else if(*(last-1)<=*pivot && *first<= *(last-1)) pivot = last-1;
308  return pivot;
309  }
310 
325  template<class IT1, class IT2>
327  const IT1& first1,
328  const IT1& last1,
329  const IT2& first2,
330  const IT2& last2,
331  const IT1& pivot)
332  {
333  typename std::iterator_traits<IT1>::value_type piv(*pivot);
334  std::swap(*pivot, *(last1-1));
335  std::swap(first2[(pivot-first1)], *(last2-1));
336  IT1 store1=first1;
337  for(IT1 it=first1; it!=last1-1; ++it){
338  if(*it<=piv){
339  std::swap(*store1, *it);
340  std::swap(first2[(store1-first1)], first2[(it-first1)]);
341  ++store1;
342  }
343  }
344  std::swap(*(last1-1), *store1);
345  std::swap(*(last2-1), first2[store1-first1]);
346  return store1;
347  }
348 
365  template<class IT1, class IT2, class IT3>
367  const IT1& first1,
368  const IT1& last1,
369  const IT2& first2,
370  const IT2& last2,
371  const IT3& first3,
372  const IT3& last3,
373  const IT1& pivot)
374  {
375  typename std::iterator_traits<IT1>::value_type piv(*pivot);
376  std::swap(*pivot, *(last1-1));
377  std::swap(first2[(pivot-first1)], *(last2-1));
378  std::swap(first3[(pivot-first1)], *(last3-1));
379  IT1 store1=first1;
380  for(IT1 it=first1; it!=last1-1; ++it){
381  if(*it<=piv){
382  std::swap(*store1, *it);
383  std::swap(first2[(store1-first1)], first2[(it-first1)]);
384  std::swap(first3[(store1-first1)], first3[(it-first1)]);
385  ++store1;
386  }
387  }
388  std::swap(*(last1-1), *store1);
389  std::swap(*(last2-1), first2[store1-first1]);
390  std::swap(*(last3-1), first3[store1-first1]);
391  return store1;
392  }
393 
404  template<class IT1, class IT2>
406  const IT1& first1,
407  const IT1& last1,
408  const IT2& first2,
409  const IT2& last2)
410  {
411  typedef typename std::iterator_traits<IT1>::difference_type DT;
412  DT DT1 = OrdinalTraits<DT>::one();
413  if(last1-first1 > DT1){
414  IT1 pivot = getPivot(first1, last1);
415  pivot = partition2(first1, last1, first2, last2, pivot);
416  quicksort2(first1, pivot, first2, first2+(pivot-first1));
417  quicksort2(pivot+1, last1, first2+(pivot-first1)+1, last2);
418  }
419  }
420 
433  template<class IT1, class IT2, class IT3>
435  const IT1& first1,
436  const IT1& last1,
437  const IT2& first2,
438  const IT2& last2,
439  const IT3& first3,
440  const IT3& last3)
441  {
442  typedef typename std::iterator_traits<IT1>::difference_type DT;
443  DT DT1 = OrdinalTraits<DT>::one();
444  if(last1-first1 > DT1){
445  IT1 pivot = getPivot(first1, last1);
446  pivot = partition3(first1, last1, first2, last2, first3, last3, pivot);
447  quicksort3(first1, pivot, first2, first2+(pivot-first1), first3, first3+(pivot-first1));
448  quicksort3(pivot+1, last1, first2+(pivot-first1)+1, last2, first3+(pivot-first1)+1, last3);
449  }
450  }
451 
458  template<class IT1, class IT2, class IT3>
459  void sh_sort3(
460  const IT1& first1,
461  const IT1& last1,
462  const IT2& first2,
463  const IT2& last2,
464  const IT3& first3,
465  const IT3& last3)
466  {
467  typedef typename std::iterator_traits<IT1>::difference_type DT;
468  DT n = last1 - first1;
469  DT m = n / 2;
470  DT z = OrdinalTraits<DT>::zero();
471  while (m > z)
472  {
473  DT max = n - m;
474  for (DT j = 0; j < max; j++)
475  {
476  for (DT k = j; k >= 0; k-=m)
477  {
478  if (first1[k+m] >= first1[k])
479  break;
480  std::swap(first1[k+m], first1[k]);
481  std::swap(first2[k+m], first2[k]);
482  std::swap(first3[k+m], first3[k]);
483  }
484  }
485  m = m/2;
486  }
487  }
488 
495  template<class IT1, class IT2>
496  void sh_sort2(
497  const IT1& first1,
498  const IT1& last1,
499  const IT2& first2,
500  const IT2& last2)
501  {
502  typedef typename std::iterator_traits<IT1>::difference_type DT;
503  DT n = last1 - first1;
504  DT m = n / 2;
505  DT z = OrdinalTraits<DT>::zero();
506  while (m > z)
507  {
508  DT max = n - m;
509  for (DT j = 0; j < max; j++)
510  {
511  for (DT k = j; k >= 0; k-=m)
512  {
513  if (first1[k+m] >= first1[k])
514  break;
515  std::swap(first1[k+m], first1[k]);
516  std::swap(first2[k+m], first2[k]);
517  }
518  }
519  m = m/2;
520  }
521  }
522 
523  } //end namespace SortDetails
524 
525 
544  template<class IT1, class IT2>
545  void sort2(const IT1 &first1, const IT1 &last1, const IT2 &first2) {
546  // Quicksort uses best-case N log N time whether or not the input
547  // data is sorted. However, the common case in Tpetra is that the
548  // input data are sorted, so we first check whether this is the
549  // case before reverting to Quicksort.
550  if(SortDetails::isAlreadySorted(first1, last1)){
551  return;
552  }
553  SortDetails::sh_sort2(first1, last1, first2, first2+(last1-first1));
554 #ifdef HAVE_TPETRA_DEBUG
555  if(!SortDetails::isAlreadySorted(first1, last1)){
556  std::cout << "Trouble: sort() did not sort !!" << std::endl;
557  return;
558  }
559 #endif
560  }
561 
562 
578  template<class IT1, class IT2, class IT3>
579  void sort3(const IT1 &first1, const IT1 &last1, const IT2 &first2,
580  const IT3 &first3)
581  {
582  // Quicksort uses best-case N log N time whether or not the input
583  // data is sorted. However, the common case in Tpetra is that the
584  // input data are sorted, so we first check whether this is the
585  // case before reverting to Quicksort.
586  if(SortDetails::isAlreadySorted(first1, last1)){
587  return;
588  }
589  SortDetails::sh_sort3(first1, last1, first2, first2+(last1-first1), first3,
590  first3+(last1-first1));
591 
592 #ifdef HAVE_TPETRA_DEBUG
593  if(!SortDetails::isAlreadySorted(first1, last1)){
594  std::cout << " Trouble sort did not actually sort... !!!!!!" <<
595  std::endl;
596  return;
597  }
598 #endif
599  }
600 
644  template<class IT1, class IT2>
645  void
646  merge2 (IT1& indResultOut, IT2& valResultOut,
647  IT1 indBeg, IT1 indEnd,
648  IT2 valBeg, IT2 valEnd)
649  {
650  if (indBeg == indEnd) {
651  indResultOut = indBeg; // It's allowed for indResultOut to alias indEnd
652  valResultOut = valBeg; // It's allowed for valResultOut to alias valEnd
653  }
654  else {
655  IT1 indResult = indBeg;
656  IT2 valResult = valBeg;
657  if (indBeg != indEnd) {
658  ++indBeg;
659  ++valBeg;
660  while (indBeg != indEnd) {
661  if (*indResult == *indBeg) { // adjacent column indices equal
662  *valResult += *valBeg; // merge entries by adding their values together
663  } else { // adjacent column indices not equal
664  *(++indResult) = *indBeg; // shift over the index
665  *(++valResult) = *valBeg; // shift over the value
666  }
667  ++indBeg;
668  ++valBeg;
669  }
670  ++indResult; // exclusive end of merged result
671  ++valResult; // exclusive end of merged result
672 
673  // mfh 24 Feb 2014: Setting these is technically correct, but
674  // since the resulting variables aren't used after this point,
675  // it may result in a build error.
676  //
677  // indEnd = indResult;
678  // valEnd = valResult;
679  }
680  indResultOut = indResult;
681  valResultOut = valResult;
682  }
683  }
684 
733  template<class IT1, class IT2, class BinaryFunction>
734  void
735  merge2 (IT1& indResultOut, IT2& valResultOut,
736  IT1 indBeg, IT1 indEnd,
737  IT2 valBeg, IT2 valEnd,
738  BinaryFunction f)
739  {
740  if (indBeg == indEnd) {
741  indResultOut = indBeg; // It's allowed for indResultOut to alias indEnd
742  valResultOut = valBeg; // It's allowed for valResultOut to alias valEnd
743  }
744  else {
745  IT1 indResult = indBeg;
746  IT2 valResult = valBeg;
747  if (indBeg != indEnd) {
748  ++indBeg;
749  ++valBeg;
750  while (indBeg != indEnd) {
751  if (*indResult == *indBeg) { // adjacent column indices equal
752  *valResult = f (*valResult, *valBeg); // merge entries via values
753  } else { // adjacent column indices not equal
754  *(++indResult) = *indBeg; // shift over the index
755  *(++valResult) = *valBeg; // shift over the value
756  }
757  ++indBeg;
758  ++valBeg;
759  }
760  ++indResult; // exclusive end of merged result
761  ++valResult; // exclusive end of merged result
762  indEnd = indResult;
763  valEnd = valResult;
764  }
765  indResultOut = indResult;
766  valResultOut = valResult;
767  }
768  }
769 
770 
799  template<class KeyInputIterType, class ValueInputIterType,
800  class KeyOutputIterType, class ValueOutputIterType,
801  class BinaryFunction>
802  void
803  keyValueMerge (KeyInputIterType keyBeg1, KeyInputIterType keyEnd1,
804  ValueInputIterType valBeg1, ValueInputIterType valEnd1,
805  KeyInputIterType keyBeg2, KeyInputIterType keyEnd2,
806  ValueInputIterType valBeg2, ValueInputIterType valEnd2,
807  KeyOutputIterType keyOut, ValueOutputIterType valOut,
808  BinaryFunction f)
809  {
810  while (keyBeg1 != keyEnd1 && keyBeg2 != keyEnd2) {
811  if (*keyBeg1 < *keyBeg2) {
812  *keyOut++ = *keyBeg1++;
813  *valOut++ = *valBeg1++;
814  } else if (*keyBeg1 > *keyBeg2) {
815  *keyOut++ = *keyBeg2++;
816  *valOut++ = *valBeg2++;
817  } else { // *keyBeg1 == *keyBeg2
818  *keyOut++ = *keyBeg1;
819  *valOut++ = f (*valBeg1, *valBeg2);
820  ++keyBeg1;
821  ++valBeg1;
822  ++keyBeg2;
823  ++valBeg2;
824  }
825  }
826  // At most one of the two sequences will be nonempty.
827  std::copy (keyBeg1, keyEnd1, keyOut);
828  std::copy (valBeg1, valEnd1, valOut);
829  std::copy (keyBeg2, keyEnd2, keyOut);
830  std::copy (valBeg2, valEnd2, valOut);
831  }
832 
833  template<class KeyInputIterType>
834  size_t
835  keyMergeCount (KeyInputIterType keyBeg1, KeyInputIterType keyEnd1,
836  KeyInputIterType keyBeg2, KeyInputIterType keyEnd2)
837  {
838  size_t count = 0;
839  while (keyBeg1 != keyEnd1 && keyBeg2 != keyEnd2) {
840  if (*keyBeg1 < *keyBeg2) {
841  ++keyBeg1;
842  } else if (*keyBeg1 > *keyBeg2) {
843  ++keyBeg2;
844  } else { // *keyBeg1 == *keyBeg2
845  ++keyBeg1;
846  ++keyBeg2;
847  }
848  ++count;
849  }
850  // At most one of the two sequences will be nonempty.
851  count += static_cast<size_t> (keyEnd1 - keyBeg1) +
852  static_cast<size_t> (keyEnd2 - keyBeg2);
853  return count;
854  }
855 
856  namespace Details {
876  bool
877  congruent (const Teuchos::Comm<int>& comm1,
878  const Teuchos::Comm<int>& comm2);
879  } // namespace Details
880 
881 } // namespace Tpetra
882 
883 #endif // TPETRA_UTIL_HPP
Namespace Tpetra contains the class and methods constituting the Tpetra library.
void quicksort3(const IT1 &first1, const IT1 &last1, const IT2 &first2, const IT2 &last2, const IT3 &first3, const IT3 &last3)
Sort the first array using Quicksort, and apply the resulting permutation to the second and third arr...
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.
IT1 partition2(const IT1 &first1, const IT1 &last1, const IT2 &first2, const IT2 &last2, const IT1 &pivot)
Partition operation for quicksort2().
void quicksort2(const IT1 &first1, const IT1 &last1, const IT2 &first2, const IT2 &last2)
Sort the first array using Quicksort, and apply the resulting permutation to the second array...
Implementation details of Tpetra.
void sh_sort3(const IT1 &first1, const IT1 &last1, const IT2 &first2, const IT2 &last2, const IT3 &first3, const IT3 &last3)
Sort the first array using shell sort, and apply the resulting permutation to the second and third ar...
void merge2(IT1 &indResultOut, IT2 &valResultOut, IT1 indBeg, IT1 indEnd, IT2 valBeg, IT2 valEnd)
Merge values in place, additively, with the same index.
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...
IT getPivot(const IT &first, const IT &last)
Determines the pivot point as part of the quicksort routine.
MapType::iterator efficientAddOrUpdate(MapType &m, const KeyArgType &k, const ValueArgType &v)
Efficiently insert or replace an entry in an std::map.
void sort2(const IT1 &first1, const IT1 &last1, const IT2 &first2)
Sort the first array, and apply the resulting permutation to the second array.
Implementation details of sort routines used by Tpetra.
IT1 partition3(const IT1 &first1, const IT1 &last1, const IT2 &first2, const IT2 &last2, const IT3 &first3, const IT3 &last3, const IT1 &pivot)
Partition operation for quicksort3().
bool isAlreadySorted(const IT1 &first, const IT1 &last)
Determines whether or not a random access sequence is already sorted.
void sh_sort2(const IT1 &first1, const IT1 &last1, const IT2 &first2, const IT2 &last2)
Sort the first array using shell sort, and apply the resulting permutation to the second array...