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
00048 namespace Anasazi {
00049
00050 template<class ScalarType, class MV, class OP>
00051 class BasicSort : public SortManager<ScalarType,MV,OP> {
00052
00053 public:
00054
00056
00067 BasicSort( const string which = "LM" ) {
00068 setSortType(which);
00069 }
00070
00072 virtual ~BasicSort() {};
00073
00075
00086 void setSortType( const string which ) {
00087 which_ = which;
00088 TEST_FOR_EXCEPTION(which_.compare("LM") && which_.compare("SM") &&
00089 which_.compare("LR") && which_.compare("SR") &&
00090 which_.compare("LI") && which_.compare("SI"), std::invalid_argument,
00091 "Anasazi::BasicSort::sort(): sorting order is not valid");
00092 };
00093
00095
00104 void sort(Eigensolver<ScalarType,MV,OP>* solver, const int n, std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &evals, std::vector<int> *perm = 0) const;
00105
00122 void sort(Eigensolver<ScalarType,MV,OP>* solver,
00123 const int n,
00124 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &r_evals,
00125 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &i_evals,
00126 std::vector<int> *perm = 0) const;
00127
00128 protected:
00129
00131
00141 string which_;
00142
00143 };
00144
00145 template<class ScalarType, class MV, class OP>
00146 void BasicSort<ScalarType,MV,OP>::sort(Eigensolver<ScalarType,MV,OP>* solver, const int n, std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &evals, std::vector<int> *perm) const
00147 {
00148 int i=0,j=0;
00149
00150 TEST_FOR_EXCEPTION(evals.size() < (unsigned int) n,
00151 std::invalid_argument, "Anasazi::BasicSort:sort(): eigenvalue vector size isn't consistent with n.");
00152 if (perm) {
00153 TEST_FOR_EXCEPTION(perm->size() < (unsigned int) n,
00154 std::invalid_argument, "Anasazi::BasicSort:sort(): permutation vector size isn't consistent with n.");
00155 }
00156
00157
00158 int tempord=0;
00159
00160 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00161 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00162
00163
00164 MagnitudeType temp;
00165
00166 Teuchos::LAPACK<int,MagnitudeType> lapack;
00167
00168
00169
00170
00171 if (perm) {
00172 for (i=0; i < n; i++) {
00173 (*perm)[i] = i;
00174 }
00175 }
00176
00177
00178
00179
00180
00181 if (!which_.compare("SM")) {
00182 for (j=1; j < n; j++) {
00183 temp = evals[j];
00184 if (perm) {
00185 tempord = (*perm)[j];
00186 }
00187 MagnitudeType temp2 = MT::magnitude(evals[j]);
00188 for (i=j-1; i >=0 && MT::magnitude(evals[i]) > temp2; i--) {
00189 evals[i+1] = evals[i];
00190 if (perm) {
00191 (*perm)[i+1]=(*perm)[i];
00192 }
00193 }
00194 evals[i+1] = temp;
00195 if (perm) {
00196 (*perm)[i+1] = tempord;
00197 }
00198 }
00199 return;
00200 }
00201
00202
00203
00204 if (!which_.compare("SR")) {
00205 for (j=1; j < n; j++) {
00206 temp = evals[j];
00207 if (perm) {
00208 tempord = (*perm)[j];
00209 }
00210 for (i=j-1; i >= 0 && evals[i] > temp; i--) {
00211 evals[i+1]=evals[i];
00212 if (perm) {
00213 (*perm)[i+1]=(*perm)[i];
00214 }
00215 }
00216 evals[i+1] = temp;
00217 if (perm) {
00218 (*perm)[i+1] = tempord;
00219 }
00220 }
00221 return;
00222 }
00223
00224
00225
00226
00227
00228 TEST_FOR_EXCEPTION(!which_.compare("SI"), SortManagerError,
00229 "Anasazi::BasicSort::sort() with one arg assumes real eigenvalues");
00230
00231
00232
00233 if (!which_.compare("LM")) {
00234 for (j=1; j < n; j++) {
00235 temp = evals[j];
00236 if (perm) {
00237 tempord = (*perm)[j];
00238 }
00239 MagnitudeType temp2 = MT::magnitude(evals[j]);
00240 for (i=j-1; i >= 0 && MT::magnitude(evals[i]) < temp2; i--) {
00241 evals[i+1]=evals[i];
00242 if (perm) {
00243 (*perm)[i+1]=(*perm)[i];
00244 }
00245 }
00246 evals[i+1] = temp;
00247 if (perm) {
00248 (*perm)[i+1] = tempord;
00249 }
00250 }
00251 return;
00252 }
00253
00254
00255
00256 if (!which_.compare("LR")) {
00257 for (j=1; j < n; j++) {
00258 temp = evals[j];
00259 if (perm) {
00260 tempord = (*perm)[j];
00261 }
00262 for (i=j-1; i >= 0 && evals[i]<temp; i--) {
00263 evals[i+1]=evals[i];
00264 if (perm) {
00265 (*perm)[i+1]=(*perm)[i];
00266 }
00267 }
00268 evals[i+1] = temp;
00269 if (perm) {
00270 (*perm)[i+1] = tempord;
00271 }
00272 }
00273 return;
00274 }
00275
00276
00277
00278
00279
00280 TEST_FOR_EXCEPTION(!which_.compare("LI"), SortManagerError,
00281 "Anasazi::BasicSort::sort() with one arg assumes real eigenvalues");
00282
00283
00284 TEST_FOR_EXCEPTION(true, std::logic_error,
00285 "Anasazi::BasicSort::sort(): sorting order is not valid");
00286 }
00287
00288
00289 template<class ScalarType, class MV, class OP>
00290 void BasicSort<ScalarType,MV,OP>::sort(Eigensolver<ScalarType,MV,OP>* solver,
00291 const int n,
00292 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &r_evals,
00293 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &i_evals,
00294 std::vector<int> *perm) const
00295 {
00296 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00297 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00298
00299 TEST_FOR_EXCEPTION(r_evals.size() < (unsigned int) n || i_evals.size() < (unsigned int) n,
00300 std::invalid_argument, "Anasazi::BasicSort:sort(): real and imaginary vector sizes aren't consistent with n.");
00301 if (perm) {
00302 TEST_FOR_EXCEPTION(perm->size() < (unsigned int) n,
00303 std::invalid_argument, "Anasazi::BasicSort:sort(): permutation vector size isn't consistent with n.");
00304 }
00305 int i=0,j=0;
00306 int tempord=0;
00307
00308 MagnitudeType temp, tempr, tempi;
00309 Teuchos::LAPACK<int,MagnitudeType> lapack;
00310
00311
00312
00313 if (perm) {
00314 for (i=0; i < n; i++) {
00315 (*perm)[i] = i;
00316 }
00317 }
00318
00319
00320
00321
00322
00323 if (!which_.compare("SM")) {
00324 for (j=1; j < n; j++) {
00325 tempr = r_evals[j]; tempi = i_evals[j];
00326 if (perm) {
00327 tempord = (*perm)[j];
00328 }
00329 temp=lapack.LAPY2(r_evals[j],i_evals[j]);
00330 for (i=j-1; i>=0 && lapack.LAPY2(r_evals[i],i_evals[i]) > temp; i--) {
00331 r_evals[i+1]=r_evals[i]; i_evals[i+1]=i_evals[i];
00332 if (perm) {
00333 (*perm)[i+1]=(*perm)[i];
00334 }
00335 }
00336 r_evals[i+1] = tempr; i_evals[i+1] = tempi;
00337 if (perm) {
00338 (*perm)[i+1] = tempord;
00339 }
00340 }
00341 return;
00342 }
00343
00344
00345
00346 if (!which_.compare("SR")) {
00347 for (j=1; j < n; j++) {
00348 tempr = r_evals[j]; tempi = i_evals[j];
00349 if (perm) {
00350 tempord = (*perm)[j];
00351 }
00352 for (i=j-1; i>=0 && r_evals[i]>tempr; i--) {
00353 r_evals[i+1]=r_evals[i]; i_evals[i+1]=i_evals[i];
00354 if (perm) {
00355 (*perm)[i+1]=(*perm)[i];
00356 }
00357 }
00358 r_evals[i+1] = tempr; i_evals[i+1] = tempi;
00359 if (perm) {
00360 (*perm)[i+1] = tempord;
00361 }
00362 }
00363 return;
00364 }
00365
00366
00367
00368 if (!which_.compare("SI")) {
00369 for (j=1; j < n; j++) {
00370 tempr = r_evals[j]; tempi = i_evals[j];
00371 if (perm) {
00372 tempord = (*perm)[j];
00373 }
00374 for (i=j-1; i>=0 && i_evals[i]>tempi; i--) {
00375 r_evals[i+1]=r_evals[i]; i_evals[i+1]=i_evals[i];
00376 if (perm) {
00377 (*perm)[i+1]=(*perm)[i];
00378 }
00379 }
00380 r_evals[i+1] = tempr; i_evals[i+1] = tempi;
00381 if (perm) {
00382 (*perm)[i+1] = tempord;
00383 }
00384 }
00385 return;
00386 }
00387
00388
00389
00390 if (!which_.compare("LM")) {
00391 for (j=1; j < n; j++) {
00392 tempr = r_evals[j]; tempi = i_evals[j];
00393 if (perm) {
00394 tempord = (*perm)[j];
00395 }
00396 temp=lapack.LAPY2(r_evals[j],i_evals[j]);
00397 for (i=j-1; i>=0 && lapack.LAPY2(r_evals[i],i_evals[i])<temp; i--) {
00398 r_evals[i+1]=r_evals[i]; i_evals[i+1]=i_evals[i];
00399 if (perm) {
00400 (*perm)[i+1]=(*perm)[i];
00401 }
00402 }
00403 r_evals[i+1] = tempr; i_evals[i+1] = tempi;
00404 if (perm) {
00405 (*perm)[i+1] = tempord;
00406 }
00407 }
00408 return;
00409 }
00410
00411
00412
00413 if (!which_.compare("LR")) {
00414 for (j=1; j < n; j++) {
00415 tempr = r_evals[j]; tempi = i_evals[j];
00416 if (perm) {
00417 tempord = (*perm)[j];
00418 }
00419 for (i=j-1; i>=0 && r_evals[i]<tempr; i--) {
00420 r_evals[i+1]=r_evals[i]; i_evals[i+1]=i_evals[i];
00421 if (perm) {
00422 (*perm)[i+1]=(*perm)[i];
00423 }
00424 }
00425 r_evals[i+1] = tempr; i_evals[i+1] = tempi;
00426 if (perm) {
00427 (*perm)[i+1] = tempord;
00428 }
00429 }
00430 return;
00431 }
00432
00433
00434
00435 if (!which_.compare("LI")) {
00436 for (j=1; j < n; j++) {
00437 tempr = r_evals[j]; tempi = i_evals[j];
00438 if (perm) {
00439 tempord = (*perm)[j];
00440 }
00441 for (i=j-1; i>=0 && i_evals[i]<tempi; i--) {
00442 r_evals[i+1]=r_evals[i]; i_evals[i+1]=i_evals[i];
00443 if (perm) {
00444 (*perm)[i+1]=(*perm)[i];
00445 }
00446 }
00447 r_evals[i+1] = tempr; i_evals[i+1] = tempi;
00448 if (perm) {
00449 (*perm)[i+1] = tempord;
00450 }
00451 }
00452 return;
00453 }
00454
00455 TEST_FOR_EXCEPTION(true, std::logic_error,
00456 "Anasazi::BasicSort::sort(): sorting order is not valid");
00457 }
00458
00459
00460 }
00461
00462 #endif // ANASAZI_BASIC_SORT_HPP
00463