Tpetra Matrix/Vector Services Version of the Day
Tpetra_Util.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //          Tpetra: Templated Linear Algebra Services Package
00005 //                 Copyright (2008) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 
00042 #ifndef TPETRA_UTIL_HPP
00043 #define TPETRA_UTIL_HPP
00044 
00071 #include "Tpetra_ConfigDefs.hpp" // for map, vector, string, and iostream
00072 #include <iterator>
00073 #include <algorithm>
00074 #include <Teuchos_Utils.hpp>
00075 #include <Teuchos_Assert.hpp>
00076 #include <sstream>
00077 
00078 #if defined(HAVE_TPETRA_THROW_EFFICIENCY_WARNINGS) || defined(HAVE_TPETRA_PRINT_EFFICIENCY_WARNINGS)
00079 
00080 
00081 
00082 
00083 
00084 
00085 
00086 
00087 
00088 
00089 
00090 
00091 
00092 
00093 
00094 
00095 
00096 
00097 
00098 
00099 
00100 
00101 
00102 
00103 
00104 
00105 
00106 
00107 
00108 
00109 
00110 
00111 
00112 #define TPETRA_EFFICIENCY_WARNING(throw_exception_test,Exception,msg)  \
00113 { \
00114   const bool tpetraEfficiencyWarningTest = (throw_exception_test); \
00115   if (tpetraEfficiencyWarningTest) { \
00116     std::ostringstream errStream; \
00117     errStream << Teuchos::typeName(*this) << ":" << std::endl; \
00118     errStream << "Efficiency warning: " << #throw_exception_test << std::endl; \
00119     errStream << msg; \
00120     std::string err = errStream.str(); \
00121     if (TPETRA_PRINTS_EFFICIENCY_WARNINGS && tpetraEfficiencyWarningTest) { \
00122       std::cerr << err << std::endl; \
00123     } \
00124     TEUCHOS_TEST_FOR_EXCEPTION(TPETRA_THROWS_EFFICIENCY_WARNINGS && tpetraEfficiencyWarningTest, Exception, err); \
00125   } \
00126 }
00127 #else
00128 
00129 
00130 
00131 
00132 
00133 
00134 
00135 
00136 
00137 
00138 
00139 
00140 
00141 
00142 
00143 
00144 
00145 
00146 
00147 
00148 
00149 
00150 
00151 
00152 
00153 
00154 
00155 
00156 
00157 
00158 
00159 
00160 
00161 #define TPETRA_EFFICIENCY_WARNING(throw_exception_test,Exception,msg)
00162 #endif
00163 
00164 // handle an abuse warning, according to HAVE_TPETRA_THROW_ABUSE_WARNINGS and HAVE_TPETRA_PRINT_ABUSE_WARNINGS
00165 #if defined(HAVE_TPETRA_THROW_ABUSE_WARNINGS) || defined(HAVE_TPETRA_PRINT_ABUSE_WARNINGS)
00166 
00167 #define TPETRA_ABUSE_WARNING(throw_exception_test,Exception,msg)                               \
00168 {                                                                                              \
00169   std::ostringstream errStream;                                                                \
00170   errStream << Teuchos::typeName(*this) << msg;                                                \
00171   std::string err = errStream.str();                                                           \
00172   if (TPETRA_PRINTS_ABUSE_WARNINGS && (throw_exception_test)) {                                \
00173     std::cerr << err << std::endl;                                                             \
00174   }                                                                                            \
00175   TEUCHOS_TEST_FOR_EXCEPTION(TPETRA_THROWS_ABUSE_WARNINGS && (throw_exception_test), Exception, err);  \
00176 }
00177 #else
00178 
00179 #define TPETRA_ABUSE_WARNING(throw_exception_test,Exception,msg)
00180 #endif
00181 
00211 #define SHARED_TEST_FOR_EXCEPTION(throw_exception_test,Exception,msg,comm) \
00212 { \
00213     using Teuchos::outArg; \
00214     const int lcl_throw_exception = (throw_exception_test) ? Teuchos::rank(comm)+1 : 0; \
00215     int gbl_throw; \
00216     Teuchos::reduceAll(comm,Teuchos::REDUCE_MAX,lcl_throw_exception,outArg(gbl_throw)); \
00217     TEUCHOS_TEST_FOR_EXCEPTION(gbl_throw,Exception,  \
00218       msg << " Failure on at least one process, including process " << gbl_throw-1 << "." << std::endl); \
00219 }
00220 
00221 #ifdef HAVE_TEUCHOS_DEBUG
00222 
00223 #define SWITCHED_TEST_FOR_EXCEPTION(throw_exception_test,Exception,msg,comm) \
00224 { \
00225     SHARED_TEST_FOR_EXCEPTION(throw_exception_test,Exception,msg,comm); \
00226 }
00227 #else
00228 
00229 #define SWITCHED_TEST_FOR_EXCEPTION(throw_exception_test,Exception,msg,comm) \
00230 { \
00231     TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test,Exception,msg); \
00232 }
00233 #endif
00234 
00235 namespace Tpetra {
00236 
00248   template<typename MapType, typename KeyArgType, typename ValueArgType>
00249   typename MapType::iterator efficientAddOrUpdate(MapType& m,
00250                           const KeyArgType & k,
00251                           const ValueArgType & v)
00252   {
00253     typename MapType::iterator lb = m.lower_bound(k);
00254     if(lb != m.end() && !(m.key_comp()(k, lb->first))) {
00255       lb->second = v;
00256       return(lb);
00257     }
00258     else {
00259       typedef typename MapType::value_type MVT;
00260       return(m.insert(lb, MVT(k, v)));
00261     }
00262   }
00263 
00271   namespace SortDetails {
00272 
00281   template<class IT1>
00282   bool isAlreadySorted(const IT1& first, const IT1& last){
00283     typedef typename std::iterator_traits<IT1>::difference_type DT;
00284     DT myit =OrdinalTraits<DT>::one();
00285     const DT sz  = last - first;
00286     for(;myit < sz; ++myit){
00287       if(first[myit] < first[myit-1]){
00288         return false;
00289       }
00290     }
00291     return true;
00292   }
00293 
00303   template<class IT>
00304   IT getPivot(const IT& first, const IT& last){
00305     IT pivot(first+(last-first)/2);
00306     if(*first<=*pivot && *(last-1)<=*first) pivot=first;
00307     else if(*(last-1)<=*pivot && *first<= *(last-1)) pivot = last-1;
00308     return pivot;
00309   }
00310 
00325   template<class IT1, class IT2>
00326   IT1 partition2(
00327     const IT1& first1,
00328     const IT1& last1,
00329     const IT2& first2,
00330     const IT2& last2,
00331     const IT1& pivot)
00332   {
00333     typename std::iterator_traits<IT1>::value_type piv(*pivot);
00334     std::swap(*pivot, *(last1-1));
00335     std::swap(first2[(pivot-first1)], *(last2-1));
00336     IT1 store1=first1;
00337     for(IT1 it=first1; it!=last1-1; ++it){
00338       if(*it<=piv){
00339         std::swap(*store1, *it);
00340         std::swap(first2[(store1-first1)], first2[(it-first1)]);
00341         ++store1;
00342       }
00343     }
00344     std::swap(*(last1-1), *store1);
00345     std::swap(*(last2-1), first2[store1-first1]);
00346     return store1;
00347   }
00348 
00365   template<class IT1, class IT2, class IT3>
00366   IT1 partition3(
00367     const IT1& first1,
00368     const IT1& last1,
00369     const IT2& first2,
00370     const IT2& last2,
00371     const IT3& first3,
00372     const IT3& last3,
00373     const IT1& pivot)
00374   {
00375     typename std::iterator_traits<IT1>::value_type piv(*pivot);
00376     std::swap(*pivot, *(last1-1));
00377     std::swap(first2[(pivot-first1)], *(last2-1));
00378     std::swap(first3[(pivot-first1)], *(last3-1));
00379     IT1 store1=first1;
00380     for(IT1 it=first1; it!=last1-1; ++it){
00381       if(*it<=piv){
00382         std::swap(*store1, *it);
00383         std::swap(first2[(store1-first1)], first2[(it-first1)]);
00384         std::swap(first3[(store1-first1)], first3[(it-first1)]);
00385         ++store1;
00386       }
00387     }
00388     std::swap(*(last1-1), *store1);
00389     std::swap(*(last2-1), first2[store1-first1]);
00390     std::swap(*(last3-1), first3[store1-first1]);
00391     return store1;
00392   }
00393 
00404   template<class IT1, class IT2>
00405   void quicksort2(
00406     const IT1& first1,
00407     const IT1& last1,
00408     const IT2& first2,
00409     const IT2& last2)
00410   {
00411     typedef typename std::iterator_traits<IT1>::difference_type DT;
00412     DT DT1 = OrdinalTraits<DT>::one();
00413     if(last1-first1 > DT1){
00414       IT1 pivot = getPivot(first1, last1);
00415       pivot = partition2(first1, last1, first2, last2, pivot);
00416       quicksort2(first1, pivot, first2, first2+(pivot-first1));
00417       quicksort2(pivot+1, last1, first2+(pivot-first1)+1, last2);
00418     }
00419   }
00420 
00433   template<class IT1, class IT2, class IT3>
00434   void quicksort3(
00435     const IT1& first1,
00436     const IT1& last1,
00437     const IT2& first2,
00438     const IT2& last2,
00439     const IT3& first3,
00440     const IT3& last3)
00441   {
00442     typedef typename std::iterator_traits<IT1>::difference_type DT;
00443     DT DT1 = OrdinalTraits<DT>::one();
00444     if(last1-first1 > DT1){
00445       IT1 pivot = getPivot(first1, last1);
00446       pivot = partition3(first1, last1, first2, last2, first3, last3, pivot);
00447       quicksort3(first1, pivot, first2, first2+(pivot-first1), first3, first3+(pivot-first1));
00448       quicksort3(pivot+1, last1, first2+(pivot-first1)+1, last2, first3+(pivot-first1)+1, last3);
00449     }
00450   }
00451 
00458   template<class IT1, class IT2, class IT3>
00459   void sh_sort3(
00460     const IT1& first1,
00461     const IT1& last1,
00462     const IT2& first2,
00463     const IT2& last2,
00464     const IT3& first3,
00465     const IT3& last3)
00466    {
00467         typedef typename std::iterator_traits<IT1>::difference_type DT;
00468         DT n = last1 - first1;
00469         DT m = n / 2;
00470         DT z = OrdinalTraits<DT>::zero();
00471         while (m > z)
00472         {
00473             DT max = n - m;
00474             for (DT j = 0; j < max; j++)
00475             {
00476                 for (DT k = j; k >= 0; k-=m)
00477                 {
00478                     if (first1[k+m] >= first1[k])
00479                         break;
00480                     std::swap(first1[k+m], first1[k]);
00481                     std::swap(first2[k+m], first2[k]);
00482                     std::swap(first3[k+m], first3[k]);
00483                 }
00484             }
00485             m = m/2;
00486         }
00487    }
00488 
00495   template<class IT1, class IT2>
00496   void sh_sort2(
00497     const IT1& first1,
00498     const IT1& last1,
00499     const IT2& first2,
00500     const IT2& last2)
00501    {
00502         typedef typename std::iterator_traits<IT1>::difference_type DT;
00503         DT n = last1 - first1;
00504         DT m = n / 2;
00505         DT z = OrdinalTraits<DT>::zero();
00506         while (m > z)
00507         {
00508             DT max = n - m;
00509             for (DT j = 0; j < max; j++)
00510             {
00511                 for (DT k = j; k >= 0; k-=m)
00512                 {
00513                     if (first1[k+m] >= first1[k])
00514                         break;
00515                     std::swap(first1[k+m], first1[k]);
00516                     std::swap(first2[k+m], first2[k]);
00517                 }
00518             }
00519             m = m/2;
00520         }
00521    }
00522 
00523   } //end namespace SortDetails
00524 
00525 
00544   template<class IT1, class IT2>
00545   void sort2(const IT1 &first1, const IT1 &last1, const IT2 &first2) {
00546     // Quicksort uses best-case N log N time whether or not the input
00547     // data is sorted.  However, the common case in Tpetra is that the
00548     // input data are sorted, so we first check whether this is the
00549     // case before reverting to Quicksort.
00550     if(SortDetails::isAlreadySorted(first1, last1)){
00551       return;
00552     }
00553     SortDetails::sh_sort2(first1, last1, first2, first2+(last1-first1));
00554 #ifdef HAVE_TPETRA_DEBUG
00555     if(!SortDetails::isAlreadySorted(first1, last1)){
00556       std::cout << "Trouble: sort() did not sort !!" << std::endl;
00557       return;
00558     }
00559 #endif
00560   }
00561 
00562 
00578   template<class IT1, class IT2, class IT3>
00579   void sort3(const IT1 &first1, const IT1 &last1, const IT2 &first2,
00580     const IT3 &first3)
00581   {
00582     // Quicksort uses best-case N log N time whether or not the input
00583     // data is sorted.  However, the common case in Tpetra is that the
00584     // input data are sorted, so we first check whether this is the
00585     // case before reverting to Quicksort.
00586     if(SortDetails::isAlreadySorted(first1, last1)){
00587       return;
00588     }
00589     SortDetails::sh_sort3(first1, last1, first2, first2+(last1-first1), first3,
00590                     first3+(last1-first1));
00591 
00592 #ifdef HAVE_TPETRA_DEBUG
00593     if(!SortDetails::isAlreadySorted(first1, last1)){
00594         std::cout << " Trouble sort did not actually sort... !!!!!!" <<
00595                         std::endl;
00596       return;
00597     }
00598 #endif
00599   }
00600 
00644   template<class IT1, class IT2>
00645   void
00646   merge2 (IT1& indResultOut, IT2& valResultOut,
00647           IT1 indBeg, IT1 indEnd,
00648           IT2 valBeg, IT2 valEnd)
00649   {
00650     if (indBeg == indEnd) {
00651       indResultOut = indBeg; // It's allowed for indResultOut to alias indEnd
00652       valResultOut = valBeg; // It's allowed for valResultOut to alias valEnd
00653     }
00654     else {
00655       IT1 indResult = indBeg;
00656       IT2 valResult = valBeg;
00657       if (indBeg != indEnd) {
00658         ++indBeg;
00659         ++valBeg;
00660         while (indBeg != indEnd) {
00661           if (*indResult == *indBeg) { // adjacent column indices equal
00662             *valResult += *valBeg; // merge entries by adding their values together
00663           } else { // adjacent column indices not equal
00664             *(++indResult) = *indBeg; // shift over the index
00665             *(++valResult) = *valBeg; // shift over the value
00666           }
00667           ++indBeg;
00668           ++valBeg;
00669         }
00670         ++indResult; // exclusive end of merged result
00671         ++valResult; // exclusive end of merged result
00672 
00673         // mfh 24 Feb 2014: Setting these is technically correct, but
00674         // since the resulting variables aren't used after this point,
00675         // it may result in a build error.
00676         //
00677         // indEnd = indResult;
00678         // valEnd = valResult;
00679       }
00680       indResultOut = indResult;
00681       valResultOut = valResult;
00682     }
00683   }
00684 
00733   template<class IT1, class IT2, class BinaryFunction>
00734   void
00735   merge2 (IT1& indResultOut, IT2& valResultOut,
00736           IT1 indBeg, IT1 indEnd,
00737           IT2 valBeg, IT2 valEnd,
00738           BinaryFunction f)
00739   {
00740     if (indBeg == indEnd) {
00741       indResultOut = indBeg; // It's allowed for indResultOut to alias indEnd
00742       valResultOut = valBeg; // It's allowed for valResultOut to alias valEnd
00743     }
00744     else {
00745       IT1 indResult = indBeg;
00746       IT2 valResult = valBeg;
00747       if (indBeg != indEnd) {
00748         ++indBeg;
00749         ++valBeg;
00750         while (indBeg != indEnd) {
00751           if (*indResult == *indBeg) { // adjacent column indices equal
00752             *valResult = f (*valResult, *valBeg); // merge entries via values
00753           } else { // adjacent column indices not equal
00754             *(++indResult) = *indBeg; // shift over the index
00755             *(++valResult) = *valBeg; // shift over the value
00756           }
00757           ++indBeg;
00758           ++valBeg;
00759         }
00760         ++indResult; // exclusive end of merged result
00761         ++valResult; // exclusive end of merged result
00762         indEnd = indResult;
00763         valEnd = valResult;
00764       }
00765       indResultOut = indResult;
00766       valResultOut = valResult;
00767     }
00768   }
00769 
00770 
00799   template<class KeyInputIterType, class ValueInputIterType,
00800            class KeyOutputIterType, class ValueOutputIterType,
00801            class BinaryFunction>
00802   void
00803   keyValueMerge (KeyInputIterType keyBeg1, KeyInputIterType keyEnd1,
00804                  ValueInputIterType valBeg1, ValueInputIterType valEnd1,
00805                  KeyInputIterType keyBeg2, KeyInputIterType keyEnd2,
00806                  ValueInputIterType valBeg2, ValueInputIterType valEnd2,
00807                  KeyOutputIterType keyOut, ValueOutputIterType valOut,
00808                  BinaryFunction f)
00809   {
00810     while (keyBeg1 != keyEnd1 && keyBeg2 != keyEnd2) {
00811       if (*keyBeg1 < *keyBeg2) {
00812         *keyOut++ = *keyBeg1++;
00813         *valOut++ = *valBeg1++;
00814       } else if (*keyBeg1 > *keyBeg2) {
00815         *keyOut++ = *keyBeg2++;
00816         *valOut++ = *valBeg2++;
00817       } else { // *keyBeg1 == *keyBeg2
00818         *keyOut++ = *keyBeg1;
00819         *valOut++ = f (*valBeg1, *valBeg2);
00820         ++keyBeg1;
00821         ++valBeg1;
00822         ++keyBeg2;
00823         ++valBeg2;
00824       }
00825     }
00826     // At most one of the two sequences will be nonempty.
00827     std::copy (keyBeg1, keyEnd1, keyOut);
00828     std::copy (valBeg1, valEnd1, valOut);
00829     std::copy (keyBeg2, keyEnd2, keyOut);
00830     std::copy (valBeg2, valEnd2, valOut);
00831   }
00832 
00833   template<class KeyInputIterType>
00834   size_t
00835   keyMergeCount (KeyInputIterType keyBeg1, KeyInputIterType keyEnd1,
00836                  KeyInputIterType keyBeg2, KeyInputIterType keyEnd2)
00837   {
00838     size_t count = 0;
00839     while (keyBeg1 != keyEnd1 && keyBeg2 != keyEnd2) {
00840       if (*keyBeg1 < *keyBeg2) {
00841         ++keyBeg1;
00842       } else if (*keyBeg1 > *keyBeg2) {
00843         ++keyBeg2;
00844       } else { // *keyBeg1 == *keyBeg2
00845         ++keyBeg1;
00846         ++keyBeg2;
00847       }
00848       ++count;
00849     }
00850     // At most one of the two sequences will be nonempty.
00851     count += static_cast<size_t> (keyEnd1 - keyBeg1) +
00852       static_cast<size_t> (keyEnd2 - keyBeg2);
00853     return count;
00854   }
00855 
00856   namespace Details {
00876     bool
00877     congruent (const Teuchos::Comm<int>& comm1,
00878                const Teuchos::Comm<int>& comm2);
00879   } // namespace Details
00880 
00881 } // namespace Tpetra
00882 
00883 #endif // TPETRA_UTIL_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines