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 
00045 #include "Tpetra_ConfigDefs.hpp" // for map, vector, string, and iostream 
00046 #include <iterator>
00047 #include <algorithm>
00048 #include <Teuchos_Utils.hpp>
00049 #include <Teuchos_Assert.hpp>
00050 #include <sstream>
00051 
00052 #if defined(HAVE_TPETRA_THROW_EFFICIENCY_WARNINGS) || defined(HAVE_TPETRA_PRINT_EFFICIENCY_WARNINGS)
00053 
00054 #define TPETRA_EFFICIENCY_WARNING(throw_exception_test,Exception,msg)                                 \
00055 {                                                                                                     \
00056   std::ostringstream errStream;                                                                       \
00057   errStream << Teuchos::typeName(*this) << msg;                                                       \
00058   std::string err = errStream.str();                                                                  \
00059   if (TPETRA_PRINTS_EFFICIENCY_WARNINGS && (throw_exception_test)) {                                  \
00060     std::cerr << err << std::endl;                                                                    \
00061   }                                                                                                   \
00062   TEUCHOS_TEST_FOR_EXCEPTION(TPETRA_THROWS_EFFICIENCY_WARNINGS && (throw_exception_test), Exception, err);    \
00063 }
00064 #else
00065 
00066 #define TPETRA_EFFICIENCY_WARNING(throw_exception_test,Exception,msg)
00067 #endif
00068 
00069 // handle an abuse warning, according to HAVE_TPETRA_THROW_ABUSE_WARNINGS and HAVE_TPETRA_PRINT_ABUSE_WARNINGS
00070 #if defined(HAVE_TPETRA_THROW_ABUSE_WARNINGS) || defined(HAVE_TPETRA_PRINT_ABUSE_WARNINGS)
00071 
00072 #define TPETRA_ABUSE_WARNING(throw_exception_test,Exception,msg)                               \
00073 {                                                                                              \
00074   std::ostringstream errStream;                                                                \
00075   errStream << Teuchos::typeName(*this) << msg;                                                \
00076   std::string err = errStream.str();                                                           \
00077   if (TPETRA_PRINTS_ABUSE_WARNINGS && (throw_exception_test)) {                                \
00078     std::cerr << err << std::endl;                                                             \
00079   }                                                                                            \
00080   TEUCHOS_TEST_FOR_EXCEPTION(TPETRA_THROWS_ABUSE_WARNINGS && (throw_exception_test), Exception, err);  \
00081 }
00082 #else
00083 
00084 #define TPETRA_ABUSE_WARNING(throw_exception_test,Exception,msg)
00085 #endif
00086 
00087 
00091 #define SHARED_TEST_FOR_EXCEPTION(throw_exception_test,Exception,msg,comm) \
00092 { \
00093     using Teuchos::outArg; \
00094     const int lcl_throw_exception = (throw_exception_test) ? Teuchos::rank(comm)+1 : 0; \
00095     int gbl_throw; \
00096     Teuchos::reduceAll(comm,Teuchos::REDUCE_MAX,lcl_throw_exception,outArg(gbl_throw)); \
00097     TEUCHOS_TEST_FOR_EXCEPTION(gbl_throw,Exception,  \
00098                        msg << " Failure on node " << gbl_throw-1 << "." << std::endl); \
00099 }
00100 
00101 #ifdef HAVE_TEUCHOS_DEBUG
00102 
00103 #define SWITCHED_TEST_FOR_EXCEPTION(throw_exception_test,Exception,msg,comm) \
00104 { \
00105     SHARED_TEST_FOR_EXCEPTION(throw_exception_test,Exception,msg,comm); \
00106 }
00107 #else 
00108 
00109 #define SWITCHED_TEST_FOR_EXCEPTION(throw_exception_test,Exception,msg,comm) \
00110 { \
00111     TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test,Exception,msg); \
00112 }
00113 #endif
00114 
00115 namespace Tpetra {
00116 
00144   template<typename MapType, typename KeyArgType, typename ValueArgType>
00145   typename MapType::iterator efficientAddOrUpdate(MapType& m, 
00146                           const KeyArgType & k, 
00147                           const ValueArgType & v) 
00148   {
00149     typename MapType::iterator lb = m.lower_bound(k);
00150     if(lb != m.end() && !(m.key_comp()(k, lb->first))) {
00151       lb->second = v;
00152       return(lb);
00153     }
00154     else {
00155       typedef typename MapType::value_type MVT;
00156       return(m.insert(lb, MVT(k, v)));
00157     }
00158   }
00159 
00167   namespace SortDetails {
00168 
00177   template<class IT1>
00178   bool isAlreadySorted(const IT1& first, const IT1& last){
00179     typedef typename std::iterator_traits<IT1>::difference_type DT;
00180     DT myit =OrdinalTraits<DT>::one();
00181     const DT sz  = last - first;
00182     for(;myit < sz; ++myit){
00183       if(first[myit] < first[myit-1]){
00184         return false;
00185       }
00186     }
00187     return true;
00188   }
00189 
00199   template<class IT>
00200   IT getPivot(const IT& first, const IT& last){
00201     IT pivot(first+(last-first)/2);
00202     if(*first<=*pivot && *(last-1)<=*first) pivot=first;
00203     else if(*(last-1)<=*pivot && *first<= *(last-1)) pivot = last-1;
00204     return pivot; 
00205   }
00206 
00221   template<class IT1, class IT2>
00222   IT1 partition2(
00223     const IT1& first1,
00224     const IT1& last1,
00225     const IT2& first2,
00226     const IT2& last2,
00227     const IT1& pivot)
00228   {
00229     typename std::iterator_traits<IT1>::value_type piv(*pivot);
00230     std::swap(*pivot, *(last1-1));
00231     std::swap(first2[(pivot-first1)], *(last2-1));
00232     IT1 store1=first1;
00233     for(IT1 it=first1; it!=last1-1; ++it){
00234       if(*it<=piv){
00235         std::swap(*store1, *it);
00236         std::swap(first2[(store1-first1)], first2[(it-first1)]);
00237         ++store1;
00238       }
00239     }
00240     std::swap(*(last1-1), *store1);
00241     std::swap(*(last2-1), first2[store1-first1]);
00242     return store1;
00243   }
00244 
00261   template<class IT1, class IT2, class IT3>
00262   IT1 partition3(
00263     const IT1& first1,
00264     const IT1& last1,
00265     const IT2& first2,
00266     const IT2& last2,
00267     const IT3& first3,
00268     const IT3& last3,
00269     const IT1& pivot)
00270   {
00271     typename std::iterator_traits<IT1>::value_type piv(*pivot);
00272     std::swap(*pivot, *(last1-1));
00273     std::swap(first2[(pivot-first1)], *(last2-1));
00274     std::swap(first3[(pivot-first1)], *(last3-1));
00275     IT1 store1=first1;
00276     for(IT1 it=first1; it!=last1-1; ++it){
00277       if(*it<=piv){
00278         std::swap(*store1, *it);
00279         std::swap(first2[(store1-first1)], first2[(it-first1)]);
00280         std::swap(first3[(store1-first1)], first3[(it-first1)]);
00281         ++store1;
00282       }
00283     }
00284     std::swap(*(last1-1), *store1);
00285     std::swap(*(last2-1), first2[store1-first1]);
00286     std::swap(*(last3-1), first3[store1-first1]);
00287     return store1;
00288   }
00289 
00300   template<class IT1, class IT2>
00301   void quicksort2(
00302     const IT1& first1,
00303     const IT1& last1,
00304     const IT2& first2,
00305     const IT2& last2)
00306   {
00307     typedef typename std::iterator_traits<IT1>::difference_type DT;
00308     DT DT1 = OrdinalTraits<DT>::one();
00309     if(last1-first1 > DT1){
00310       IT1 pivot = getPivot(first1, last1);
00311       pivot = partition2(first1, last1, first2, last2, pivot);
00312       quicksort2(first1, pivot, first2, first2+(pivot-first1));
00313       quicksort2(pivot+1, last1, first2+(pivot-first1)+1, last2);
00314     }
00315   }
00316 
00329   template<class IT1, class IT2, class IT3>
00330   void quicksort3(
00331     const IT1& first1,
00332     const IT1& last1,
00333     const IT2& first2,
00334     const IT2& last2,
00335     const IT3& first3,
00336     const IT3& last3)
00337   {
00338     typedef typename std::iterator_traits<IT1>::difference_type DT;
00339     DT DT1 = OrdinalTraits<DT>::one();
00340     if(last1-first1 > DT1){
00341       IT1 pivot = getPivot(first1, last1);
00342       pivot = partition3(first1, last1, first2, last2, first3, last3, pivot);
00343       quicksort3(first1, pivot, first2, first2+(pivot-first1), first3, first3+(pivot-first1));
00344       quicksort3(pivot+1, last1, first2+(pivot-first1)+1, last2, first3+(pivot-first1)+1, last3);
00345     }
00346   }
00347 
00348 
00349   } //end namespace SortDetails
00350 
00351 
00370   template<class IT1, class IT2>
00371   void sort2(const IT1 &first1, const IT1 &last1, const IT2 &first2) {
00372     // Quicksort uses best-case N log N time whether or not the input
00373     // data is sorted.  However, the common case in Tpetra is that the
00374     // input data are sorted, so we first check whether this is the
00375     // case before reverting to Quicksort.
00376     if(SortDetails::isAlreadySorted(first1, last1)){
00377       return;
00378     }
00379     SortDetails::quicksort2(first1, last1, first2, first2+(last1-first1));
00380   }
00381 
00382   
00397   template<class IT1, class IT2, class IT3>
00398   void sort3(const IT1 &first1, const IT1 &last1, const IT2 &first2, const IT3 &first3)
00399   {
00400     // Quicksort uses best-case N log N time whether or not the input
00401     // data is sorted.  However, the common case in Tpetra is that the
00402     // input data are sorted, so we first check whether this is the
00403     // case before reverting to Quicksort.
00404     if(SortDetails::isAlreadySorted(first1, last1)){
00405       return;
00406     }
00407     SortDetails::quicksort3(first1, last1, first2, first2+(last1-first1), first3, first3+(last1-first1));
00408   }
00409 
00410 } // namespace Tpetra
00411 
00412 
00413 #endif // TPETRA_UTIL_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines