Anasazi  Version of the Day
AnasaziBasicSort.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Anasazi: Block Eigensolvers Package
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
33 #ifndef ANASAZI_BASIC_SORT_HPP
34 #define ANASAZI_BASIC_SORT_HPP
35 
43 #include "AnasaziConfigDefs.hpp"
44 #include "AnasaziSortManager.hpp"
45 #include "Teuchos_LAPACK.hpp"
46 #include "Teuchos_ScalarTraits.hpp"
47 #include "Teuchos_ParameterList.hpp"
48 
49 namespace Anasazi {
50 
51  template<class MagnitudeType>
52  class BasicSort : public SortManager<MagnitudeType> {
53 
54  public:
55 
61  BasicSort( Teuchos::ParameterList &pl );
62 
67  BasicSort( const std::string &which = "LM" );
68 
70  virtual ~BasicSort();
71 
73 
82  void setSortType( const std::string &which );
83 
98  void sort(std::vector<MagnitudeType> &evals, Teuchos::RCP<std::vector<int> > perm = Teuchos::null, int n = -1) const;
99 
118  void sort(std::vector<MagnitudeType> &r_evals,
119  std::vector<MagnitudeType> &i_evals,
120  Teuchos::RCP<std::vector<int> > perm = Teuchos::null,
121  int n = -1) const;
122 
123  protected:
124 
125  // enum for sort type
126  enum SType {
127  LM, SM,
128  LR, SR,
129  LI, SI
130  };
131  SType which_;
132 
133  // sorting methods
134  template <class LTorGT>
135  struct compMag {
136  // for real-only LM,SM
137  bool operator()(MagnitudeType, MagnitudeType);
138  // for real-only LM,SM with permutation
139  template <class First, class Second>
140  bool operator()(std::pair<First,Second>, std::pair<First,Second>);
141  };
142 
143  template <class LTorGT>
144  struct compMag2 {
145  // for real-imag LM,SM
146  bool operator()(std::pair<MagnitudeType,MagnitudeType>, std::pair<MagnitudeType,MagnitudeType>);
147  // for real-imag LM,SM with permutation
148  template <class First, class Second>
149  bool operator()(std::pair<First,Second>, std::pair<First,Second>);
150  };
151 
152  template <class LTorGT>
153  struct compAlg {
154  // for real-imag LR,SR,LI,SI
155  bool operator()(MagnitudeType, MagnitudeType);
156  template <class First, class Second>
157  bool operator()(std::pair<First,Second>, std::pair<First,Second>);
158  };
159 
160  template <typename pair_type>
161  struct sel1st
162  {
163  const typename pair_type::first_type &operator()(const pair_type &v) const;
164  };
165 
166  template <typename pair_type>
167  struct sel2nd
168  {
169  const typename pair_type::second_type &operator()(const pair_type &v) const;
170  };
171  };
172 
173 
175  // IMPLEMENTATION
177 
178  template<class MagnitudeType>
179  BasicSort<MagnitudeType>::BasicSort(Teuchos::ParameterList &pl)
180  {
181  std::string which = "LM";
182  which = pl.get("Sort Strategy",which);
183  setSortType(which);
184  }
185 
186  template<class MagnitudeType>
187  BasicSort<MagnitudeType>::BasicSort(const std::string &which)
188  {
189  setSortType(which);
190  }
191 
192  template<class MagnitudeType>
194  {}
195 
196  template<class MagnitudeType>
197  void BasicSort<MagnitudeType>::setSortType(const std::string &which)
198  {
199  // make upper case
200  std::string whichlc(which);
201  std::transform(which.begin(),which.end(),whichlc.begin(),(int(*)(int)) std::toupper);
202  if (whichlc == "LM") {
203  which_ = LM;
204  }
205  else if (whichlc == "SM") {
206  which_ = SM;
207  }
208  else if (whichlc == "LR") {
209  which_ = LR;
210  }
211  else if (whichlc == "SR") {
212  which_ = SR;
213  }
214  else if (whichlc == "LI") {
215  which_ = LI;
216  }
217  else if (whichlc == "SI") {
218  which_ = SI;
219  }
220  else {
221  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Anasazi::BasicSort::setSortType(): sorting order is not valid");
222  }
223  }
224 
225  template<class MagnitudeType>
226  void BasicSort<MagnitudeType>::sort(std::vector<MagnitudeType> &evals, Teuchos::RCP<std::vector<int> > perm, int n) const
227  {
228  TEUCHOS_TEST_FOR_EXCEPTION(n < -1, std::invalid_argument, "Anasazi::BasicSort::sort(r): n must be n >= 0 or n == -1.");
229  if (n == -1) {
230  n = evals.size();
231  }
232  TEUCHOS_TEST_FOR_EXCEPTION(evals.size() < (unsigned int) n,
233  std::invalid_argument, "Anasazi::BasicSort::sort(r): eigenvalue vector size isn't consistent with n.");
234  if (perm != Teuchos::null) {
235  TEUCHOS_TEST_FOR_EXCEPTION(perm->size() < (unsigned int) n,
236  std::invalid_argument, "Anasazi::BasicSort::sort(r): permutation vector size isn't consistent with n.");
237  }
238 
239  typedef std::greater<MagnitudeType> greater_mt;
240  typedef std::less<MagnitudeType> less_mt;
241 
242  if (perm == Teuchos::null) {
243  //
244  // if permutation index is not required, just sort using the values
245  //
246  if (which_ == LM ) {
247  std::sort(evals.begin(),evals.begin()+n,compMag<greater_mt>());
248  }
249  else if (which_ == SM) {
250  std::sort(evals.begin(),evals.begin()+n,compMag<less_mt>());
251  }
252  else if (which_ == LR) {
253  std::sort(evals.begin(),evals.begin()+n,compAlg<greater_mt>());
254  }
255  else if (which_ == SR) {
256  std::sort(evals.begin(),evals.begin()+n,compAlg<less_mt>());
257  }
258  else {
259  TEUCHOS_TEST_FOR_EXCEPTION(true, SortManagerError, "Anasazi::BasicSort::sort(r): LI or SI sorting invalid for real scalar types." );
260  }
261  }
262  else {
263  //
264  // if permutation index is required, we must sort the two at once
265  // in this case, we arrange a pair structure: <value,index>
266  // default comparison operator for pair<t1,t2> is lexographic:
267  // compare first t1, then t2
268  // this works fine for us here.
269  //
270 
271  // copy the values and indices into the pair structure
272  std::vector< std::pair<MagnitudeType,int> > pairs(n);
273  for (int i=0; i<n; i++) {
274  pairs[i] = std::make_pair(evals[i],i);
275  }
276 
277  // sort the pair structure
278  if (which_ == LM) {
279  std::sort(pairs.begin(),pairs.begin()+n,compMag<greater_mt>());
280  }
281  else if (which_ == SM) {
282  std::sort(pairs.begin(),pairs.begin()+n,compMag<less_mt>());
283  }
284  else if (which_ == LR) {
285  std::sort(pairs.begin(),pairs.begin()+n,compAlg<greater_mt>());
286  }
287  else if (which_ == SR) {
288  std::sort(pairs.begin(),pairs.begin()+n,compAlg<less_mt>());
289  }
290  else {
291  TEUCHOS_TEST_FOR_EXCEPTION(true, SortManagerError, "Anasazi::BasicSort::sort(r): LI or SI sorting invalid for real scalar types." );
292  }
293 
294  // copy the values and indices out of the pair structure
295  std::transform(pairs.begin(),pairs.end(),evals.begin(),sel1st< std::pair<MagnitudeType,int> >());
296  std::transform(pairs.begin(),pairs.end(),perm->begin(),sel2nd< std::pair<MagnitudeType,int> >());
297  }
298  }
299 
300 
301  template<class T1, class T2>
302  class MakePairOp {
303  public:
304  std::pair<T1,T2> operator()( const T1 &t1, const T2 &t2 )
305  { return std::make_pair(t1, t2); }
306  };
307 
308 
309  template<class MagnitudeType>
310  void BasicSort<MagnitudeType>::sort(std::vector<MagnitudeType> &r_evals,
311  std::vector<MagnitudeType> &i_evals,
312  Teuchos::RCP< std::vector<int> > perm,
313  int n) const
314  {
315 
316  //typedef typename std::vector<MagnitudeType>::iterator r_eval_iter_t; // unused
317  //typedef typename std::vector<MagnitudeType>::iterator i_eval_iter_t; // unused
318 
319  TEUCHOS_TEST_FOR_EXCEPTION(n < -1, std::invalid_argument, "Anasazi::BasicSort::sort(r,i): n must be n >= 0 or n == -1.");
320  if (n == -1) {
321  n = r_evals.size() < i_evals.size() ? r_evals.size() : i_evals.size();
322  }
323  TEUCHOS_TEST_FOR_EXCEPTION(r_evals.size() < (unsigned int) n || i_evals.size() < (unsigned int) n,
324  std::invalid_argument, "Anasazi::BasicSort::sort(r,i): eigenvalue vector size isn't consistent with n.");
325  if (perm != Teuchos::null) {
326  TEUCHOS_TEST_FOR_EXCEPTION(perm->size() < (unsigned int) n,
327  std::invalid_argument, "Anasazi::BasicSort::sort(r,i): permutation vector size isn't consistent with n.");
328  }
329 
330  typedef std::greater<MagnitudeType> greater_mt;
331  typedef std::less<MagnitudeType> less_mt;
332 
333  //
334  // put values into pairs
335  //
336  if (perm == Teuchos::null) {
337  //
338  // not permuting, so we don't need indices in the pairs
339  //
340  std::vector< std::pair<MagnitudeType,MagnitudeType> > pairs(n);
341  // for LM,SM, the order doesn't matter
342  // for LI,SI, the imaginary goes first
343  // for LR,SR, the real goes in first
344  if (which_ == LR || which_ == SR || which_ == LM || which_ == SM) {
345  std::transform(
346  r_evals.begin(), r_evals.begin()+n,
347  i_evals.begin(), pairs.begin(),
348  MakePairOp<MagnitudeType,MagnitudeType>());
349  }
350  else {
351  std::transform(
352  i_evals.begin(), i_evals.begin()+n,
353  r_evals.begin(), pairs.begin(),
354  MakePairOp<MagnitudeType,MagnitudeType>());
355  }
356 
357  if (which_ == LR || which_ == LI) {
358  std::sort(pairs.begin(),pairs.end(),compAlg<greater_mt>());
359  }
360  else if (which_ == SR || which_ == SI) {
361  std::sort(pairs.begin(),pairs.end(),compAlg<less_mt>());
362  }
363  else if (which_ == LM) {
364  std::sort(pairs.begin(),pairs.end(),compMag2<greater_mt>());
365  }
366  else {
367  std::sort(pairs.begin(),pairs.end(),compMag2<less_mt>());
368  }
369 
370  // extract the values
371  // for LM,SM,LR,SR: order is (real,imag)
372  // for LI,SI: order is (imag,real)
373  if (which_ == LR || which_ == SR || which_ == LM || which_ == SM) {
374  std::transform(pairs.begin(),pairs.end(),r_evals.begin(),sel1st< std::pair<MagnitudeType,MagnitudeType> >());
375  std::transform(pairs.begin(),pairs.end(),i_evals.begin(),sel2nd< std::pair<MagnitudeType,MagnitudeType> >());
376  }
377  else {
378  std::transform(pairs.begin(),pairs.end(),r_evals.begin(),sel2nd< std::pair<MagnitudeType,MagnitudeType> >());
379  std::transform(pairs.begin(),pairs.end(),i_evals.begin(),sel1st< std::pair<MagnitudeType,MagnitudeType> >());
380  }
381  }
382  else {
383  //
384  // permuting, we need indices in the pairs
385  //
386  std::vector< std::pair< std::pair<MagnitudeType,MagnitudeType>, int > > pairs(n);
387  // for LM,SM, the order doesn't matter
388  // for LI,SI, the imaginary goes first
389  // for LR,SR, the real goes in first
390  if (which_ == LR || which_ == SR || which_ == LM || which_ == SM) {
391  for (int i=0; i<n; i++) {
392  pairs[i] = std::make_pair(std::make_pair(r_evals[i],i_evals[i]),i);
393  }
394  }
395  else {
396  for (int i=0; i<n; i++) {
397  pairs[i] = std::make_pair(std::make_pair(i_evals[i],r_evals[i]),i);
398  }
399  }
400 
401  if (which_ == LR || which_ == LI) {
402  std::sort(pairs.begin(),pairs.end(),compAlg<greater_mt>());
403  }
404  else if (which_ == SR || which_ == SI) {
405  std::sort(pairs.begin(),pairs.end(),compAlg<less_mt>());
406  }
407  else if (which_ == LM) {
408  std::sort(pairs.begin(),pairs.end(),compMag2<greater_mt>());
409  }
410  else {
411  std::sort(pairs.begin(),pairs.end(),compMag2<less_mt>());
412  }
413 
414  // extract the values
415  // for LM,SM,LR,SR: order is (real,imag)
416  // for LI,SI: order is (imag,real)
417  if (which_ == LR || which_ == SR || which_ == LM || which_ == SM) {
418  for (int i=0; i<n; i++) {
419  r_evals[i] = pairs[i].first.first;
420  i_evals[i] = pairs[i].first.second;
421  (*perm)[i] = pairs[i].second;
422  }
423  }
424  else {
425  for (int i=0; i<n; i++) {
426  i_evals[i] = pairs[i].first.first;
427  r_evals[i] = pairs[i].first.second;
428  (*perm)[i] = pairs[i].second;
429  }
430  }
431  }
432  }
433 
434 
435  template<class MagnitudeType>
436  template<class LTorGT>
437  bool BasicSort<MagnitudeType>::compMag<LTorGT>::operator()(MagnitudeType v1, MagnitudeType v2)
438  {
439  typedef Teuchos::ScalarTraits<MagnitudeType> MTT;
440  LTorGT comp;
441  return comp( MTT::magnitude(v1), MTT::magnitude(v2) );
442  }
443 
444  template<class MagnitudeType>
445  template<class LTorGT>
446  bool BasicSort<MagnitudeType>::compMag2<LTorGT>::operator()(std::pair<MagnitudeType,MagnitudeType> v1, std::pair<MagnitudeType,MagnitudeType> v2)
447  {
448  MagnitudeType m1 = v1.first*v1.first + v1.second*v1.second;
449  MagnitudeType m2 = v2.first*v2.first + v2.second*v2.second;
450  LTorGT comp;
451  return comp( m1, m2 );
452  }
453 
454  template<class MagnitudeType>
455  template<class LTorGT>
456  bool BasicSort<MagnitudeType>::compAlg<LTorGT>::operator()(MagnitudeType v1, MagnitudeType v2)
457  {
458  LTorGT comp;
459  return comp( v1, v2 );
460  }
461 
462  template<class MagnitudeType>
463  template<class LTorGT>
464  template<class First, class Second>
465  bool BasicSort<MagnitudeType>::compMag<LTorGT>::operator()(std::pair<First,Second> v1, std::pair<First,Second> v2) {
466  return (*this)(v1.first,v2.first);
467  }
468 
469  template<class MagnitudeType>
470  template<class LTorGT>
471  template<class First, class Second>
472  bool BasicSort<MagnitudeType>::compMag2<LTorGT>::operator()(std::pair<First,Second> v1, std::pair<First,Second> v2) {
473  return (*this)(v1.first,v2.first);
474  }
475 
476  template<class MagnitudeType>
477  template<class LTorGT>
478  template<class First, class Second>
479  bool BasicSort<MagnitudeType>::compAlg<LTorGT>::operator()(std::pair<First,Second> v1, std::pair<First,Second> v2) {
480  return (*this)(v1.first,v2.first);
481  }
482 
483  template <class MagnitudeType>
484  template <typename pair_type>
485  const typename pair_type::first_type &
487  {
488  return v.first;
489  }
490 
491  template <class MagnitudeType>
492  template <typename pair_type>
493  const typename pair_type::second_type &
495  {
496  return v.second;
497  }
498 
499 } // namespace Anasazi
500 
501 #endif // ANASAZI_BASIC_SORT_HPP
502 
BasicSort(Teuchos::ParameterList &pl)
Parameter list driven constructor.
An implementation of the Anasazi::SortManager that performs a collection of common sorting techniques...
void setSortType(const std::string &which)
Set sort type.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
SortManagerError is thrown when the Anasazi::SortManager is unable to sort the numbers, due to some failure of the sort method or error in calling it.
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...
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.
Anasazi&#39;s templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
virtual ~BasicSort()
Destructor.