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., 59 Temple Place, Suite 330, Boston, MA 02111-1307
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       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     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     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       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       switch (which_) {
00247         case LM:
00248           std::sort(evals.begin(),evals.begin()+n,compMag<greater_mt>());
00249           break;
00250         case SM:
00251           std::sort(evals.begin(),evals.begin()+n,compMag<less_mt>());
00252           break;
00253         case LR:
00254           std::sort(evals.begin(),evals.begin()+n,compAlg<greater_mt>());
00255           break;
00256         case SR:
00257           std::sort(evals.begin(),evals.begin()+n,compAlg<less_mt>());
00258           break;
00259         case SI:
00260         case LI:
00261           TEST_FOR_EXCEPTION(true, SortManagerError, "Anasazi::BasicSort::sort(r): LI or SI sorting invalid for real scalar types." );
00262           break;
00263       }
00264     }
00265     else {
00266       // 
00267       // if permutation index is required, we must sort the two at once
00268       // in this case, we arrange a pair structure: <value,index>
00269       // default comparison operator for pair<t1,t2> is lexographic:
00270       //    compare first t1, then t2
00271       // this works fine for us here.
00272       //
00273 
00274       // copy the values and indices into the pair structure
00275       std::vector< std::pair<MagnitudeType,int> > pairs(n);
00276       for (int i=0; i<n; i++) {
00277         pairs[i] = std::make_pair(evals[i],i);
00278       }
00279 
00280       // sort the pair structure
00281       switch (which_) {
00282         case LM:
00283           std::sort(pairs.begin(),pairs.begin()+n,compMag<greater_mt>());
00284           break;
00285         case SM:
00286           std::sort(pairs.begin(),pairs.begin()+n,compMag<less_mt>());
00287           break;
00288         case LR:
00289           std::sort(pairs.begin(),pairs.begin()+n,compAlg<greater_mt>());
00290           break;
00291         case SR:
00292           std::sort(pairs.begin(),pairs.begin()+n,compAlg<less_mt>());
00293           break;
00294         case SI:
00295         case LI:
00296           TEST_FOR_EXCEPTION(true, SortManagerError, "Anasazi::BasicSort::sort(r): LI or SI sorting invalid for real scalar types." );
00297           break;
00298       }
00299 
00300       // copy the values and indices out of the pair structure
00301       std::transform(pairs.begin(),pairs.end(),evals.begin(),sel1st< std::pair<MagnitudeType,int> >());
00302       std::transform(pairs.begin(),pairs.end(),perm->begin(),sel2nd< std::pair<MagnitudeType,int> >());
00303     }
00304   }
00305 
00306 
00307   template<class T1, class T2>
00308   class MakePairOp {
00309   public:
00310     std::pair<T1,T2> operator()( const T1 &t1, const T2 &t2 )
00311       { return std::make_pair(t1, t2); }
00312   };
00313 
00314 
00315   template<class MagnitudeType>
00316   void BasicSort<MagnitudeType>::sort(std::vector<MagnitudeType> &r_evals, 
00317                                       std::vector<MagnitudeType> &i_evals, 
00318                                       Teuchos::RCP< std::vector<int> > perm,
00319                                       int n) const 
00320   {
00321 
00322     typedef typename std::vector<MagnitudeType>::iterator r_eval_iter_t;
00323     typedef typename std::vector<MagnitudeType>::iterator i_eval_iter_t;
00324 
00325     TEST_FOR_EXCEPTION(n < -1, std::invalid_argument, "Anasazi::BasicSort::sort(r,i): n must be n >= 0 or n == -1.");
00326     if (n == -1) {
00327       n = r_evals.size() < i_evals.size() ? r_evals.size() : i_evals.size();
00328     }
00329     TEST_FOR_EXCEPTION(r_evals.size() < (unsigned int) n || i_evals.size() < (unsigned int) n,
00330                        std::invalid_argument, "Anasazi::BasicSort::sort(r,i): eigenvalue vector size isn't consistent with n.");
00331     if (perm != Teuchos::null) {
00332       TEST_FOR_EXCEPTION(perm->size() < (unsigned int) n,
00333                          std::invalid_argument, "Anasazi::BasicSort::sort(r,i): permutation vector size isn't consistent with n.");
00334     }
00335 
00336     typedef std::greater<MagnitudeType> greater_mt;
00337     typedef std::less<MagnitudeType>    less_mt;
00338 
00339     //
00340     // put values into pairs
00341     //
00342     if (perm == Teuchos::null) {
00343       //
00344       // not permuting, so we don't need indices in the pairs
00345       // 
00346       std::vector< std::pair<MagnitudeType,MagnitudeType> > pairs(n);
00347       // for LM,SM, the order doesn't matter
00348       // for LI,SI, the imaginary goes first
00349       // for LR,SR, the real goes in first
00350       switch (which_) {
00351         case LR:
00352         case SR:
00353         case LM:
00354         case SM:
00355           std::transform(
00356             r_evals.begin(), r_evals.begin()+n,
00357             i_evals.begin(), pairs.begin(),
00358             MakePairOp<MagnitudeType,MagnitudeType>());
00359           break;
00360         case LI:
00361         case SI:
00362           std::transform(
00363             i_evals.begin(), i_evals.begin()+n,
00364             r_evals.begin(), pairs.begin(),
00365             MakePairOp<MagnitudeType,MagnitudeType>());
00366           break;
00367       }
00368 
00369       switch (which_) {
00370         case LR:
00371         case LI:
00372           std::sort(pairs.begin(),pairs.end(),compAlg<greater_mt>());
00373           break;
00374         case SR:
00375         case SI:
00376           std::sort(pairs.begin(),pairs.end(),compAlg<less_mt>());
00377           break;
00378         case LM:
00379           std::sort(pairs.begin(),pairs.end(),compMag2<greater_mt>());
00380           break;
00381         case SM:
00382           std::sort(pairs.begin(),pairs.end(),compMag2<less_mt>());
00383           break;
00384       }
00385 
00386       // extract the values
00387       // for LM,SM,LR,SR: order is (real,imag)
00388       // for LI,SI: order is (imag,real)
00389       switch (which_) {
00390         case LR:
00391         case SR:
00392         case LM:
00393         case SM:
00394           std::transform(pairs.begin(),pairs.end(),r_evals.begin(),sel1st< std::pair<MagnitudeType,MagnitudeType> >());
00395           std::transform(pairs.begin(),pairs.end(),i_evals.begin(),sel2nd< std::pair<MagnitudeType,MagnitudeType> >());
00396           break;
00397         case LI:
00398         case SI:
00399           std::transform(pairs.begin(),pairs.end(),r_evals.begin(),sel2nd< std::pair<MagnitudeType,MagnitudeType> >());
00400           std::transform(pairs.begin(),pairs.end(),i_evals.begin(),sel1st< std::pair<MagnitudeType,MagnitudeType> >());
00401           break;
00402       }
00403     }
00404     else {
00405       //
00406       // permuting, we need indices in the pairs
00407       // 
00408       std::vector< std::pair< std::pair<MagnitudeType,MagnitudeType>, int > > pairs(n);
00409       // for LM,SM, the order doesn't matter
00410       // for LI,SI, the imaginary goes first
00411       // for LR,SR, the real goes in first
00412       switch (which_) {
00413         case LR:
00414         case SR:
00415         case LM:
00416         case SM:
00417           for (int i=0; i<n; i++) {
00418             pairs[i] = std::make_pair(std::make_pair(r_evals[i],i_evals[i]),i);
00419           }
00420           break;
00421         case LI:
00422         case SI:
00423           for (int i=0; i<n; i++) {
00424             pairs[i] = std::make_pair(std::make_pair(i_evals[i],r_evals[i]),i);
00425           }
00426           break;
00427       }
00428 
00429       switch (which_) {
00430         case LR:
00431         case LI:
00432           std::sort(pairs.begin(),pairs.end(),compAlg<greater_mt>());
00433           break;
00434         case SR:
00435         case SI:
00436           std::sort(pairs.begin(),pairs.end(),compAlg<less_mt>());
00437           break;
00438         case LM:
00439           std::sort(pairs.begin(),pairs.end(),compMag2<greater_mt>());
00440           break;
00441         case SM:
00442           std::sort(pairs.begin(),pairs.end(),compMag2<less_mt>());
00443           break;
00444       }
00445 
00446       // extract the values
00447       // for LM,SM,LR,SR: order is (real,imag)
00448       // for LI,SI: order is (imag,real)
00449       switch (which_) {
00450         case LR:
00451         case SR:
00452         case LM:
00453         case SM:
00454           for (int i=0; i<n; i++) {
00455             r_evals[i] = pairs[i].first.first;
00456             i_evals[i] = pairs[i].first.second;
00457             (*perm)[i] = pairs[i].second;
00458           }
00459           break;
00460         case LI:
00461         case SI:
00462           for (int i=0; i<n; i++) {
00463             i_evals[i] = pairs[i].first.first;
00464             r_evals[i] = pairs[i].first.second;
00465             (*perm)[i] = pairs[i].second;
00466           }
00467           break;
00468       }
00469     }
00470   }
00471 
00472 
00473   template<class MagnitudeType>
00474   template<class LTorGT>
00475   bool BasicSort<MagnitudeType>::compMag<LTorGT>::operator()(MagnitudeType v1, MagnitudeType v2)
00476   {
00477     typedef Teuchos::ScalarTraits<MagnitudeType> MTT;
00478     LTorGT comp;
00479     return comp( MTT::magnitude(v1), MTT::magnitude(v2) );
00480   }
00481 
00482   template<class MagnitudeType>
00483   template<class LTorGT>
00484   bool BasicSort<MagnitudeType>::compMag2<LTorGT>::operator()(std::pair<MagnitudeType,MagnitudeType> v1, std::pair<MagnitudeType,MagnitudeType> v2)
00485   {
00486     MagnitudeType m1 = v1.first*v1.first + v1.second*v1.second;
00487     MagnitudeType m2 = v2.first*v2.first + v2.second*v2.second;
00488     LTorGT comp;
00489     return comp( m1, m2 );
00490   }
00491 
00492   template<class MagnitudeType>
00493   template<class LTorGT>
00494   bool BasicSort<MagnitudeType>::compAlg<LTorGT>::operator()(MagnitudeType v1, MagnitudeType v2) 
00495   {
00496     LTorGT comp;
00497     return comp( v1, v2 );
00498   }
00499 
00500   template<class MagnitudeType>
00501   template<class LTorGT>
00502   template<class First, class Second>
00503   bool BasicSort<MagnitudeType>::compMag<LTorGT>::operator()(std::pair<First,Second> v1, std::pair<First,Second> v2) {
00504     return (*this)(v1.first,v2.first);
00505   }
00506 
00507   template<class MagnitudeType>
00508   template<class LTorGT>
00509   template<class First, class Second>
00510   bool BasicSort<MagnitudeType>::compMag2<LTorGT>::operator()(std::pair<First,Second> v1, std::pair<First,Second> v2) {
00511     return (*this)(v1.first,v2.first);
00512   }
00513 
00514   template<class MagnitudeType>
00515   template<class LTorGT>
00516   template<class First, class Second>
00517   bool BasicSort<MagnitudeType>::compAlg<LTorGT>::operator()(std::pair<First,Second> v1, std::pair<First,Second> v2) {
00518     return (*this)(v1.first,v2.first);
00519   }
00520 
00521   template <class MagnitudeType>
00522   template <typename pair_type>
00523   const typename pair_type::first_type &
00524   BasicSort<MagnitudeType>::sel1st<pair_type>::operator()(const pair_type &v) const 
00525   {
00526     return v.first;
00527   }
00528 
00529   template <class MagnitudeType>
00530   template <typename pair_type>
00531   const typename pair_type::second_type &
00532   BasicSort<MagnitudeType>::sel2nd<pair_type>::operator()(const pair_type &v) const 
00533   {
00534     return v.second;
00535   }
00536 
00537 } // namespace Anasazi
00538 
00539 #endif // ANASAZI_BASIC_SORT_HPP
00540 

Generated on Thu Dec 17 11:02:20 2009 for Anasazi by  doxygen 1.5.9