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 #ifndef ANASAZI_SOLVER_UTILS_HPP
00030 #define ANASAZI_SOLVER_UTILS_HPP
00031
00048 #include "AnasaziConfigDefs.hpp"
00049 #include "AnasaziMultiVecTraits.hpp"
00050 #include "AnasaziOperatorTraits.hpp"
00051 #include "Teuchos_ScalarTraits.hpp"
00052
00053 #include "AnasaziOutputManager.hpp"
00054 #include "Teuchos_BLAS.hpp"
00055 #include "Teuchos_LAPACK.hpp"
00056 #include "Teuchos_SerialDenseMatrix.hpp"
00057
00058 namespace Anasazi {
00059
00060 template<class ScalarType, class MV, class OP>
00061 class SolverUtils
00062 {
00063 public:
00064 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00065 typedef typename Teuchos::ScalarTraits<ScalarType> SCT;
00066
00068
00069
00071 SolverUtils();
00072
00074 virtual ~SolverUtils() {};
00075
00077
00079
00080
00082 static void permuteVectors(const int n, const std::vector<int> &perm, MV &Q, std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType >* resids = 0);
00083
00085 static void permuteVectors(const std::vector<int> &perm, Teuchos::SerialDenseMatrix<int,ScalarType> &Q);
00086
00088
00090
00091
00093
00114 static void applyHouse(int k, MV &V, const Teuchos::SerialDenseMatrix<int,ScalarType> &H, const std::vector<ScalarType> &tau, Teuchos::RCP<MV> workMV = Teuchos::null);
00115
00117
00119
00120
00122
00148 static int directSolver(int size, const Teuchos::SerialDenseMatrix<int,ScalarType> &KK,
00149 Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > MM,
00150 Teuchos::SerialDenseMatrix<int,ScalarType> &EV,
00151 std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &theta,
00152 int &nev, int esType = 0);
00154
00156
00157
00159
00161 static typename Teuchos::ScalarTraits<ScalarType>::magnitudeType errorEquality(const MV &X, const MV &MX, Teuchos::RCP<const OP> M = Teuchos::null);
00162
00164
00165 private:
00166
00168
00169
00170 typedef MultiVecTraits<ScalarType,MV> MVT;
00171 typedef OperatorTraits<ScalarType,MV,OP> OPT;
00172
00174 };
00175
00176
00177
00178
00179
00180
00181
00182 template<class ScalarType, class MV, class OP>
00183 SolverUtils<ScalarType, MV, OP>::SolverUtils() {}
00184
00185
00186
00187
00188
00189
00190
00191
00193
00194 template<class ScalarType, class MV, class OP>
00195 void SolverUtils<ScalarType, MV, OP>::permuteVectors(
00196 const int n,
00197 const std::vector<int> &perm,
00198 MV &Q,
00199 std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType >* resids)
00200 {
00201
00202
00203
00204 int i, j;
00205 std::vector<int> permcopy(perm), swapvec(n-1);
00206 std::vector<int> index(1);
00207 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00208 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00209
00210 TEST_FOR_EXCEPTION(n > MVT::GetNumberVecs(Q), std::invalid_argument, "Anasazi::SolverUtils::permuteVectors(): argument n larger than width of input multivector.");
00211
00212
00213
00214
00215
00216 for (i=0; i<n-1; i++) {
00217
00218
00219 for (j=i; j<n; j++) {
00220 if (permcopy[j] == i) {
00221
00222 break;
00223 }
00224 TEST_FOR_EXCEPTION(j == n-1, std::invalid_argument, "Anasazi::SolverUtils::permuteVectors(): permutation index invalid.");
00225 }
00226
00227
00228 std::swap( permcopy[j], permcopy[i] );
00229
00230 swapvec[i] = j;
00231 }
00232
00233
00234 for (i=n-2; i>=0; i--) {
00235 j = swapvec[i];
00236
00237
00238
00239
00240 if (resids) {
00241 std::swap( (*resids)[i], (*resids)[j] );
00242 }
00243
00244
00245 index[0] = j;
00246 Teuchos::RCP<MV> tmpQ = MVT::CloneCopy( Q, index );
00247 Teuchos::RCP<MV> tmpQj = MVT::CloneView( Q, index );
00248 index[0] = i;
00249 Teuchos::RCP<MV> tmpQi = MVT::CloneView( Q, index );
00250 MVT::MvAddMv( one, *tmpQi, zero, *tmpQi, *tmpQj );
00251 MVT::MvAddMv( one, *tmpQ, zero, *tmpQ, *tmpQi );
00252 }
00253 }
00254
00255
00257
00258 template<class ScalarType, class MV, class OP>
00259 void SolverUtils<ScalarType, MV, OP>::permuteVectors(
00260 const std::vector<int> &perm,
00261 Teuchos::SerialDenseMatrix<int,ScalarType> &Q)
00262 {
00263
00264
00265 Teuchos::BLAS<int,ScalarType> blas;
00266 const int n = perm.size();
00267 const int m = Q.numRows();
00268
00269 TEST_FOR_EXCEPTION(n != Q.numCols(), std::invalid_argument, "Anasazi::SolverUtils::permuteVectors(): size of permutation vector not equal to number of columns.");
00270
00271
00272 Teuchos::SerialDenseMatrix<int,ScalarType> copyQ( Q );
00273 for (int i=0; i<n; i++) {
00274 blas.COPY(m, copyQ[perm[i]], 1, Q[i], 1);
00275 }
00276 }
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286 template<class ScalarType, class MV, class OP>
00287 void SolverUtils<ScalarType, MV, OP>::applyHouse(int k, MV &V, const Teuchos::SerialDenseMatrix<int,ScalarType> &H, const std::vector<ScalarType> &tau, Teuchos::RCP<MV> workMV) {
00288
00289 const int n = MVT::GetNumberVecs(V);
00290 const ScalarType ONE = SCT::one();
00291 const ScalarType ZERO = SCT::zero();
00292
00293
00294 if (MVT::GetNumberVecs(V) == 0 || MVT::GetVecLength(V) == 0 || k == 0) {
00295 return;
00296 }
00297
00298 if (workMV == Teuchos::null) {
00299
00300 workMV = MVT::Clone(V,1);
00301 }
00302 else if (MVT::GetNumberVecs(*workMV) > 1) {
00303 std::vector<int> first(1);
00304 first[0] = 0;
00305 workMV = MVT::CloneView(*workMV,first);
00306 }
00307 else {
00308 TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*workMV) < 1,std::invalid_argument,"Anasazi::SolverUtils::applyHouse(): work multivector was empty.");
00309 }
00310
00311
00312 TEST_FOR_EXCEPTION( H.numCols() != k, std::invalid_argument,"Anasazi::SolverUtils::applyHouse(): H must have at least k columns.");
00313 TEST_FOR_EXCEPTION( (int)tau.size() != k, std::invalid_argument,"Anasazi::SolverUtils::applyHouse(): tau must have at least k entries.");
00314 TEST_FOR_EXCEPTION( H.numRows() != MVT::GetNumberVecs(V), std::invalid_argument,"Anasazi::SolverUtils::applyHouse(): Size of H,V are inconsistent.");
00315
00316
00317
00318 for (int i=0; i<k; i++) {
00319
00320
00321 std::vector<int> activeind(n-i);
00322 for (int j=0; j<n-i; j++) activeind[j] = j+i;
00323 Teuchos::RCP<MV> actV = MVT::CloneView(V,activeind);
00324
00325
00326
00327
00328
00329 Teuchos::SerialDenseMatrix<int,ScalarType> v(Teuchos::Copy,H,n-i,1,i,i);
00330
00331
00332 v(0,0) = ONE;
00333
00334
00335
00336
00337 MVT::MvTimesMatAddMv(-tau[i],*actV,v,ZERO,*workMV);
00338
00339
00340
00341 Teuchos::SerialDenseMatrix<int,ScalarType> vT(v,Teuchos::CONJ_TRANS);
00342 MVT::MvTimesMatAddMv(ONE,*workMV,vT,ONE,*actV);
00343
00344 actV = Teuchos::null;
00345 }
00346 }
00347
00348
00349
00350
00351
00352
00353
00354
00355 template<class ScalarType, class MV, class OP>
00356 int SolverUtils<ScalarType, MV, OP>::directSolver(
00357 int size,
00358 const Teuchos::SerialDenseMatrix<int,ScalarType> &KK,
00359 Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > MM,
00360 Teuchos::SerialDenseMatrix<int,ScalarType> &EV,
00361 std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &theta,
00362 int &nev, int esType)
00363 {
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406 Teuchos::LAPACK<int,ScalarType> lapack;
00407 Teuchos::BLAS<int,ScalarType> blas;
00408
00409 int rank = 0;
00410 int info = 0;
00411
00412 if (size < nev || size < 0) {
00413 return -1;
00414 }
00415 if (KK.numCols() < size || KK.numRows() < size) {
00416 return -2;
00417 }
00418 if ((esType == 0 || esType == 1)) {
00419 if (MM == Teuchos::null) {
00420 return -3;
00421 }
00422 else if (MM->numCols() < size || MM->numRows() < size) {
00423 return -3;
00424 }
00425 }
00426 if (EV.numCols() < size || EV.numRows() < size) {
00427 return -4;
00428 }
00429 if (theta.size() < (unsigned int) size) {
00430 return -5;
00431 }
00432 if (nev <= 0) {
00433 return -6;
00434 }
00435
00436
00437 std::string lapack_name = "hetrd";
00438 std::string lapack_opts = "u";
00439 int NB = lapack.ILAENV(1, lapack_name, lapack_opts, size, -1, -1, -1);
00440 int lwork = size*(NB+1);
00441 std::vector<ScalarType> work(lwork);
00442 std::vector<MagnitudeType> rwork(3*size-2);
00443
00444
00445 std::vector<MagnitudeType> tt( size );
00446 typedef typename std::vector<MagnitudeType>::iterator MTIter;
00447
00448 MagnitudeType tol = SCT::magnitude(SCT::squareroot(SCT::eps()));
00449
00450 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00451 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00452
00453 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > KKcopy, MMcopy;
00454 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > U;
00455
00456 switch (esType) {
00457 default:
00458 case 0:
00459
00460
00461
00462 for (rank = size; rank > 0; --rank) {
00463
00464 U = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(rank,rank) );
00465
00466
00467
00468 KKcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, KK, rank, rank ) );
00469 MMcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *MM, rank, rank ) );
00470
00471
00472
00473 info = 0;
00474 lapack.HEGV(1, 'V', 'U', rank, KKcopy->values(), KKcopy->stride(),
00475 MMcopy->values(), MMcopy->stride(), &tt[0], &work[0], lwork,
00476 &rwork[0], &info);
00477
00478
00479
00480 if (info < 0) {
00481 std::cerr << std::endl;
00482 std::cerr << "Anasazi::SolverUtils::directSolver(): In HEGV, argument " << -info << "has an illegal value.\n";
00483 std::cerr << std::endl;
00484 return -20;
00485 }
00486 if (info > 0) {
00487 if (info > rank)
00488 rank = info - rank;
00489 continue;
00490 }
00491
00492
00493
00494 MMcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *MM, rank, rank ) );
00495 for (int i = 0; i < rank; ++i) {
00496 for (int j = 0; j < i; ++j) {
00497 (*MMcopy)(i,j) = SCT::conjugate((*MM)(j,i));
00498 }
00499 }
00500
00501 TEST_FOR_EXCEPTION(
00502 U->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*MMcopy,*KKcopy,zero) != 0,
00503 std::logic_error, "Anasazi::SolverUtils::directSolver() call to Teuchos::SerialDenseMatrix::multiply() returned an error.");
00504
00505 TEST_FOR_EXCEPTION(
00506 MMcopy->multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,one,*KKcopy,*U,zero) != 0,
00507 std::logic_error, "Anasazi::SolverUtils::directSolver() call to Teuchos::SerialDenseMatrix::multiply() returned an error.");
00508 MagnitudeType maxNorm = SCT::magnitude(zero);
00509 MagnitudeType maxOrth = SCT::magnitude(zero);
00510 for (int i = 0; i < rank; ++i) {
00511 for (int j = i; j < rank; ++j) {
00512 if (j == i)
00513 maxNorm = SCT::magnitude((*MMcopy)(i,j) - one) > maxNorm
00514 ? SCT::magnitude((*MMcopy)(i,j) - one) : maxNorm;
00515 else
00516 maxOrth = SCT::magnitude((*MMcopy)(i,j)) > maxOrth
00517 ? SCT::magnitude((*MMcopy)(i,j)) : maxOrth;
00518 }
00519 }
00520
00521
00522
00523
00524
00525
00526
00527
00528 if ((maxNorm <= tol) && (maxOrth <= tol)) {
00529 break;
00530 }
00531 }
00532
00533
00534
00535
00536
00537 nev = (rank < nev) ? rank : nev;
00538 EV.putScalar( zero );
00539 std::copy(tt.begin(),tt.begin()+nev,theta.begin());
00540 for (int i = 0; i < nev; ++i) {
00541 blas.COPY( rank, (*KKcopy)[i], 1, EV[i], 1 );
00542 }
00543 break;
00544
00545 case 1:
00546
00547
00548
00549
00550
00551 KKcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, KK, size, size ) );
00552 MMcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *MM, size, size ) );
00553
00554
00555
00556 info = 0;
00557 lapack.HEGV(1, 'V', 'U', size, KKcopy->values(), KKcopy->stride(),
00558 MMcopy->values(), MMcopy->stride(), &tt[0], &work[0], lwork,
00559 &rwork[0], &info);
00560
00561
00562
00563 if (info < 0) {
00564 std::cerr << std::endl;
00565 std::cerr << "Anasazi::SolverUtils::directSolver(): In HEGV, argument " << -info << "has an illegal value.\n";
00566 std::cerr << std::endl;
00567 return -20;
00568 }
00569 if (info > 0) {
00570 if (info > size)
00571 nev = 0;
00572 else {
00573 std::cerr << std::endl;
00574 std::cerr << "Anasazi::SolverUtils::directSolver(): In HEGV, DPOTRF or DHEEV returned an error code (" << info << ").\n";
00575 std::cerr << std::endl;
00576 return -20;
00577 }
00578 }
00579
00580
00581
00582 std::copy(tt.begin(),tt.begin()+nev,theta.begin());
00583 for (int i = 0; i < nev; ++i) {
00584 blas.COPY( size, (*KKcopy)[i], 1, EV[i], 1 );
00585 }
00586 break;
00587
00588 case 10:
00589
00590
00591
00592
00593
00594 KKcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, KK, size, size ) );
00595
00596
00597
00598 lapack.HEEV('V', 'U', size, KKcopy->values(), KKcopy->stride(), &tt[0], &work[0], lwork, &rwork[0], &info);
00599
00600
00601 if (info != 0) {
00602 std::cerr << std::endl;
00603 if (info < 0) {
00604 std::cerr << "Anasazi::SolverUtils::directSolver(): In DHEEV, argument " << -info << " has an illegal value\n";
00605 }
00606 else {
00607 std::cerr << "Anasazi::SolverUtils::directSolver(): In DHEEV, the algorithm failed to converge (" << info << ").\n";
00608 }
00609 std::cerr << std::endl;
00610 info = -20;
00611 break;
00612 }
00613
00614
00615
00616 std::copy(tt.begin(),tt.begin()+nev,theta.begin());
00617 for (int i = 0; i < nev; ++i) {
00618 blas.COPY( size, (*KKcopy)[i], 1, EV[i], 1 );
00619 }
00620 break;
00621 }
00622
00623 return info;
00624 }
00625
00626
00627
00628
00629
00630
00631
00632
00633 template<class ScalarType, class MV, class OP>
00634 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
00635 SolverUtils<ScalarType, MV, OP>::errorEquality(const MV &X, const MV &MX, Teuchos::RCP<const OP> M)
00636 {
00637
00638
00639
00640
00641 MagnitudeType maxDiff = SCT::magnitude(SCT::zero());
00642
00643 int xc = MVT::GetNumberVecs(X);
00644 int mxc = MVT::GetNumberVecs(MX);
00645
00646 TEST_FOR_EXCEPTION(xc != mxc,std::invalid_argument,"Anasazi::SolverUtils::errorEquality(): input multivecs have different number of columns.");
00647 if (xc == 0) {
00648 return maxDiff;
00649 }
00650
00651 MagnitudeType maxCoeffX = SCT::magnitude(SCT::zero());
00652 std::vector<MagnitudeType> tmp( xc );
00653 MVT::MvNorm(MX, tmp);
00654
00655 for (int i = 0; i < xc; ++i) {
00656 maxCoeffX = (tmp[i] > maxCoeffX) ? tmp[i] : maxCoeffX;
00657 }
00658
00659 std::vector<int> index( 1 );
00660 Teuchos::RCP<MV> MtimesX;
00661 if (M != Teuchos::null) {
00662 MtimesX = MVT::Clone( X, xc );
00663 OPT::Apply( *M, X, *MtimesX );
00664 }
00665 else {
00666 MtimesX = MVT::CloneCopy(X);
00667 }
00668 MVT::MvAddMv( -1.0, MX, 1.0, *MtimesX, *MtimesX );
00669 MVT::MvNorm( *MtimesX, tmp );
00670
00671 for (int i = 0; i < xc; ++i) {
00672 maxDiff = (tmp[i] > maxDiff) ? tmp[i] : maxDiff;
00673 }
00674
00675 return (maxCoeffX == 0.0) ? maxDiff : maxDiff/maxCoeffX;
00676
00677 }
00678
00679 }
00680
00681 #endif // ANASAZI_SOLVER_UTILS_HPP
00682