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
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "Teuchos_LAPACK.hpp"
00043 #include "Teuchos_ParameterList.hpp"
00044 #include "LOCA_GlobalData.H"
00045 #include "LOCA_EigenvalueSort_Strategies.H"
00046
00047 NOX::Abstract::Group::ReturnType
00048 LOCA::EigenvalueSort::LargestMagnitude::sort(int n, double* evals,
00049 std::vector<int>* perm) const
00050 {
00051 int i, j, tempord;
00052 double temp, temp2;
00053 Teuchos::LAPACK<int,double> lapack;
00054
00055
00056 if (perm) {
00057 for (i=0; i < n; i++) {
00058 (*perm)[i] = i;
00059 }
00060 }
00061
00062
00063
00064
00065 for (j=1; j < n; ++j) {
00066 temp = evals[j];
00067 if (perm)
00068 tempord = (*perm)[j];
00069 temp2 = evals[j]*evals[j];
00070 for (i=j-1; i>=0 && (evals[i]*evals[i])<temp2; --i) {
00071 evals[i+1]=evals[i];
00072 if (perm)
00073 (*perm)[i+1]=(*perm)[i];
00074 }
00075 evals[i+1] = temp;
00076 if (perm)
00077 (*perm)[i+1] = tempord;
00078 }
00079 return NOX::Abstract::Group::Ok;
00080 }
00081
00082 NOX::Abstract::Group::ReturnType
00083 LOCA::EigenvalueSort::LargestMagnitude::sort(int n, double* r_evals,
00084 double* i_evals,
00085 std::vector<int>* perm) const
00086 {
00087 int i, j, tempord;
00088 double temp, tempr, tempi;
00089 Teuchos::LAPACK<int,double> lapack;
00090
00091
00092
00093 if (perm) {
00094 for (i=0; i < n; i++) {
00095 (*perm)[i] = i;
00096 }
00097 }
00098
00099
00100
00101
00102 for (j=1; j < n; ++j) {
00103 tempr = r_evals[j]; tempi = i_evals[j];
00104 if (perm)
00105 tempord = (*perm)[j];
00106 temp=lapack.LAPY2(r_evals[j],i_evals[j]);
00107 for (i=j-1; i>=0 && lapack.LAPY2(r_evals[i],i_evals[i])<temp; --i) {
00108 r_evals[i+1]=r_evals[i]; i_evals[i+1]=i_evals[i];
00109 if (perm)
00110 (*perm)[i+1]=(*perm)[i];
00111 }
00112 r_evals[i+1] = tempr; i_evals[i+1] = tempi;
00113 if (perm)
00114 (*perm)[i+1] = tempord;
00115 }
00116 return NOX::Abstract::Group::Ok;
00117 }
00118
00119 NOX::Abstract::Group::ReturnType
00120 LOCA::EigenvalueSort::SmallestMagnitude::sort(int n, double* evals,
00121 std::vector<int>* perm) const
00122 {
00123 int i, j, tempord;
00124 double temp, temp2;
00125 Teuchos::LAPACK<int,double> lapack;
00126
00127
00128 if (perm) {
00129 for (i=0; i < n; i++) {
00130 (*perm)[i] = i;
00131 }
00132 }
00133
00134
00135
00136
00137 for (j=1; j < n; ++j) {
00138 temp = evals[j];
00139 if (perm)
00140 tempord = (*perm)[j];
00141 temp2 = evals[j]*evals[j];
00142 for (i=j-1; i>=0 && (evals[i]*evals[i])>temp2; --i) {
00143 evals[i+1]=evals[i];
00144 if (perm)
00145 (*perm)[i+1]=(*perm)[i];
00146 }
00147 evals[i+1] = temp;
00148 if (perm)
00149 (*perm)[i+1] = tempord;
00150 }
00151 return NOX::Abstract::Group::Ok;
00152 }
00153
00154 NOX::Abstract::Group::ReturnType
00155 LOCA::EigenvalueSort::SmallestMagnitude::sort(int n, double* r_evals,
00156 double* i_evals,
00157 std::vector<int>* perm) const
00158 {
00159 int i, j, tempord;
00160 double temp, tempr, tempi;
00161 Teuchos::LAPACK<int,double> lapack;
00162
00163
00164
00165 if (perm) {
00166 for (i=0; i < n; i++) {
00167 (*perm)[i] = i;
00168 }
00169 }
00170
00171
00172
00173
00174 for (j=1; j < n; ++j) {
00175 tempr = r_evals[j]; tempi = i_evals[j];
00176 if (perm)
00177 tempord = (*perm)[j];
00178 temp=lapack.LAPY2(r_evals[j],i_evals[j]);
00179 for (i=j-1; i>=0 && lapack.LAPY2(r_evals[i],i_evals[i])>temp; --i) {
00180 r_evals[i+1]=r_evals[i]; i_evals[i+1]=i_evals[i];
00181 if (perm)
00182 (*perm)[i+1]=(*perm)[i];
00183 }
00184 r_evals[i+1] = tempr; i_evals[i+1] = tempi;
00185 if (perm)
00186 (*perm)[i+1] = tempord;
00187 }
00188 return NOX::Abstract::Group::Ok;
00189 }
00190
00191
00192 NOX::Abstract::Group::ReturnType
00193 LOCA::EigenvalueSort::LargestReal::sort(int n, double* evals,
00194 std::vector<int>* perm) const
00195 {
00196 int i, j, tempord;
00197 double temp;
00198 Teuchos::LAPACK<int,double> lapack;
00199
00200
00201 if (perm) {
00202 for (i=0; i < n; i++) {
00203 (*perm)[i] = i;
00204 }
00205 }
00206
00207
00208
00209
00210 for (j=1; j < n; ++j) {
00211 temp = evals[j];
00212 if (perm)
00213 tempord = (*perm)[j];
00214 for (i=j-1; i>=0 && evals[i]<temp; --i) {
00215 evals[i+1]=evals[i];
00216 if (perm)
00217 (*perm)[i+1]=(*perm)[i];
00218 }
00219 evals[i+1] = temp;
00220 if (perm)
00221 (*perm)[i+1] = tempord;
00222 }
00223 return NOX::Abstract::Group::Ok;
00224 }
00225
00226 NOX::Abstract::Group::ReturnType
00227 LOCA::EigenvalueSort::LargestReal::sort(int n, double* r_evals,
00228 double* i_evals,
00229 std::vector<int>* perm) const
00230 {
00231 int i, j, tempord;
00232 double tempr, tempi;
00233 Teuchos::LAPACK<int,double> lapack;
00234
00235
00236 if (perm) {
00237 for (i=0; i < n; i++) {
00238 (*perm)[i] = i;
00239 }
00240 }
00241
00242
00243
00244
00245 for (j=1; j < n; ++j) {
00246 tempr = r_evals[j]; tempi = i_evals[j];
00247 if (perm)
00248 tempord = (*perm)[j];
00249 for (i=j-1; i>=0 && r_evals[i]<tempr; --i) {
00250 r_evals[i+1]=r_evals[i]; i_evals[i+1]=i_evals[i];
00251 if (perm)
00252 (*perm)[i+1]=(*perm)[i];
00253 }
00254 r_evals[i+1] = tempr; i_evals[i+1] = tempi;
00255 if (perm)
00256 (*perm)[i+1] = tempord;
00257 }
00258 return NOX::Abstract::Group::Ok;
00259 }
00260
00261 NOX::Abstract::Group::ReturnType
00262 LOCA::EigenvalueSort::SmallestReal::sort(int n, double* evals,
00263 std::vector<int>* perm) const
00264 {
00265 int i, j, tempord;
00266 double temp;
00267 Teuchos::LAPACK<int,double> lapack;
00268
00269
00270 if (perm) {
00271 for (i=0; i < n; i++) {
00272 (*perm)[i] = i;
00273 }
00274 }
00275
00276
00277
00278
00279 for (j=1; j < n; ++j) {
00280 temp = evals[j];
00281 if (perm)
00282 tempord = (*perm)[j];
00283 for (i=j-1; i>=0 && evals[i]>temp; --i) {
00284 evals[i+1]=evals[i];
00285 if (perm)
00286 (*perm)[i+1]=(*perm)[i];
00287 }
00288 evals[i+1] = temp;
00289 if (perm)
00290 (*perm)[i+1] = tempord;
00291 }
00292 return NOX::Abstract::Group::Ok;
00293 }
00294
00295 NOX::Abstract::Group::ReturnType
00296 LOCA::EigenvalueSort::SmallestReal::sort(int n, double* r_evals,
00297 double* i_evals,
00298 std::vector<int>* perm) const
00299 {
00300 int i, j, tempord;
00301 double tempr, tempi;
00302 Teuchos::LAPACK<int,double> lapack;
00303
00304
00305 if (perm) {
00306 for (i=0; i < n; i++) {
00307 (*perm)[i] = i;
00308 }
00309 }
00310
00311
00312
00313
00314 for (j=1; j < n; ++j) {
00315 tempr = r_evals[j]; tempi = i_evals[j];
00316 if (perm)
00317 tempord = (*perm)[j];
00318 for (i=j-1; i>=0 && r_evals[i]>tempr; --i) {
00319 r_evals[i+1]=r_evals[i]; i_evals[i+1]=i_evals[i];
00320 if (perm)
00321 (*perm)[i+1]=(*perm)[i];
00322 }
00323 r_evals[i+1] = tempr; i_evals[i+1] = tempi;
00324 if (perm)
00325 (*perm)[i+1] = tempord;
00326 }
00327 return NOX::Abstract::Group::Ok;
00328 }
00329
00330 NOX::Abstract::Group::ReturnType
00331 LOCA::EigenvalueSort::LargestImaginary::sort(int n, double* evals,
00332 std::vector<int>* perm) const
00333 {
00334 return NOX::Abstract::Group::NotDefined;
00335 }
00336
00337 NOX::Abstract::Group::ReturnType
00338 LOCA::EigenvalueSort::LargestImaginary::sort(int n, double* r_evals,
00339 double* i_evals,
00340 std::vector<int>* perm) const
00341 {
00342 int i, j, tempord;
00343 double tempr, tempi;
00344 Teuchos::LAPACK<int,double> lapack;
00345
00346
00347 if (perm) {
00348 for (i=0; i < n; i++) {
00349 (*perm)[i] = i;
00350 }
00351 }
00352
00353
00354
00355
00356 for (j=1; j < n; ++j) {
00357 tempr = r_evals[j]; tempi = i_evals[j];
00358 if (perm)
00359 tempord = (*perm)[j];
00360 for (i=j-1; i>=0 && i_evals[i]<tempi; --i) {
00361 r_evals[i+1]=r_evals[i]; i_evals[i+1]=i_evals[i];
00362 if (perm)
00363 (*perm)[i+1]=(*perm)[i];
00364 }
00365 r_evals[i+1] = tempr; i_evals[i+1] = tempi;
00366 if (perm)
00367 (*perm)[i+1] = tempord;
00368 }
00369 return NOX::Abstract::Group::Ok;
00370 }
00371
00372 NOX::Abstract::Group::ReturnType
00373 LOCA::EigenvalueSort::SmallestImaginary::sort(int n, double* evals,
00374 std::vector<int>* perm) const
00375 {
00376 return NOX::Abstract::Group::NotDefined;
00377 }
00378
00379 NOX::Abstract::Group::ReturnType
00380 LOCA::EigenvalueSort::SmallestImaginary::sort(int n, double* r_evals,
00381 double* i_evals,
00382 std::vector<int>* perm) const
00383 {
00384 int i, j, tempord;
00385 double tempr, tempi;
00386 Teuchos::LAPACK<int,double> lapack;
00387
00388
00389 if (perm) {
00390 for (i=0; i < n; i++) {
00391 (*perm)[i] = i;
00392 }
00393 }
00394
00395
00396
00397
00398 for (j=1; j < n; ++j) {
00399 tempr = r_evals[j]; tempi = i_evals[j];
00400 if (perm)
00401 tempord = (*perm)[j];
00402 for (i=j-1; i>=0 && i_evals[i]>tempi; --i) {
00403 r_evals[i+1]=r_evals[i]; i_evals[i+1]=i_evals[i];
00404 if (perm)
00405 (*perm)[i+1]=(*perm)[i];
00406 }
00407 r_evals[i+1] = tempr; i_evals[i+1] = tempi;
00408 if (perm)
00409 (*perm)[i+1] = tempord;
00410 }
00411 return NOX::Abstract::Group::Ok;
00412 }
00413
00414 LOCA::EigenvalueSort::LargestRealInverseCayley::LargestRealInverseCayley(
00415 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00416 const Teuchos::RCP<Teuchos::ParameterList>& eigenParams) :
00417 sigma(0.0),
00418 mu(0.0)
00419 {
00420 sigma = eigenParams->get("Cayley Pole",0.0);
00421 mu = eigenParams->get("Cayley Zero",0.0);
00422 }
00423
00424 NOX::Abstract::Group::ReturnType
00425 LOCA::EigenvalueSort::LargestRealInverseCayley::sort(
00426 int n, double* evals,
00427 std::vector<int>* perm) const
00428 {
00429 int i, j, tempord;
00430 double temp, templambda;
00431 Teuchos::LAPACK<int,double> lapack;
00432
00433
00434 if (perm) {
00435 for (i=0; i < n; i++) {
00436 (*perm)[i] = i;
00437 }
00438 }
00439
00440
00441
00442
00443
00444 for (j=1; j < n; ++j) {
00445 temp = evals[j];
00446 tempord = (*perm)[j];
00447 templambda=realLambda(evals[j],0);
00448 for (i=j-1; i>=0 && realLambda(evals[i],0)<templambda; --i) {
00449 evals[i+1]=evals[i];
00450 (*perm)[i+1]=(*perm)[i];
00451 }
00452 evals[i+1] = temp; (*perm)[i+1] = tempord;
00453 }
00454 return NOX::Abstract::Group::Ok;
00455 }
00456
00457 NOX::Abstract::Group::ReturnType
00458 LOCA::EigenvalueSort::LargestRealInverseCayley::sort(
00459 int n, double* r_evals,
00460 double* i_evals,
00461 std::vector<int>* perm) const
00462 {
00463 int i, j, tempord;
00464 double temp, tempr, tempi;
00465 Teuchos::LAPACK<int,double> lapack;
00466
00467
00468 if (perm) {
00469 for (i=0; i < n; i++) {
00470 (*perm)[i] = i;
00471 }
00472 }
00473
00474
00475
00476
00477
00478 for (j=1; j < n; ++j) {
00479 tempr = r_evals[j]; tempi = i_evals[j];
00480 tempord = (*perm)[j];
00481 temp=realLambda(r_evals[j],i_evals[j]);
00482 for (i=j-1; i>=0 && realLambda(r_evals[i],i_evals[i])<temp; --i) {
00483 r_evals[i+1]=r_evals[i]; i_evals[i+1]=i_evals[i];
00484 (*perm)[i+1]=(*perm)[i];
00485 }
00486 r_evals[i+1] = tempr; i_evals[i+1] = tempi; (*perm)[i+1] = tempord;
00487 }
00488 return NOX::Abstract::Group::Ok;
00489 }
00490
00491 double
00492 LOCA::EigenvalueSort::LargestRealInverseCayley::realLambda(double er,
00493 double ei) const
00494 {
00495
00496 double reLambda = (sigma*(er*er+ei*ei) - (sigma+mu)*er + mu) /
00497 ( (er-1.0)*(er-1.0) + ei*ei);
00498 if (reLambda > sigma)
00499 return -1.0e6;
00500 else
00501 return reLambda;
00502 }