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_BLOCK_KRYLOV_SCHUR_SOLMGR_HPP
00030 #define ANASAZI_BLOCK_KRYLOV_SCHUR_SOLMGR_HPP
00031
00036 #include "AnasaziConfigDefs.hpp"
00037 #include "AnasaziTypes.hpp"
00038
00039 #include "AnasaziEigenproblem.hpp"
00040 #include "AnasaziSolverManager.hpp"
00041 #include "AnasaziSolverUtils.hpp"
00042
00043 #include "AnasaziBlockKrylovSchur.hpp"
00044 #include "AnasaziBasicSort.hpp"
00045 #include "AnasaziSVQBOrthoManager.hpp"
00046 #include "AnasaziBasicOrthoManager.hpp"
00047 #include "AnasaziStatusTestResNorm.hpp"
00048 #include "AnasaziStatusTestWithOrdering.hpp"
00049 #include "AnasaziStatusTestCombo.hpp"
00050 #include "AnasaziStatusTestOutput.hpp"
00051 #include "AnasaziBasicOutputManager.hpp"
00052 #include "Teuchos_BLAS.hpp"
00053 #include "Teuchos_LAPACK.hpp"
00054 #include "Teuchos_TimeMonitor.hpp"
00055
00056
00077 namespace Anasazi {
00078
00079
00106 template<class ScalarType, class MV, class OP>
00107 class BlockKrylovSchurSolMgr : public SolverManager<ScalarType,MV,OP> {
00108
00109 private:
00110 typedef MultiVecTraits<ScalarType,MV> MVT;
00111 typedef OperatorTraits<ScalarType,MV,OP> OPT;
00112 typedef Teuchos::ScalarTraits<ScalarType> SCT;
00113 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00114 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00115
00116 public:
00117
00119
00120
00138 BlockKrylovSchurSolMgr( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00139 Teuchos::ParameterList &pl );
00140
00142 virtual ~BlockKrylovSchurSolMgr() {};
00144
00146
00147
00149 const Eigenproblem<ScalarType,MV,OP>& getProblem() const {
00150 return *_problem;
00151 }
00152
00154 int getNumIters() const {
00155 return _numIters;
00156 }
00157
00160 std::vector<Value<ScalarType> > getRitzValues() const {
00161 std::vector<Value<ScalarType> > ret( _ritzValues );
00162 return ret;
00163 }
00164
00171 Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00172 return tuple(_timerSolve, _timerRestarting);
00173 }
00174
00176
00178
00179
00198 ReturnType solve();
00199
00201 void setGlobalStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global);
00202
00204 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getGlobalStatusTest() const;
00205
00207 void setDebugStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug);
00208
00210 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getDebugStatusTest() const;
00211
00213
00214 private:
00215 Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > _problem;
00216 Teuchos::RCP<SortManager<MagnitudeType> > _sort;
00217
00218 std::string _whch, _ortho;
00219 MagnitudeType _ortho_kappa;
00220
00221 MagnitudeType _convtol;
00222 int _maxRestarts;
00223 bool _relconvtol,_conjSplit;
00224 int _blockSize, _numBlocks, _stepSize, _nevBlocks, _xtra_nevBlocks;
00225 int _numIters;
00226 int _verbosity;
00227 bool _inSituRestart;
00228
00229 std::vector<Value<ScalarType> > _ritzValues;
00230
00231 Teuchos::RCP<Teuchos::Time> _timerSolve, _timerRestarting;
00232
00233 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > globalTest_;
00234 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > debugTest_;
00235
00236 int _printNum;
00237 };
00238
00239
00240
00241 template<class ScalarType, class MV, class OP>
00242 BlockKrylovSchurSolMgr<ScalarType,MV,OP>::BlockKrylovSchurSolMgr(
00243 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00244 Teuchos::ParameterList &pl ) :
00245 _problem(problem),
00246 _whch("LM"),
00247 _ortho("SVQB"),
00248 _ortho_kappa(-1.0),
00249 _convtol(0),
00250 _maxRestarts(20),
00251 _relconvtol(true),
00252 _conjSplit(false),
00253 _blockSize(0),
00254 _numBlocks(0),
00255 _stepSize(0),
00256 _nevBlocks(0),
00257 _xtra_nevBlocks(0),
00258 _numIters(0),
00259 _verbosity(Anasazi::Errors),
00260 _inSituRestart(false),
00261 _timerSolve(Teuchos::TimeMonitor::getNewTimer("BlockKrylovSchurSolMgr::solve()")),
00262 _timerRestarting(Teuchos::TimeMonitor::getNewTimer("BlockKrylovSchurSolMgr restarting")),
00263 _printNum(-1)
00264 {
00265 TEST_FOR_EXCEPTION(_problem == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
00266 TEST_FOR_EXCEPTION(!_problem->isProblemSet(), std::invalid_argument, "Problem not set.");
00267 TEST_FOR_EXCEPTION(_problem->getInitVec() == Teuchos::null, std::invalid_argument, "Problem does not contain initial vectors to clone from.");
00268
00269 const int nev = _problem->getNEV();
00270
00271
00272 _convtol = pl.get("Convergence Tolerance",MT::prec());
00273 _relconvtol = pl.get("Relative Convergence Tolerance",_relconvtol);
00274
00275
00276 _maxRestarts = pl.get("Maximum Restarts",_maxRestarts);
00277
00278
00279 _blockSize = pl.get("Block Size",1);
00280 TEST_FOR_EXCEPTION(_blockSize <= 0, std::invalid_argument,
00281 "Anasazi::BlockKrylovSchurSolMgr: \"Block Size\" must be strictly positive.");
00282
00283
00284 _xtra_nevBlocks = pl.get("Extra NEV Blocks",0);
00285 if (nev%_blockSize) {
00286 _nevBlocks = nev/_blockSize + _xtra_nevBlocks + 1;
00287 } else {
00288 _nevBlocks = nev/_blockSize + _xtra_nevBlocks;
00289 }
00290
00291 _numBlocks = pl.get("Num Blocks",3*_nevBlocks);
00292 TEST_FOR_EXCEPTION(_numBlocks <= _nevBlocks, std::invalid_argument,
00293 "Anasazi::BlockKrylovSchurSolMgr: \"Num Blocks\" must be strictly positive and large enough to compute the requested eigenvalues.");
00294
00295 TEST_FOR_EXCEPTION(_numBlocks*_blockSize > MVT::GetVecLength(*_problem->getInitVec()),
00296 std::invalid_argument,
00297 "Anasazi::BlockKrylovSchurSolMgr: Potentially impossible orthogonality requests. Reduce basis size.");
00298
00299
00300 if (_maxRestarts) {
00301 _stepSize = pl.get("Step Size", (_maxRestarts+1)*(_numBlocks+1));
00302 } else {
00303 _stepSize = pl.get("Step Size", _numBlocks+1);
00304 }
00305 TEST_FOR_EXCEPTION(_stepSize < 1, std::invalid_argument,
00306 "Anasazi::BlockKrylovSchurSolMgr: \"Step Size\" must be strictly positive.");
00307
00308
00309 if (pl.isParameter("Sort Manager")) {
00310 _sort = Teuchos::getParameter<Teuchos::RCP<Anasazi::SortManager<MagnitudeType> > >(pl,"Sort Manager");
00311 } else {
00312
00313 _whch = pl.get("Which",_whch);
00314 TEST_FOR_EXCEPTION(_whch != "SM" && _whch != "LM" && _whch != "SR" && _whch != "LR" && _whch != "SI" && _whch != "LI",
00315 std::invalid_argument, "Invalid sorting string.");
00316 _sort = Teuchos::rcp( new BasicSort<MagnitudeType>(_whch) );
00317 }
00318
00319
00320 _ortho = pl.get("Orthogonalization",_ortho);
00321 if (_ortho != "DGKS" && _ortho != "SVQB") {
00322 _ortho = "SVQB";
00323 }
00324
00325
00326 _ortho_kappa = pl.get("Orthogonalization Constant",_ortho_kappa);
00327
00328
00329 if (pl.isParameter("Verbosity")) {
00330 if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
00331 _verbosity = pl.get("Verbosity", _verbosity);
00332 } else {
00333 _verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
00334 }
00335 }
00336
00337
00338 if (pl.isParameter("In Situ Restarting")) {
00339 if (Teuchos::isParameterType<bool>(pl,"In Situ Restarting")) {
00340 _inSituRestart = pl.get("In Situ Restarting",_inSituRestart);
00341 } else {
00342 _inSituRestart = (bool)Teuchos::getParameter<int>(pl,"In Situ Restarting");
00343 }
00344 }
00345
00346 _printNum = pl.get<int>("Print Number of Ritz Values",-1);
00347 }
00348
00349
00350
00351 template<class ScalarType, class MV, class OP>
00352 ReturnType
00353 BlockKrylovSchurSolMgr<ScalarType,MV,OP>::solve() {
00354
00355 const int nev = _problem->getNEV();
00356 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00357 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00358
00359 Teuchos::BLAS<int,ScalarType> blas;
00360 Teuchos::LAPACK<int,ScalarType> lapack;
00361 typedef SolverUtils<ScalarType,MV,OP> msutils;
00362
00364
00365 Teuchos::RCP<BasicOutputManager<ScalarType> > printer = Teuchos::rcp( new BasicOutputManager<ScalarType>(_verbosity) );
00366
00368
00369
00370
00371 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convtest;
00372 if (globalTest_ == Teuchos::null) {
00373 convtest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(_convtol,nev,StatusTestResNorm<ScalarType,MV,OP>::RITZRES_2NORM,_relconvtol) );
00374 }
00375 else {
00376 convtest = globalTest_;
00377 }
00378 Teuchos::RCP<StatusTestWithOrdering<ScalarType,MV,OP> > ordertest
00379 = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,_sort,nev) );
00380
00381 Teuchos::Array<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
00382 alltests.push_back(ordertest);
00383
00384 if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
00385
00386 Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > combotest
00387 = Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP>( StatusTestCombo<ScalarType,MV,OP>::OR, alltests) );
00388
00389 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
00390 if ( printer->isVerbosity(Debug) ) {
00391 outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed+Failed+Undefined ) );
00392 }
00393 else {
00394 outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed ) );
00395 }
00396
00398
00399 Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho;
00400 if (_ortho=="SVQB") {
00401 ortho = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(_problem->getM()) );
00402 } else if (_ortho=="DGKS") {
00403 if (_ortho_kappa <= 0) {
00404 ortho = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>(_problem->getM()) );
00405 }
00406 else {
00407 ortho = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>(_problem->getM(),_ortho_kappa) );
00408 }
00409 } else {
00410 TEST_FOR_EXCEPTION(_ortho!="SVQB"&&_ortho!="DGKS",std::logic_error,"Anasazi::BlockKrylovSchurSolMgr::solve(): Invalid orthogonalization type.");
00411 }
00412
00414
00415 Teuchos::ParameterList plist;
00416 plist.set("Block Size",_blockSize);
00417 plist.set("Num Blocks",_numBlocks);
00418 plist.set("Step Size",_stepSize);
00419 if (_printNum == -1) {
00420 plist.set("Print Number of Ritz Values",_nevBlocks*_blockSize);
00421 }
00422 else {
00423 plist.set("Print Number of Ritz Values",_printNum);
00424 }
00425
00427
00428 Teuchos::RCP<BlockKrylovSchur<ScalarType,MV,OP> > bks_solver
00429 = Teuchos::rcp( new BlockKrylovSchur<ScalarType,MV,OP>(_problem,_sort,printer,outputtest,ortho,plist) );
00430
00431 Teuchos::RCP< const MV > probauxvecs = _problem->getAuxVecs();
00432 if (probauxvecs != Teuchos::null) {
00433 bks_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
00434 }
00435
00436
00437
00438
00439
00440
00441 Teuchos::RCP<MV> workMV;
00442 if (_maxRestarts > 0) {
00443 if (_inSituRestart==true) {
00444
00445 workMV = MVT::Clone( *_problem->getInitVec(), 1 );
00446 }
00447 else {
00448 if (_problem->isHermitian()) {
00449 workMV = MVT::Clone( *_problem->getInitVec(), _nevBlocks*_blockSize + _blockSize );
00450 } else {
00451 workMV = MVT::Clone( *_problem->getInitVec(), _nevBlocks*_blockSize+1 + _blockSize );
00452 }
00453 }
00454 } else {
00455 workMV = Teuchos::null;
00456 }
00457
00458
00459 Eigensolution<ScalarType,MV> sol;
00460 sol.numVecs = 0;
00461 _problem->setSolution(sol);
00462
00463 int numRestarts = 0;
00464 int cur_nevBlocks = 0;
00465
00466
00467 {
00468 Teuchos::TimeMonitor slvtimer(*_timerSolve);
00469
00470
00471 while (1) {
00472 try {
00473 bks_solver->iterate();
00474
00476
00477
00478
00480 if (ordertest->getStatus() == Passed ) {
00481
00482
00483
00484 break;
00485 }
00487
00488
00489
00491
00492
00493
00494
00495
00496
00497
00498 else if ( (bks_solver->getCurSubspaceDim() == bks_solver->getMaxSubspaceDim()) ||
00499 (!_problem->isHermitian() && !_conjSplit && (bks_solver->getCurSubspaceDim()+1 == bks_solver->getMaxSubspaceDim())) ) {
00500
00501 Teuchos::TimeMonitor restimer(*_timerRestarting);
00502
00503 if ( numRestarts >= _maxRestarts ) {
00504 break;
00505 }
00506 numRestarts++;
00507
00508 printer->stream(Debug) << " Performing restart number " << numRestarts << " of " << _maxRestarts << std::endl << std::endl;
00509
00510
00511 if (!bks_solver->isSchurCurrent())
00512 bks_solver->computeSchurForm( true );
00513
00514
00515 _ritzValues = bks_solver->getRitzValues();
00516
00517
00518 BlockKrylovSchurState<ScalarType,MV> oldState = bks_solver->getState();
00519
00520
00521 int curDim = oldState.curDim;
00522
00523
00524 std::vector<int> ritzIndex = bks_solver->getRitzIndex();
00525 if (ritzIndex[_nevBlocks*_blockSize-1]==1) {
00526 _conjSplit = true;
00527 cur_nevBlocks = _nevBlocks*_blockSize+1;
00528 } else {
00529 _conjSplit = false;
00530 cur_nevBlocks = _nevBlocks*_blockSize;
00531 }
00532
00533
00534
00535
00536 Teuchos::SerialDenseMatrix<int,ScalarType> Qnev(Teuchos::View, *(oldState.Q), curDim, cur_nevBlocks);
00537
00538
00539 std::vector<int> curind( curDim );
00540 for (int i=0; i<curDim; i++) { curind[i] = i; }
00541 Teuchos::RCP<const MV> basistemp = MVT::CloneView( *(oldState.V), curind );
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551 Teuchos::RCP<MV> newF;
00552 if (_inSituRestart) {
00553
00554
00555 Teuchos::RCP<MV> solverbasis = Teuchos::rcp_const_cast<MV>(oldState.V);
00556 Teuchos::SerialDenseMatrix<int,ScalarType> copyQnev(Qnev);
00557
00558
00559 std::vector<ScalarType> tau(cur_nevBlocks), work(cur_nevBlocks);
00560 int info;
00561 lapack.GEQRF(curDim,cur_nevBlocks,copyQnev.values(),copyQnev.stride(),&tau[0],&work[0],work.size(),&info);
00562 TEST_FOR_EXCEPTION(info != 0,std::logic_error,
00563 "Anasazi::BlockDavidsonSolMgr::solve(): error calling GEQRF during restarting.");
00564
00565 std::vector<ScalarType> d(cur_nevBlocks);
00566 for (int j=0; j<copyQnev.numCols(); j++) {
00567 d[j] = copyQnev(j,j);
00568 }
00569 if (printer->isVerbosity(Debug)) {
00570 Teuchos::SerialDenseMatrix<int,ScalarType> R(Teuchos::Copy,copyQnev,cur_nevBlocks,cur_nevBlocks);
00571 for (int j=0; j<R.numCols(); j++) {
00572 R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
00573 for (int i=j+1; i<R.numRows(); i++) {
00574 R(i,j) = zero;
00575 }
00576 }
00577 printer->stream(Debug) << "||Triangular factor of Su - I||: " << R.normFrobenius() << std::endl;
00578 }
00579
00580
00581
00582
00583 curind.resize(curDim);
00584 for (int i=0; i<curDim; i++) curind[i] = i;
00585 Teuchos::RCP<MV> oldV = MVT::CloneView(*solverbasis,curind);
00586 msutils::applyHouse(cur_nevBlocks,*oldV,copyQnev,tau,workMV);
00587
00588 oldV = Teuchos::null;
00589
00590
00591 curind.resize(cur_nevBlocks);
00592 for (int i=0; i<cur_nevBlocks; i++) { curind[i] = i; }
00593 oldV = MVT::CloneView( *solverbasis, curind );
00594 MVT::MvScale(*oldV,d);
00595 oldV = Teuchos::null;
00596
00597 curind.resize(_blockSize);
00598 for (int i=0; i<_blockSize; i++) { curind[i] = cur_nevBlocks + i; }
00599 newF = MVT::CloneView( *solverbasis, curind );
00600 }
00601 else {
00602
00603 curind.resize(cur_nevBlocks);
00604 for (int i=0; i<cur_nevBlocks; i++) { curind[i] = i; }
00605 Teuchos::RCP<MV> tmp_newV = MVT::CloneView(*workMV, curind );
00606
00607 MVT::MvTimesMatAddMv( one, *basistemp, Qnev, zero, *tmp_newV );
00608 tmp_newV = Teuchos::null;
00609
00610 curind.resize(_blockSize);
00611 for (int i=0; i<_blockSize; i++) { curind[i] = cur_nevBlocks + i; }
00612 newF = MVT::CloneView( *workMV, curind );
00613 }
00614
00615
00616 curind.resize(_blockSize);
00617 for (int i=0; i<_blockSize; i++) { curind[i] = curDim + i; }
00618 Teuchos::RCP<const MV> oldF = MVT::CloneView( *(oldState.V), curind );
00619 for (int i=0; i<_blockSize; i++) { curind[i] = i; }
00620 MVT::SetBlock( *oldF, curind, *newF );
00621 newF = Teuchos::null;
00622
00623
00624
00625
00626
00627 Teuchos::SerialDenseMatrix<int,ScalarType> oldS(Teuchos::View, *(oldState.S), cur_nevBlocks+_blockSize, cur_nevBlocks);
00628 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > newH =
00629 Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( oldS ) );
00630
00631
00632 Teuchos::SerialDenseMatrix<int,ScalarType> oldB(Teuchos::View, *(oldState.H), _blockSize, _blockSize, curDim, curDim-_blockSize);
00633
00634
00635 Teuchos::SerialDenseMatrix<int,ScalarType> subQ(Teuchos::View, *(oldState.Q), _blockSize, cur_nevBlocks, curDim-_blockSize);
00636
00637
00638 Teuchos::SerialDenseMatrix<int,ScalarType> newB(Teuchos::View, *newH, _blockSize, cur_nevBlocks, cur_nevBlocks);
00639
00640
00641 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, _blockSize, cur_nevBlocks, _blockSize, one,
00642 oldB.values(), oldB.stride(), subQ.values(), subQ.stride(), zero, newB.values(), newB.stride() );
00643
00644
00645
00646
00647 BlockKrylovSchurState<ScalarType,MV> newstate;
00648 if (_inSituRestart) {
00649 newstate.V = oldState.V;
00650 } else {
00651 newstate.V = workMV;
00652 }
00653 newstate.H = newH;
00654 newstate.curDim = cur_nevBlocks;
00655 bks_solver->initialize(newstate);
00656
00657 }
00659
00660
00661
00662
00664 else {
00665 TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::BlockKrylovSchurSolMgr::solve(): Invalid return from bks_solver::iterate().");
00666 }
00667 }
00668 catch (AnasaziError err) {
00669 printer->stream(Errors)
00670 << "Anasazi::BlockKrylovSchurSolMgr::solve() caught unexpected exception from Anasazi::BlockKrylovSchur::iterate() at iteration " << bks_solver->getNumIters() << std::endl
00671 << err.what() << std::endl
00672 << "Anasazi::BlockKrylovSchurSolMgr::solve() returning Unconverged with no solutions." << std::endl;
00673 return Unconverged;
00674 }
00675 }
00676
00677
00678
00679 workMV = Teuchos::null;
00680
00681
00682 _ritzValues = bks_solver->getRitzValues();
00683
00684 sol.numVecs = ordertest->howMany();
00685 std::vector<int> whichVecs = ordertest->whichVecs();
00686
00687
00688 if (sol.numVecs > 0) {
00689
00690 bool keepMore = false;
00691 int numEvecs = sol.numVecs;
00692 if (whichVecs[sol.numVecs-1] > (sol.numVecs-1)) {
00693 keepMore = true;
00694 numEvecs = whichVecs[sol.numVecs-1]+1;
00695 }
00696
00697
00698 std::vector<int> tmpIndex = bks_solver->getRitzIndex();
00699 if (tmpIndex[whichVecs[sol.numVecs-1]]==1) {
00700 sol.numVecs++;
00701 numEvecs++;
00702 }
00703
00704
00705 bks_solver->setNumRitzVectors(numEvecs);
00706 bks_solver->computeRitzVectors();
00707
00708
00709
00710
00711 if (!keepMore) {
00712 sol.index = bks_solver->getRitzIndex();
00713 sol.Evals = bks_solver->getRitzValues();
00714 sol.Evecs = MVT::CloneCopy( *(bks_solver->getRitzVectors()) );
00715 }
00716
00717
00718
00719 sol.Evals.resize(sol.numVecs);
00720 sol.index.resize(sol.numVecs);
00721
00722
00723 if (keepMore) {
00724 std::vector<Anasazi::Value<ScalarType> > tmpEvals = bks_solver->getRitzValues();
00725 for (int vec_i=0; vec_i<sol.numVecs; ++vec_i) {
00726 sol.index[vec_i] = tmpIndex[whichVecs[vec_i]];
00727 sol.Evals[vec_i] = tmpEvals[whichVecs[vec_i]];
00728 }
00729 sol.Evecs = MVT::CloneCopy( *(bks_solver->getRitzVectors()), whichVecs );
00730 }
00731
00732
00733 sol.Espace = sol.Evecs;
00734 }
00735 }
00736
00737
00738 bks_solver->currentStatus(printer->stream(FinalSummary));
00739
00740
00741 Teuchos::TimeMonitor::summarize(printer->stream(TimingDetails));
00742
00743 _problem->setSolution(sol);
00744 printer->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
00745
00746
00747 _numIters = bks_solver->getNumIters();
00748
00749 if (sol.numVecs < nev) {
00750 return Unconverged;
00751 }
00752 return Converged;
00753 }
00754
00755
00756 template <class ScalarType, class MV, class OP>
00757 void
00758 BlockKrylovSchurSolMgr<ScalarType,MV,OP>::setGlobalStatusTest(
00759 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global)
00760 {
00761 globalTest_ = global;
00762 }
00763
00764 template <class ScalarType, class MV, class OP>
00765 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
00766 BlockKrylovSchurSolMgr<ScalarType,MV,OP>::getGlobalStatusTest() const
00767 {
00768 return globalTest_;
00769 }
00770
00771 template <class ScalarType, class MV, class OP>
00772 void
00773 BlockKrylovSchurSolMgr<ScalarType,MV,OP>::setDebugStatusTest(
00774 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug)
00775 {
00776 debugTest_ = debug;
00777 }
00778
00779 template <class ScalarType, class MV, class OP>
00780 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
00781 BlockKrylovSchurSolMgr<ScalarType,MV,OP>::getDebugStatusTest() const
00782 {
00783 return debugTest_;
00784 }
00785
00786 }
00787
00788 #endif