Anasazi Version of the Day
Loading...
Searching...
No Matches
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 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40// @HEADER
41
46#ifndef ANASAZI_BASIC_SORT_HPP
47#define ANASAZI_BASIC_SORT_HPP
48
56#include "AnasaziConfigDefs.hpp"
58#include "Teuchos_LAPACK.hpp"
59#include "Teuchos_ScalarTraits.hpp"
60#include "Teuchos_ParameterList.hpp"
61
62namespace Anasazi {
63
64 template<class MagnitudeType>
65 class BasicSort : public SortManager<MagnitudeType> {
66
67 public:
68
74 BasicSort( Teuchos::ParameterList &pl );
75
80 BasicSort( const std::string &which = "LM" );
81
83 virtual ~BasicSort();
84
86
95 void setSortType( const std::string &which );
96
111 void sort(std::vector<MagnitudeType> &evals, Teuchos::RCP<std::vector<int> > perm = Teuchos::null, int n = -1) const;
112
131 void sort(std::vector<MagnitudeType> &r_evals,
132 std::vector<MagnitudeType> &i_evals,
133 Teuchos::RCP<std::vector<int> > perm = Teuchos::null,
134 int n = -1) const;
135
136 protected:
137
138 // enum for sort type
139 enum SType {
140 LM, SM,
141 LR, SR,
142 LI, SI
143 };
144 SType which_;
145
146 // sorting methods
147 template <class LTorGT>
148 struct compMag {
149 // for real-only LM,SM
150 bool operator()(MagnitudeType, MagnitudeType);
151 // for real-only LM,SM with permutation
152 template <class First, class Second>
153 bool operator()(std::pair<First,Second>, std::pair<First,Second>);
154 };
155
156 template <class LTorGT>
157 struct compMag2 {
158 // for real-imag LM,SM
159 bool operator()(std::pair<MagnitudeType,MagnitudeType>, std::pair<MagnitudeType,MagnitudeType>);
160 // for real-imag LM,SM with permutation
161 template <class First, class Second>
162 bool operator()(std::pair<First,Second>, std::pair<First,Second>);
163 };
164
165 template <class LTorGT>
166 struct compAlg {
167 // for real-imag LR,SR,LI,SI
168 bool operator()(MagnitudeType, MagnitudeType);
169 template <class First, class Second>
170 bool operator()(std::pair<First,Second>, std::pair<First,Second>);
171 };
172
173 template <typename pair_type>
174 struct sel1st
175 {
176 const typename pair_type::first_type &operator()(const pair_type &v) const;
177 };
178
179 template <typename pair_type>
180 struct sel2nd
181 {
182 const typename pair_type::second_type &operator()(const pair_type &v) const;
183 };
184 };
185
186
188 // IMPLEMENTATION
190
191 template<class MagnitudeType>
192 BasicSort<MagnitudeType>::BasicSort(Teuchos::ParameterList &pl)
193 {
194 std::string which = "LM";
195 which = pl.get("Sort Strategy",which);
196 setSortType(which);
197 }
198
199 template<class MagnitudeType>
200 BasicSort<MagnitudeType>::BasicSort(const std::string &which)
201 {
202 setSortType(which);
203 }
204
205 template<class MagnitudeType>
208
209 template<class MagnitudeType>
210 void BasicSort<MagnitudeType>::setSortType(const std::string &which)
211 {
212 // make upper case
213 std::string whichlc(which);
214 std::transform(which.begin(),which.end(),whichlc.begin(),(int(*)(int)) std::toupper);
215 if (whichlc == "LM") {
216 which_ = LM;
217 }
218 else if (whichlc == "SM") {
219 which_ = SM;
220 }
221 else if (whichlc == "LR") {
222 which_ = LR;
223 }
224 else if (whichlc == "SR") {
225 which_ = SR;
226 }
227 else if (whichlc == "LI") {
228 which_ = LI;
229 }
230 else if (whichlc == "SI") {
231 which_ = SI;
232 }
233 else {
234 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Anasazi::BasicSort::setSortType(): sorting order is not valid");
235 }
236 }
237
238 template<class MagnitudeType>
239 void BasicSort<MagnitudeType>::sort(std::vector<MagnitudeType> &evals, Teuchos::RCP<std::vector<int> > perm, int n) const
240 {
241 TEUCHOS_TEST_FOR_EXCEPTION(n < -1, std::invalid_argument, "Anasazi::BasicSort::sort(r): n must be n >= 0 or n == -1.");
242 if (n == -1) {
243 n = evals.size();
244 }
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.");
250 }
251
252 typedef std::greater<MagnitudeType> greater_mt;
253 typedef std::less<MagnitudeType> less_mt;
254
255 if (perm == Teuchos::null) {
256 //
257 // if permutation index is not required, just sort using the values
258 //
259 if (which_ == LM ) {
260 std::sort(evals.begin(),evals.begin()+n,compMag<greater_mt>());
261 }
262 else if (which_ == SM) {
263 std::sort(evals.begin(),evals.begin()+n,compMag<less_mt>());
264 }
265 else if (which_ == LR) {
266 std::sort(evals.begin(),evals.begin()+n,compAlg<greater_mt>());
267 }
268 else if (which_ == SR) {
269 std::sort(evals.begin(),evals.begin()+n,compAlg<less_mt>());
270 }
271 else {
272 TEUCHOS_TEST_FOR_EXCEPTION(true, SortManagerError, "Anasazi::BasicSort::sort(r): LI or SI sorting invalid for real scalar types." );
273 }
274 }
275 else {
276 //
277 // if permutation index is required, we must sort the two at once
278 // in this case, we arrange a pair structure: <value,index>
279 // default comparison operator for pair<t1,t2> is lexographic:
280 // compare first t1, then t2
281 // this works fine for us here.
282 //
283
284 // copy the values and indices into the pair structure
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);
288 }
289
290 // sort the pair structure
291 if (which_ == LM) {
292 std::sort(pairs.begin(),pairs.begin()+n,compMag<greater_mt>());
293 }
294 else if (which_ == SM) {
295 std::sort(pairs.begin(),pairs.begin()+n,compMag<less_mt>());
296 }
297 else if (which_ == LR) {
298 std::sort(pairs.begin(),pairs.begin()+n,compAlg<greater_mt>());
299 }
300 else if (which_ == SR) {
301 std::sort(pairs.begin(),pairs.begin()+n,compAlg<less_mt>());
302 }
303 else {
304 TEUCHOS_TEST_FOR_EXCEPTION(true, SortManagerError, "Anasazi::BasicSort::sort(r): LI or SI sorting invalid for real scalar types." );
305 }
306
307 // copy the values and indices out of the pair structure
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> >());
310 }
311 }
312
313
314 template<class T1, class T2>
315 class MakePairOp {
316 public:
317 std::pair<T1,T2> operator()( const T1 &t1, const T2 &t2 )
318 { return std::make_pair(t1, t2); }
319 };
320
321
322 template<class MagnitudeType>
323 void BasicSort<MagnitudeType>::sort(std::vector<MagnitudeType> &r_evals,
324 std::vector<MagnitudeType> &i_evals,
325 Teuchos::RCP< std::vector<int> > perm,
326 int n) const
327 {
328
329 //typedef typename std::vector<MagnitudeType>::iterator r_eval_iter_t; // unused
330 //typedef typename std::vector<MagnitudeType>::iterator i_eval_iter_t; // unused
331
332 TEUCHOS_TEST_FOR_EXCEPTION(n < -1, std::invalid_argument, "Anasazi::BasicSort::sort(r,i): n must be n >= 0 or n == -1.");
333 if (n == -1) {
334 n = r_evals.size() < i_evals.size() ? r_evals.size() : i_evals.size();
335 }
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.");
341 }
342
343 typedef std::greater<MagnitudeType> greater_mt;
344 typedef std::less<MagnitudeType> less_mt;
345
346 //
347 // put values into pairs
348 //
349 if (perm == Teuchos::null) {
350 //
351 // not permuting, so we don't need indices in the pairs
352 //
353 std::vector< std::pair<MagnitudeType,MagnitudeType> > pairs(n);
354 // for LM,SM, the order doesn't matter
355 // for LI,SI, the imaginary goes first
356 // for LR,SR, the real goes in first
357 if (which_ == LR || which_ == SR || which_ == LM || which_ == SM) {
358 std::transform(
359 r_evals.begin(), r_evals.begin()+n,
360 i_evals.begin(), pairs.begin(),
361 MakePairOp<MagnitudeType,MagnitudeType>());
362 }
363 else {
364 std::transform(
365 i_evals.begin(), i_evals.begin()+n,
366 r_evals.begin(), pairs.begin(),
367 MakePairOp<MagnitudeType,MagnitudeType>());
368 }
369
370 if (which_ == LR || which_ == LI) {
371 std::sort(pairs.begin(),pairs.end(),compAlg<greater_mt>());
372 }
373 else if (which_ == SR || which_ == SI) {
374 std::sort(pairs.begin(),pairs.end(),compAlg<less_mt>());
375 }
376 else if (which_ == LM) {
377 std::sort(pairs.begin(),pairs.end(),compMag2<greater_mt>());
378 }
379 else {
380 std::sort(pairs.begin(),pairs.end(),compMag2<less_mt>());
381 }
382
383 // extract the values
384 // for LM,SM,LR,SR: order is (real,imag)
385 // for LI,SI: order is (imag,real)
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> >());
389 }
390 else {
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> >());
393 }
394 }
395 else {
396 //
397 // permuting, we need indices in the pairs
398 //
399 std::vector< std::pair< std::pair<MagnitudeType,MagnitudeType>, int > > pairs(n);
400 // for LM,SM, the order doesn't matter
401 // for LI,SI, the imaginary goes first
402 // for LR,SR, the real goes in first
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);
406 }
407 }
408 else {
409 for (int i=0; i<n; i++) {
410 pairs[i] = std::make_pair(std::make_pair(i_evals[i],r_evals[i]),i);
411 }
412 }
413
414 if (which_ == LR || which_ == LI) {
415 std::sort(pairs.begin(),pairs.end(),compAlg<greater_mt>());
416 }
417 else if (which_ == SR || which_ == SI) {
418 std::sort(pairs.begin(),pairs.end(),compAlg<less_mt>());
419 }
420 else if (which_ == LM) {
421 std::sort(pairs.begin(),pairs.end(),compMag2<greater_mt>());
422 }
423 else {
424 std::sort(pairs.begin(),pairs.end(),compMag2<less_mt>());
425 }
426
427 // extract the values
428 // for LM,SM,LR,SR: order is (real,imag)
429 // for LI,SI: order is (imag,real)
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;
435 }
436 }
437 else {
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;
442 }
443 }
444 }
445 }
446
447
448 template<class MagnitudeType>
449 template<class LTorGT>
450 bool BasicSort<MagnitudeType>::compMag<LTorGT>::operator()(MagnitudeType v1, MagnitudeType v2)
451 {
452 typedef Teuchos::ScalarTraits<MagnitudeType> MTT;
453 LTorGT comp;
454 return comp( MTT::magnitude(v1), MTT::magnitude(v2) );
455 }
456
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)
460 {
461 MagnitudeType m1 = v1.first*v1.first + v1.second*v1.second;
462 MagnitudeType m2 = v2.first*v2.first + v2.second*v2.second;
463 LTorGT comp;
464 return comp( m1, m2 );
465 }
466
467 template<class MagnitudeType>
468 template<class LTorGT>
469 bool BasicSort<MagnitudeType>::compAlg<LTorGT>::operator()(MagnitudeType v1, MagnitudeType v2)
470 {
471 LTorGT comp;
472 return comp( v1, v2 );
473 }
474
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);
480 }
481
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);
487 }
488
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);
494 }
495
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
500 {
501 return v.first;
502 }
503
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
508 {
509 return v.second;
510 }
511
512} // namespace Anasazi
513
514#endif // ANASAZI_BASIC_SORT_HPP
515
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.