00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
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
00126 enum SType {
00127 LM, SM,
00128 LR, SR,
00129 LI, SI
00130 };
00131 SType which_;
00132
00133
00134 template <class LTorGT>
00135 struct compMag {
00136
00137 bool operator()(MagnitudeType, MagnitudeType);
00138
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
00146 bool operator()(std::pair<MagnitudeType,MagnitudeType>, std::pair<MagnitudeType,MagnitudeType>);
00147
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
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
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
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
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
00268
00269
00270
00271
00272
00273
00274
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
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
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
00341
00342 if (perm == Teuchos::null) {
00343
00344
00345
00346 std::vector< std::pair<MagnitudeType,MagnitudeType> > pairs(n);
00347
00348
00349
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
00387
00388
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
00407
00408 std::vector< std::pair< std::pair<MagnitudeType,MagnitudeType>, int > > pairs(n);
00409
00410
00411
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
00447
00448
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 }
00538
00539 #endif // ANASAZI_BASIC_SORT_HPP
00540