Anasazi Version of the Day
AnasaziBasicSort.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: Block Eigensolvers Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 
00033 #ifndef ANASAZI_BASIC_SORT_HPP
00034 #define ANASAZI_BASIC_SORT_HPP
00035 
00043 #include "AnasaziConfigDefs.hpp"
00044 #include "AnasaziSortManager.hpp"
00045 #include "Teuchos_LAPACK.hpp"
00046 #include "Teuchos_ScalarTraits.hpp"
00047 #include "Teuchos_ParameterList.hpp"
00048 
00049 namespace Anasazi {
00050 
00051   template<class MagnitudeType>
00052   class BasicSort : public SortManager<MagnitudeType> {
00053 
00054   public:
00055 
00061     BasicSort( Teuchos::ParameterList &pl );
00062 
00067     BasicSort( const std::string &which = "LM" );
00068 
00070     virtual ~BasicSort();
00071 
00073 
00082     void setSortType( const std::string &which );
00083 
00098     void sort(std::vector<MagnitudeType> &evals, Teuchos::RCP<std::vector<int> > perm = Teuchos::null, int n = -1) const;
00099 
00118     void sort(std::vector<MagnitudeType> &r_evals,
00119               std::vector<MagnitudeType> &i_evals,
00120               Teuchos::RCP<std::vector<int> > perm = Teuchos::null,
00121               int n = -1) const;
00122 
00123   protected:
00124 
00125     // enum for sort type
00126     enum SType {
00127       LM, SM,
00128       LR, SR,
00129       LI, SI
00130     };
00131     SType which_;
00132 
00133     // sorting methods
00134     template <class LTorGT>
00135     struct compMag {
00136       // for real-only LM,SM
00137       bool operator()(MagnitudeType, MagnitudeType);
00138       // for real-only LM,SM with permutation
00139       template <class First, class Second>
00140         bool operator()(std::pair<First,Second>, std::pair<First,Second>);
00141     };
00142 
00143     template <class LTorGT>
00144     struct compMag2 {
00145       // for real-imag LM,SM
00146       bool operator()(std::pair<MagnitudeType,MagnitudeType>, std::pair<MagnitudeType,MagnitudeType>);
00147       // for real-imag LM,SM with permutation
00148       template <class First, class Second>
00149         bool operator()(std::pair<First,Second>, std::pair<First,Second>);
00150     };
00151 
00152     template <class LTorGT>
00153     struct compAlg {
00154       // for real-imag LR,SR,LI,SI
00155       bool operator()(MagnitudeType, MagnitudeType);
00156       template <class First, class Second>
00157         bool operator()(std::pair<First,Second>, std::pair<First,Second>);
00158     };
00159 
00160     template <typename pair_type>
00161     struct sel1st
00162     {
00163       const typename pair_type::first_type &operator()(const pair_type &v) const;
00164     };
00165 
00166     template <typename pair_type>
00167     struct sel2nd
00168     {
00169       const typename pair_type::second_type &operator()(const pair_type &v) const;
00170     };
00171   };
00172 
00173 
00175   //  IMPLEMENTATION
00177 
00178   template<class MagnitudeType>
00179   BasicSort<MagnitudeType>::BasicSort(Teuchos::ParameterList &pl)
00180   {
00181     std::string which = "LM";
00182     which = pl.get("Sort Strategy",which);
00183     setSortType(which);
00184   }
00185 
00186   template<class MagnitudeType>
00187   BasicSort<MagnitudeType>::BasicSort(const std::string &which)
00188   {
00189     setSortType(which);
00190   }
00191 
00192   template<class MagnitudeType>
00193   BasicSort<MagnitudeType>::~BasicSort()
00194   {}
00195 
00196   template<class MagnitudeType>
00197   void BasicSort<MagnitudeType>::setSortType(const std::string &which)
00198   {
00199     // make upper case
00200     std::string whichlc(which);
00201     std::transform(which.begin(),which.end(),whichlc.begin(),(int(*)(int)) std::toupper);
00202     if (whichlc == "LM") {
00203       which_ = LM;
00204     }
00205     else if (whichlc == "SM") {
00206       which_ = SM;
00207     }
00208     else if (whichlc == "LR") {
00209       which_ = LR;
00210     }
00211     else if (whichlc == "SR") {
00212       which_ = SR;
00213     }
00214     else if (whichlc == "LI") {
00215       which_ = LI;
00216     }
00217     else if (whichlc == "SI") {
00218       which_ = SI;
00219     }
00220     else {
00221       TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Anasazi::BasicSort::setSortType(): sorting order is not valid");
00222     }
00223   }
00224 
00225   template<class MagnitudeType>
00226   void BasicSort<MagnitudeType>::sort(std::vector<MagnitudeType> &evals, Teuchos::RCP<std::vector<int> > perm, int n) const
00227   {
00228     TEUCHOS_TEST_FOR_EXCEPTION(n < -1, std::invalid_argument, "Anasazi::BasicSort::sort(r): n must be n >= 0 or n == -1.");
00229     if (n == -1) {
00230       n = evals.size();
00231     }
00232     TEUCHOS_TEST_FOR_EXCEPTION(evals.size() < (unsigned int) n,
00233                        std::invalid_argument, "Anasazi::BasicSort::sort(r): eigenvalue vector size isn't consistent with n.");
00234     if (perm != Teuchos::null) {
00235       TEUCHOS_TEST_FOR_EXCEPTION(perm->size() < (unsigned int) n,
00236                          std::invalid_argument, "Anasazi::BasicSort::sort(r): permutation vector size isn't consistent with n.");
00237     }
00238 
00239     typedef std::greater<MagnitudeType> greater_mt;
00240     typedef std::less<MagnitudeType>    less_mt;
00241 
00242     if (perm == Teuchos::null) {
00243       //
00244       // if permutation index is not required, just sort using the values
00245       //
00246       if (which_ == LM ) {
00247         std::sort(evals.begin(),evals.begin()+n,compMag<greater_mt>());
00248       }
00249       else if (which_ == SM) {
00250         std::sort(evals.begin(),evals.begin()+n,compMag<less_mt>());
00251       }
00252       else if (which_ == LR) {
00253         std::sort(evals.begin(),evals.begin()+n,compAlg<greater_mt>());
00254       }
00255       else if (which_ == SR) {
00256         std::sort(evals.begin(),evals.begin()+n,compAlg<less_mt>());
00257       }
00258       else {
00259         TEUCHOS_TEST_FOR_EXCEPTION(true, SortManagerError, "Anasazi::BasicSort::sort(r): LI or SI sorting invalid for real scalar types." );
00260       }
00261     }
00262     else {
00263       //
00264       // if permutation index is required, we must sort the two at once
00265       // in this case, we arrange a pair structure: <value,index>
00266       // default comparison operator for pair<t1,t2> is lexographic:
00267       //    compare first t1, then t2
00268       // this works fine for us here.
00269       //
00270 
00271       // copy the values and indices into the pair structure
00272       std::vector< std::pair<MagnitudeType,int> > pairs(n);
00273       for (int i=0; i<n; i++) {
00274         pairs[i] = std::make_pair(evals[i],i);
00275       }
00276 
00277       // sort the pair structure
00278       if (which_ == LM) {
00279         std::sort(pairs.begin(),pairs.begin()+n,compMag<greater_mt>());
00280       }
00281       else if (which_ == SM) {
00282         std::sort(pairs.begin(),pairs.begin()+n,compMag<less_mt>());
00283       }
00284       else if (which_ == LR) {
00285         std::sort(pairs.begin(),pairs.begin()+n,compAlg<greater_mt>());
00286       }
00287       else if (which_ == SR) {
00288         std::sort(pairs.begin(),pairs.begin()+n,compAlg<less_mt>());
00289       }
00290       else {
00291         TEUCHOS_TEST_FOR_EXCEPTION(true, SortManagerError, "Anasazi::BasicSort::sort(r): LI or SI sorting invalid for real scalar types." );
00292       }
00293 
00294       // copy the values and indices out of the pair structure
00295       std::transform(pairs.begin(),pairs.end(),evals.begin(),sel1st< std::pair<MagnitudeType,int> >());
00296       std::transform(pairs.begin(),pairs.end(),perm->begin(),sel2nd< std::pair<MagnitudeType,int> >());
00297     }
00298   }
00299 
00300 
00301   template<class T1, class T2>
00302   class MakePairOp {
00303   public:
00304     std::pair<T1,T2> operator()( const T1 &t1, const T2 &t2 )
00305       { return std::make_pair(t1, t2); }
00306   };
00307 
00308 
00309   template<class MagnitudeType>
00310   void BasicSort<MagnitudeType>::sort(std::vector<MagnitudeType> &r_evals,
00311                                       std::vector<MagnitudeType> &i_evals,
00312                                       Teuchos::RCP< std::vector<int> > perm,
00313                                       int n) const
00314   {
00315 
00316     //typedef typename std::vector<MagnitudeType>::iterator r_eval_iter_t; // unused
00317     //typedef typename std::vector<MagnitudeType>::iterator i_eval_iter_t; // unused
00318 
00319     TEUCHOS_TEST_FOR_EXCEPTION(n < -1, std::invalid_argument, "Anasazi::BasicSort::sort(r,i): n must be n >= 0 or n == -1.");
00320     if (n == -1) {
00321       n = r_evals.size() < i_evals.size() ? r_evals.size() : i_evals.size();
00322     }
00323     TEUCHOS_TEST_FOR_EXCEPTION(r_evals.size() < (unsigned int) n || i_evals.size() < (unsigned int) n,
00324                        std::invalid_argument, "Anasazi::BasicSort::sort(r,i): eigenvalue vector size isn't consistent with n.");
00325     if (perm != Teuchos::null) {
00326       TEUCHOS_TEST_FOR_EXCEPTION(perm->size() < (unsigned int) n,
00327                          std::invalid_argument, "Anasazi::BasicSort::sort(r,i): permutation vector size isn't consistent with n.");
00328     }
00329 
00330     typedef std::greater<MagnitudeType> greater_mt;
00331     typedef std::less<MagnitudeType>    less_mt;
00332 
00333     //
00334     // put values into pairs
00335     //
00336     if (perm == Teuchos::null) {
00337       //
00338       // not permuting, so we don't need indices in the pairs
00339       //
00340       std::vector< std::pair<MagnitudeType,MagnitudeType> > pairs(n);
00341       // for LM,SM, the order doesn't matter
00342       // for LI,SI, the imaginary goes first
00343       // for LR,SR, the real goes in first
00344       if (which_ == LR || which_ == SR || which_ == LM || which_ == SM) {
00345         std::transform(
00346           r_evals.begin(), r_evals.begin()+n,
00347           i_evals.begin(), pairs.begin(),
00348           MakePairOp<MagnitudeType,MagnitudeType>());
00349       }
00350       else {
00351         std::transform(
00352           i_evals.begin(), i_evals.begin()+n,
00353           r_evals.begin(), pairs.begin(),
00354           MakePairOp<MagnitudeType,MagnitudeType>());
00355       }
00356 
00357       if (which_ == LR || which_ == LI) {
00358         std::sort(pairs.begin(),pairs.end(),compAlg<greater_mt>());
00359       }
00360       else if (which_ == SR || which_ == SI) {
00361         std::sort(pairs.begin(),pairs.end(),compAlg<less_mt>());
00362       }
00363       else if (which_ == LM) {
00364         std::sort(pairs.begin(),pairs.end(),compMag2<greater_mt>());
00365       }
00366       else {
00367         std::sort(pairs.begin(),pairs.end(),compMag2<less_mt>());
00368       }
00369 
00370       // extract the values
00371       // for LM,SM,LR,SR: order is (real,imag)
00372       // for LI,SI: order is (imag,real)
00373       if (which_ == LR || which_ == SR || which_ == LM || which_ == SM) {
00374         std::transform(pairs.begin(),pairs.end(),r_evals.begin(),sel1st< std::pair<MagnitudeType,MagnitudeType> >());
00375         std::transform(pairs.begin(),pairs.end(),i_evals.begin(),sel2nd< std::pair<MagnitudeType,MagnitudeType> >());
00376       }
00377       else {
00378         std::transform(pairs.begin(),pairs.end(),r_evals.begin(),sel2nd< std::pair<MagnitudeType,MagnitudeType> >());
00379         std::transform(pairs.begin(),pairs.end(),i_evals.begin(),sel1st< std::pair<MagnitudeType,MagnitudeType> >());
00380       }
00381     }
00382     else {
00383       //
00384       // permuting, we need indices in the pairs
00385       //
00386       std::vector< std::pair< std::pair<MagnitudeType,MagnitudeType>, int > > pairs(n);
00387       // for LM,SM, the order doesn't matter
00388       // for LI,SI, the imaginary goes first
00389       // for LR,SR, the real goes in first
00390       if (which_ == LR || which_ == SR || which_ == LM || which_ == SM) {
00391         for (int i=0; i<n; i++) {
00392           pairs[i] = std::make_pair(std::make_pair(r_evals[i],i_evals[i]),i);
00393         }
00394       }
00395       else {
00396         for (int i=0; i<n; i++) {
00397           pairs[i] = std::make_pair(std::make_pair(i_evals[i],r_evals[i]),i);
00398         }
00399       }
00400 
00401       if (which_ == LR || which_ == LI) {
00402         std::sort(pairs.begin(),pairs.end(),compAlg<greater_mt>());
00403       }
00404       else if (which_ == SR || which_ == SI) {
00405         std::sort(pairs.begin(),pairs.end(),compAlg<less_mt>());
00406       }
00407       else if (which_ == LM) {
00408         std::sort(pairs.begin(),pairs.end(),compMag2<greater_mt>());
00409       }
00410       else {
00411         std::sort(pairs.begin(),pairs.end(),compMag2<less_mt>());
00412       }
00413 
00414       // extract the values
00415       // for LM,SM,LR,SR: order is (real,imag)
00416       // for LI,SI: order is (imag,real)
00417       if (which_ == LR || which_ == SR || which_ == LM || which_ == SM) {
00418         for (int i=0; i<n; i++) {
00419           r_evals[i] = pairs[i].first.first;
00420           i_evals[i] = pairs[i].first.second;
00421           (*perm)[i] = pairs[i].second;
00422         }
00423       }
00424       else {
00425         for (int i=0; i<n; i++) {
00426           i_evals[i] = pairs[i].first.first;
00427           r_evals[i] = pairs[i].first.second;
00428           (*perm)[i] = pairs[i].second;
00429         }
00430       }
00431     }
00432   }
00433 
00434 
00435   template<class MagnitudeType>
00436   template<class LTorGT>
00437   bool BasicSort<MagnitudeType>::compMag<LTorGT>::operator()(MagnitudeType v1, MagnitudeType v2)
00438   {
00439     typedef Teuchos::ScalarTraits<MagnitudeType> MTT;
00440     LTorGT comp;
00441     return comp( MTT::magnitude(v1), MTT::magnitude(v2) );
00442   }
00443 
00444   template<class MagnitudeType>
00445   template<class LTorGT>
00446   bool BasicSort<MagnitudeType>::compMag2<LTorGT>::operator()(std::pair<MagnitudeType,MagnitudeType> v1, std::pair<MagnitudeType,MagnitudeType> v2)
00447   {
00448     MagnitudeType m1 = v1.first*v1.first + v1.second*v1.second;
00449     MagnitudeType m2 = v2.first*v2.first + v2.second*v2.second;
00450     LTorGT comp;
00451     return comp( m1, m2 );
00452   }
00453 
00454   template<class MagnitudeType>
00455   template<class LTorGT>
00456   bool BasicSort<MagnitudeType>::compAlg<LTorGT>::operator()(MagnitudeType v1, MagnitudeType v2)
00457   {
00458     LTorGT comp;
00459     return comp( v1, v2 );
00460   }
00461 
00462   template<class MagnitudeType>
00463   template<class LTorGT>
00464   template<class First, class Second>
00465   bool BasicSort<MagnitudeType>::compMag<LTorGT>::operator()(std::pair<First,Second> v1, std::pair<First,Second> v2) {
00466     return (*this)(v1.first,v2.first);
00467   }
00468 
00469   template<class MagnitudeType>
00470   template<class LTorGT>
00471   template<class First, class Second>
00472   bool BasicSort<MagnitudeType>::compMag2<LTorGT>::operator()(std::pair<First,Second> v1, std::pair<First,Second> v2) {
00473     return (*this)(v1.first,v2.first);
00474   }
00475 
00476   template<class MagnitudeType>
00477   template<class LTorGT>
00478   template<class First, class Second>
00479   bool BasicSort<MagnitudeType>::compAlg<LTorGT>::operator()(std::pair<First,Second> v1, std::pair<First,Second> v2) {
00480     return (*this)(v1.first,v2.first);
00481   }
00482 
00483   template <class MagnitudeType>
00484   template <typename pair_type>
00485   const typename pair_type::first_type &
00486   BasicSort<MagnitudeType>::sel1st<pair_type>::operator()(const pair_type &v) const
00487   {
00488     return v.first;
00489   }
00490 
00491   template <class MagnitudeType>
00492   template <typename pair_type>
00493   const typename pair_type::second_type &
00494   BasicSort<MagnitudeType>::sel2nd<pair_type>::operator()(const pair_type &v) const
00495   {
00496     return v.second;
00497   }
00498 
00499 } // namespace Anasazi
00500 
00501 #endif // ANASAZI_BASIC_SORT_HPP
00502 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends