46#ifndef ANASAZI_BASIC_SORT_HPP
47#define ANASAZI_BASIC_SORT_HPP
58#include "Teuchos_LAPACK.hpp"
59#include "Teuchos_ScalarTraits.hpp"
60#include "Teuchos_ParameterList.hpp"
64 template<
class MagnitudeType>
80 BasicSort(
const std::string &which =
"LM" );
111 void sort(std::vector<MagnitudeType> &evals, Teuchos::RCP<std::vector<int> > perm = Teuchos::null,
int n = -1)
const;
131 void sort(std::vector<MagnitudeType> &r_evals,
132 std::vector<MagnitudeType> &i_evals,
133 Teuchos::RCP<std::vector<int> > perm = Teuchos::null,
147 template <
class LTorGT>
150 bool operator()(MagnitudeType, MagnitudeType);
152 template <
class First,
class Second>
153 bool operator()(std::pair<First,Second>, std::pair<First,Second>);
156 template <
class LTorGT>
159 bool operator()(std::pair<MagnitudeType,MagnitudeType>, std::pair<MagnitudeType,MagnitudeType>);
161 template <
class First,
class Second>
162 bool operator()(std::pair<First,Second>, std::pair<First,Second>);
165 template <
class LTorGT>
168 bool operator()(MagnitudeType, MagnitudeType);
169 template <
class First,
class Second>
170 bool operator()(std::pair<First,Second>, std::pair<First,Second>);
173 template <
typename pair_type>
176 const typename pair_type::first_type &operator()(
const pair_type &v)
const;
179 template <
typename pair_type>
182 const typename pair_type::second_type &operator()(
const pair_type &v)
const;
191 template<
class MagnitudeType>
194 std::string which =
"LM";
195 which = pl.get(
"Sort Strategy",which);
199 template<
class MagnitudeType>
205 template<
class MagnitudeType>
209 template<
class MagnitudeType>
213 std::string whichlc(which);
214 std::transform(which.begin(),which.end(),whichlc.begin(),(
int(*)(
int)) std::toupper);
215 if (whichlc ==
"LM") {
218 else if (whichlc ==
"SM") {
221 else if (whichlc ==
"LR") {
224 else if (whichlc ==
"SR") {
227 else if (whichlc ==
"LI") {
230 else if (whichlc ==
"SI") {
234 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Anasazi::BasicSort::setSortType(): sorting order is not valid");
238 template<
class MagnitudeType>
241 TEUCHOS_TEST_FOR_EXCEPTION(n < -1, std::invalid_argument,
"Anasazi::BasicSort::sort(r): n must be n >= 0 or n == -1.");
245 TEUCHOS_TEST_FOR_EXCEPTION(evals.size() < (
unsigned int) n,
246 std::invalid_argument,
"Anasazi::BasicSort::sort(r): eigenvalue vector size isn't consistent with n.");
247 if (perm != Teuchos::null) {
248 TEUCHOS_TEST_FOR_EXCEPTION(perm->size() < (
unsigned int) n,
249 std::invalid_argument,
"Anasazi::BasicSort::sort(r): permutation vector size isn't consistent with n.");
252 typedef std::greater<MagnitudeType> greater_mt;
253 typedef std::less<MagnitudeType> less_mt;
255 if (perm == Teuchos::null) {
260 std::sort(evals.begin(),evals.begin()+n,compMag<greater_mt>());
262 else if (which_ == SM) {
263 std::sort(evals.begin(),evals.begin()+n,compMag<less_mt>());
265 else if (which_ == LR) {
266 std::sort(evals.begin(),evals.begin()+n,compAlg<greater_mt>());
268 else if (which_ == SR) {
269 std::sort(evals.begin(),evals.begin()+n,compAlg<less_mt>());
272 TEUCHOS_TEST_FOR_EXCEPTION(
true,
SortManagerError,
"Anasazi::BasicSort::sort(r): LI or SI sorting invalid for real scalar types." );
285 std::vector< std::pair<MagnitudeType,int> > pairs(n);
286 for (
int i=0; i<n; i++) {
287 pairs[i] = std::make_pair(evals[i],i);
292 std::sort(pairs.begin(),pairs.begin()+n,compMag<greater_mt>());
294 else if (which_ == SM) {
295 std::sort(pairs.begin(),pairs.begin()+n,compMag<less_mt>());
297 else if (which_ == LR) {
298 std::sort(pairs.begin(),pairs.begin()+n,compAlg<greater_mt>());
300 else if (which_ == SR) {
301 std::sort(pairs.begin(),pairs.begin()+n,compAlg<less_mt>());
304 TEUCHOS_TEST_FOR_EXCEPTION(
true,
SortManagerError,
"Anasazi::BasicSort::sort(r): LI or SI sorting invalid for real scalar types." );
308 std::transform(pairs.begin(),pairs.end(),evals.begin(),sel1st< std::pair<MagnitudeType,int> >());
309 std::transform(pairs.begin(),pairs.end(),perm->begin(),sel2nd< std::pair<MagnitudeType,int> >());
314 template<
class T1,
class T2>
317 std::pair<T1,T2> operator()(
const T1 &t1,
const T2 &t2 )
318 {
return std::make_pair(t1, t2); }
322 template<
class MagnitudeType>
324 std::vector<MagnitudeType> &i_evals,
325 Teuchos::RCP< std::vector<int> > perm,
332 TEUCHOS_TEST_FOR_EXCEPTION(n < -1, std::invalid_argument,
"Anasazi::BasicSort::sort(r,i): n must be n >= 0 or n == -1.");
334 n = r_evals.size() < i_evals.size() ? r_evals.size() : i_evals.size();
336 TEUCHOS_TEST_FOR_EXCEPTION(r_evals.size() < (
unsigned int) n || i_evals.size() < (
unsigned int) n,
337 std::invalid_argument,
"Anasazi::BasicSort::sort(r,i): eigenvalue vector size isn't consistent with n.");
338 if (perm != Teuchos::null) {
339 TEUCHOS_TEST_FOR_EXCEPTION(perm->size() < (
unsigned int) n,
340 std::invalid_argument,
"Anasazi::BasicSort::sort(r,i): permutation vector size isn't consistent with n.");
343 typedef std::greater<MagnitudeType> greater_mt;
344 typedef std::less<MagnitudeType> less_mt;
349 if (perm == Teuchos::null) {
353 std::vector< std::pair<MagnitudeType,MagnitudeType> > pairs(n);
357 if (which_ == LR || which_ == SR || which_ == LM || which_ == SM) {
359 r_evals.begin(), r_evals.begin()+n,
360 i_evals.begin(), pairs.begin(),
361 MakePairOp<MagnitudeType,MagnitudeType>());
365 i_evals.begin(), i_evals.begin()+n,
366 r_evals.begin(), pairs.begin(),
367 MakePairOp<MagnitudeType,MagnitudeType>());
370 if (which_ == LR || which_ == LI) {
371 std::sort(pairs.begin(),pairs.end(),compAlg<greater_mt>());
373 else if (which_ == SR || which_ == SI) {
374 std::sort(pairs.begin(),pairs.end(),compAlg<less_mt>());
376 else if (which_ == LM) {
377 std::sort(pairs.begin(),pairs.end(),compMag2<greater_mt>());
380 std::sort(pairs.begin(),pairs.end(),compMag2<less_mt>());
386 if (which_ == LR || which_ == SR || which_ == LM || which_ == SM) {
387 std::transform(pairs.begin(),pairs.end(),r_evals.begin(),sel1st< std::pair<MagnitudeType,MagnitudeType> >());
388 std::transform(pairs.begin(),pairs.end(),i_evals.begin(),sel2nd< std::pair<MagnitudeType,MagnitudeType> >());
391 std::transform(pairs.begin(),pairs.end(),r_evals.begin(),sel2nd< std::pair<MagnitudeType,MagnitudeType> >());
392 std::transform(pairs.begin(),pairs.end(),i_evals.begin(),sel1st< std::pair<MagnitudeType,MagnitudeType> >());
399 std::vector< std::pair< std::pair<MagnitudeType,MagnitudeType>,
int > > pairs(n);
403 if (which_ == LR || which_ == SR || which_ == LM || which_ == SM) {
404 for (
int i=0; i<n; i++) {
405 pairs[i] = std::make_pair(std::make_pair(r_evals[i],i_evals[i]),i);
409 for (
int i=0; i<n; i++) {
410 pairs[i] = std::make_pair(std::make_pair(i_evals[i],r_evals[i]),i);
414 if (which_ == LR || which_ == LI) {
415 std::sort(pairs.begin(),pairs.end(),compAlg<greater_mt>());
417 else if (which_ == SR || which_ == SI) {
418 std::sort(pairs.begin(),pairs.end(),compAlg<less_mt>());
420 else if (which_ == LM) {
421 std::sort(pairs.begin(),pairs.end(),compMag2<greater_mt>());
424 std::sort(pairs.begin(),pairs.end(),compMag2<less_mt>());
430 if (which_ == LR || which_ == SR || which_ == LM || which_ == SM) {
431 for (
int i=0; i<n; i++) {
432 r_evals[i] = pairs[i].first.first;
433 i_evals[i] = pairs[i].first.second;
434 (*perm)[i] = pairs[i].second;
438 for (
int i=0; i<n; i++) {
439 i_evals[i] = pairs[i].first.first;
440 r_evals[i] = pairs[i].first.second;
441 (*perm)[i] = pairs[i].second;
448 template<
class MagnitudeType>
449 template<
class LTorGT>
452 typedef Teuchos::ScalarTraits<MagnitudeType> MTT;
454 return comp( MTT::magnitude(v1), MTT::magnitude(v2) );
457 template<
class MagnitudeType>
458 template<
class LTorGT>
459 bool BasicSort<MagnitudeType>::compMag2<LTorGT>::operator()(std::pair<MagnitudeType,MagnitudeType> v1, std::pair<MagnitudeType,MagnitudeType> v2)
461 MagnitudeType m1 = v1.first*v1.first + v1.second*v1.second;
462 MagnitudeType m2 = v2.first*v2.first + v2.second*v2.second;
464 return comp( m1, m2 );
467 template<
class MagnitudeType>
468 template<
class LTorGT>
469 bool BasicSort<MagnitudeType>::compAlg<LTorGT>::operator()(MagnitudeType v1, MagnitudeType v2)
472 return comp( v1, v2 );
475 template<
class MagnitudeType>
476 template<
class LTorGT>
477 template<
class First,
class Second>
478 bool BasicSort<MagnitudeType>::compMag<LTorGT>::operator()(std::pair<First,Second> v1, std::pair<First,Second> v2) {
479 return (*
this)(v1.first,v2.first);
482 template<
class MagnitudeType>
483 template<
class LTorGT>
484 template<
class First,
class Second>
485 bool BasicSort<MagnitudeType>::compMag2<LTorGT>::operator()(std::pair<First,Second> v1, std::pair<First,Second> v2) {
486 return (*
this)(v1.first,v2.first);
489 template<
class MagnitudeType>
490 template<
class LTorGT>
491 template<
class First,
class Second>
492 bool BasicSort<MagnitudeType>::compAlg<LTorGT>::operator()(std::pair<First,Second> v1, std::pair<First,Second> v2) {
493 return (*
this)(v1.first,v2.first);
496 template <
class MagnitudeType>
497 template <
typename pair_type>
498 const typename pair_type::first_type &
499 BasicSort<MagnitudeType>::sel1st<pair_type>::operator()(
const pair_type &v)
const
504 template <
class MagnitudeType>
505 template <
typename pair_type>
506 const typename pair_type::second_type &
507 BasicSort<MagnitudeType>::sel2nd<pair_type>::operator()(
const pair_type &v)
const
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Virtual base class which defines the interface between an eigensolver and a class whose job is the so...
An implementation of the Anasazi::SortManager that performs a collection of common sorting techniques...
virtual ~BasicSort()
Destructor.
BasicSort(Teuchos::ParameterList &pl)
Parameter list driven constructor.
void sort(std::vector< MagnitudeType > &evals, Teuchos::RCP< std::vector< int > > perm=Teuchos::null, int n=-1) const
Sort real eigenvalues, optionally returning the permutation vector.
void setSortType(const std::string &which)
Set sort type.
SortManagerError is thrown when the Anasazi::SortManager is unable to sort the numbers,...
Anasazi's templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.