299 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
300 Teuchos::ParameterList &pl ) :
316 inSituRestart_(false),
320#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
321 ,timerSolve_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockKrylovSchurSolMgr::solve()")),
322 timerRestarting_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockKrylovSchurSolMgr restarting"))
325 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument,
"Problem not given to solver manager.");
326 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument,
"Problem not set.");
327 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null, std::invalid_argument,
"Problem does not contain initial vectors to clone from.");
329 const int nev = problem_->getNEV();
332 convtol_ = pl.get(
"Convergence Tolerance",MT::prec());
333 relconvtol_ = pl.get(
"Relative Convergence Tolerance",relconvtol_);
336 maxRestarts_ = pl.get(
"Maximum Restarts",maxRestarts_);
339 blockSize_ = pl.get(
"Block Size",1);
340 TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
341 "Anasazi::BlockKrylovSchurSolMgr: \"Block Size\" must be strictly positive.");
344 xtra_nevBlocks_ = pl.get(
"Extra NEV Blocks",0);
345 if (nev%blockSize_) {
346 nevBlocks_ = nev/blockSize_ + 1;
348 nevBlocks_ = nev/blockSize_;
354 if (pl.isParameter(
"Dynamic Extra NEV")) {
355 if (Teuchos::isParameterType<bool>(pl,
"Dynamic Extra NEV")) {
356 dynXtraNev_ = pl.get(
"Dynamic Extra NEV",dynXtraNev_);
358 dynXtraNev_ = ( Teuchos::getParameter<int>(pl,
"Dynamic Extra NEV") != 0 );
362 numBlocks_ = pl.get(
"Num Blocks",3*nevBlocks_);
363 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= nevBlocks_, std::invalid_argument,
364 "Anasazi::BlockKrylovSchurSolMgr: \"Num Blocks\" must be strictly positive and large enough to compute the requested eigenvalues.");
366 TEUCHOS_TEST_FOR_EXCEPTION(
static_cast<ptrdiff_t
>(numBlocks_)*blockSize_ > MVT::GetGlobalLength(*problem_->getInitVec()),
367 std::invalid_argument,
368 "Anasazi::BlockKrylovSchurSolMgr: Potentially impossible orthogonality requests. Reduce basis size.");
372 stepSize_ = pl.get(
"Step Size", (maxRestarts_+1)*(numBlocks_+1));
374 stepSize_ = pl.get(
"Step Size", numBlocks_+1);
376 TEUCHOS_TEST_FOR_EXCEPTION(stepSize_ < 1, std::invalid_argument,
377 "Anasazi::BlockKrylovSchurSolMgr: \"Step Size\" must be strictly positive.");
380 if (pl.isParameter(
"Sort Manager")) {
381 sort_ = Teuchos::getParameter<Teuchos::RCP<Anasazi::SortManager<MagnitudeType> > >(pl,
"Sort Manager");
384 whch_ = pl.get(
"Which",whch_);
385 TEUCHOS_TEST_FOR_EXCEPTION(whch_ !=
"SM" && whch_ !=
"LM" && whch_ !=
"SR" && whch_ !=
"LR" && whch_ !=
"SI" && whch_ !=
"LI",
386 std::invalid_argument,
"Invalid sorting string.");
387 sort_ = Teuchos::rcp(
new BasicSort<MagnitudeType>(whch_) );
391 ortho_ = pl.get(
"Orthogonalization",ortho_);
392 if (ortho_ !=
"DGKS" && ortho_ !=
"SVQB" && ortho_ !=
"ICGS") {
397 ortho_kappa_ = pl.get(
"Orthogonalization Constant",ortho_kappa_);
401 osProc_ = pl.get(
"Output Processor", osProc_);
404 if (pl.isParameter(
"Output Stream")) {
405 osp_ = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,
"Output Stream");
412 if (pl.isParameter(
"Verbosity")) {
413 if (Teuchos::isParameterType<int>(pl,
"Verbosity")) {
414 verbosity_ = pl.get(
"Verbosity", verbosity_);
416 verbosity_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,
"Verbosity");
421 if (pl.isParameter(
"In Situ Restarting")) {
422 if (Teuchos::isParameterType<bool>(pl,
"In Situ Restarting")) {
423 inSituRestart_ = pl.get(
"In Situ Restarting",inSituRestart_);
425 inSituRestart_ = ( Teuchos::getParameter<int>(pl,
"In Situ Restarting") != 0 );
429 printNum_ = pl.get<
int>(
"Print Number of Ritz Values",printNum_);
438 const int nev = problem_->getNEV();
439 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
440 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
442 Teuchos::BLAS<int,ScalarType> blas;
443 Teuchos::LAPACK<int,ScalarType> lapack;
447 Teuchos::RCP<OutputManager<ScalarType> > printer = Teuchos::rcp(
new OutputManager<ScalarType>(verbosity_,osp_) );
453 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convtest;
454 if (globalTest_ == Teuchos::null) {
455 convtest = Teuchos::rcp(
new StatusTestResNorm<ScalarType,MV,OP>(convtol_,nev,RITZRES_2NORM,relconvtol_) );
458 convtest = globalTest_;
460 Teuchos::RCP<StatusTestWithOrdering<ScalarType,MV,OP> > ordertest
461 = Teuchos::rcp(
new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sort_,nev) );
463 Teuchos::Array<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
464 alltests.push_back(ordertest);
466 if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
468 Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > combotest
471 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
472 if ( printer->isVerbosity(
Debug) ) {
473 outputtest = Teuchos::rcp(
new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,
Passed+
Failed+
Undefined ) );
476 outputtest = Teuchos::rcp(
new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,
Passed ) );
481 Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho;
482 if (ortho_==
"SVQB") {
483 ortho = Teuchos::rcp(
new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
484 }
else if (ortho_==
"DGKS") {
485 if (ortho_kappa_ <= 0) {
486 ortho = Teuchos::rcp(
new BasicOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
488 ortho = Teuchos::rcp(
new BasicOrthoManager<ScalarType,MV,OP>(problem_->getM(),ortho_kappa_) );
490 }
else if (ortho_==
"ICGS") {
491 ortho = Teuchos::rcp(
new ICGSOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
493 TEUCHOS_TEST_FOR_EXCEPTION(ortho_!=
"SVQB"&&ortho_!=
"DGKS"&&ortho_!=
"ICGS",std::logic_error,
"Anasazi::BlockKrylovSchurSolMgr::solve(): Invalid orthogonalization type.");
498 Teuchos::ParameterList plist;
499 plist.set(
"Block Size",blockSize_);
500 plist.set(
"Num Blocks",numBlocks_);
501 plist.set(
"Step Size",stepSize_);
502 if (printNum_ == -1) {
503 plist.set(
"Print Number of Ritz Values",nevBlocks_*blockSize_);
506 plist.set(
"Print Number of Ritz Values",printNum_);
511 Teuchos::RCP<BlockKrylovSchur<ScalarType,MV,OP> > bks_solver
512 = Teuchos::rcp(
new BlockKrylovSchur<ScalarType,MV,OP>(problem_,sort_,printer,outputtest,ortho,plist) );
514 Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
515 if (probauxvecs != Teuchos::null) {
516 bks_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
526 int maxXtraBlocks = 0;
527 if ( dynXtraNev_ ) maxXtraBlocks = ( ( bks_solver->getMaxSubspaceDim() - nev ) / blockSize_ ) / 2;
529 Teuchos::RCP<MV> workMV;
530 if (maxRestarts_ > 0) {
531 if (inSituRestart_==
true) {
533 workMV = MVT::Clone( *problem_->getInitVec(), 1 );
536 if (problem_->isHermitian()) {
537 workMV = MVT::Clone( *problem_->getInitVec(), (nevBlocks_+maxXtraBlocks)*blockSize_ + blockSize_ );
540 workMV = MVT::Clone( *problem_->getInitVec(), (nevBlocks_+maxXtraBlocks)*blockSize_ + 1 + blockSize_ );
544 workMV = Teuchos::null;
548 Eigensolution<ScalarType,MV> sol;
550 problem_->setSolution(sol);
553 int cur_nevBlocks = 0;
557#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
558 Teuchos::TimeMonitor slvtimer(*timerSolve_);
564 bks_solver->iterate();
571 if ( ordertest->getStatus() ==
Passed ) {
589 else if ( (bks_solver->getCurSubspaceDim() == bks_solver->getMaxSubspaceDim()) ||
590 (!problem_->isHermitian() && !conjSplit_ && (bks_solver->getCurSubspaceDim()+1 == bks_solver->getMaxSubspaceDim())) ) {
593 if (!bks_solver->isSchurCurrent()) {
594 bks_solver->computeSchurForm(
true );
597 outputtest->checkStatus( &*bks_solver );
601 if ( numRestarts >= maxRestarts_ || ordertest->getStatus() ==
Passed) {
606#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
607 Teuchos::TimeMonitor restimer(*timerRestarting_);
611 int numConv = ordertest->howMany();
612 cur_nevBlocks = nevBlocks_*blockSize_;
615 int moreNevBlocks = std::min( maxXtraBlocks, std::max( numConv/blockSize_, xtra_nevBlocks_) );
617 cur_nevBlocks += moreNevBlocks * blockSize_;
618 else if ( xtra_nevBlocks_ )
619 cur_nevBlocks += xtra_nevBlocks_ * blockSize_;
627 printer->stream(
Debug) <<
" Performing restart number " << numRestarts <<
" of " << maxRestarts_ << std::endl << std::endl;
628 printer->stream(
Debug) <<
" - Current NEV blocks is " << cur_nevBlocks <<
", the minimum is " << nevBlocks_*blockSize_ << std::endl;
631 ritzValues_ = bks_solver->getRitzValues();
634 BlockKrylovSchurState<ScalarType,MV> oldState = bks_solver->getState();
637 int curDim = oldState.curDim;
640 std::vector<int> ritzIndex = bks_solver->getRitzIndex();
641 if (ritzIndex[cur_nevBlocks-1]==1) {
650 if (problem_->isHermitian() && conjSplit_)
653 <<
" Eigenproblem is Hermitian, complex eigenvalues have been detected, and eigenvalues of interest split a conjugate pair!!!"
655 <<
" Block Krylov-Schur eigensolver cannot guarantee correct behavior in this situation, please turn Hermitian flag off!!!"
662 Teuchos::SerialDenseMatrix<int,ScalarType> Qnev(Teuchos::View, *(oldState.Q), curDim, cur_nevBlocks);
665 std::vector<int> curind( curDim );
666 for (
int i=0; i<curDim; i++) { curind[i] = i; }
667 Teuchos::RCP<const MV> basistemp = MVT::CloneView( *(oldState.V), curind );
677 Teuchos::RCP<MV> newF;
678 if (inSituRestart_) {
681 Teuchos::RCP<MV> solverbasis = Teuchos::rcp_const_cast<MV>(oldState.V);
682 Teuchos::SerialDenseMatrix<int,ScalarType> copyQnev(Teuchos::Copy, Qnev);
685 std::vector<ScalarType> tau(cur_nevBlocks), work(cur_nevBlocks);
687 lapack.GEQRF(curDim,cur_nevBlocks,copyQnev.values(),copyQnev.stride(),&tau[0],&work[0],work.size(),&info);
688 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,
689 "Anasazi::BlockKrylovSchurSolMgr::solve(): error calling GEQRF during restarting.");
691 std::vector<ScalarType> d(cur_nevBlocks);
692 for (
int j=0; j<copyQnev.numCols(); j++) {
693 d[j] = copyQnev(j,j);
695 if (printer->isVerbosity(
Debug)) {
696 Teuchos::SerialDenseMatrix<int,ScalarType> R(Teuchos::Copy,copyQnev,cur_nevBlocks,cur_nevBlocks);
697 for (
int j=0; j<R.numCols(); j++) {
698 R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
699 for (
int i=j+1; i<R.numRows(); i++) {
703 printer->stream(
Debug) <<
"||Triangular factor of Su - I||: " << R.normFrobenius() << std::endl;
709 curind.resize(curDim);
710 for (
int i=0; i<curDim; i++) curind[i] = i;
712 Teuchos::RCP<MV> oldV = MVT::CloneViewNonConst(*solverbasis,curind);
717 curind.resize(cur_nevBlocks);
718 for (
int i=0; i<cur_nevBlocks; i++) { curind[i] = i; }
720 Teuchos::RCP<MV> newV = MVT::CloneViewNonConst( *solverbasis, curind );
721 MVT::MvScale(*newV,d);
724 curind.resize(blockSize_);
725 for (
int i=0; i<blockSize_; i++) { curind[i] = cur_nevBlocks + i; }
726 newF = MVT::CloneViewNonConst( *solverbasis, curind );
730 curind.resize(cur_nevBlocks);
731 for (
int i=0; i<cur_nevBlocks; i++) { curind[i] = i; }
732 Teuchos::RCP<MV> tmp_newV = MVT::CloneViewNonConst(*workMV, curind );
734 MVT::MvTimesMatAddMv( one, *basistemp, Qnev, zero, *tmp_newV );
735 tmp_newV = Teuchos::null;
737 curind.resize(blockSize_);
738 for (
int i=0; i<blockSize_; i++) { curind[i] = cur_nevBlocks + i; }
739 newF = MVT::CloneViewNonConst( *workMV, curind );
743 curind.resize(blockSize_);
744 for (
int i=0; i<blockSize_; i++) { curind[i] = curDim + i; }
745 Teuchos::RCP<const MV> oldF = MVT::CloneView( *(oldState.V), curind );
746 for (
int i=0; i<blockSize_; i++) { curind[i] = i; }
747 MVT::SetBlock( *oldF, curind, *newF );
748 newF = Teuchos::null;
754 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > newH =
755 Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::Copy, *(oldState.S), cur_nevBlocks+blockSize_, cur_nevBlocks) );
758 Teuchos::SerialDenseMatrix<int,ScalarType> oldB(Teuchos::View, *(oldState.H), blockSize_, blockSize_, curDim, curDim-blockSize_);
761 Teuchos::SerialDenseMatrix<int,ScalarType> subQ(Teuchos::View, *(oldState.Q), blockSize_, cur_nevBlocks, curDim-blockSize_);
764 Teuchos::SerialDenseMatrix<int,ScalarType> newB(Teuchos::View, *newH, blockSize_, cur_nevBlocks, cur_nevBlocks);
767 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, blockSize_, cur_nevBlocks, blockSize_, one,
768 oldB.values(), oldB.stride(), subQ.values(), subQ.stride(), zero, newB.values(), newB.stride() );
773 BlockKrylovSchurState<ScalarType,MV> newstate;
774 if (inSituRestart_) {
775 newstate.V = oldState.V;
780 newstate.curDim = cur_nevBlocks;
781 bks_solver->initialize(newstate);
791 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Anasazi::BlockKrylovSchurSolMgr::solve(): Invalid return from bks_solver::iterate().");
796 <<
"Anasazi::BlockKrylovSchurSolMgr::solve() caught unexpected exception from Anasazi::BlockKrylovSchur::iterate() at iteration " << bks_solver->getNumIters() << std::endl
797 << err.what() << std::endl
798 <<
"Anasazi::BlockKrylovSchurSolMgr::solve() returning Unconverged with no solutions." << std::endl;
805 workMV = Teuchos::null;
808 ritzValues_ = bks_solver->getRitzValues();
810 sol.numVecs = ordertest->howMany();
811 printer->stream(
Debug) <<
"ordertest->howMany() : " << sol.numVecs << std::endl;
812 std::vector<int> whichVecs = ordertest->whichVecs();
815 if (sol.numVecs > 0) {
818 std::vector<int> tmpIndex = bks_solver->getRitzIndex();
819 for (
int i=0; i<(int)ritzValues_.size(); ++i) {
820 printer->stream(
Debug) << ritzValues_[i].realpart <<
" + i " << ritzValues_[i].imagpart <<
", Index = " << tmpIndex[i] << std::endl;
822 printer->stream(
Debug) <<
"Number of converged eigenpairs (before) = " << sol.numVecs << std::endl;
823 for (
int i=0; i<sol.numVecs; ++i) {
824 printer->stream(
Debug) <<
"whichVecs[" << i <<
"] = " << whichVecs[i] <<
", tmpIndex[" << whichVecs[i] <<
"] = " << tmpIndex[whichVecs[i]] << std::endl;
826 if (tmpIndex[whichVecs[sol.numVecs-1]]==1) {
827 printer->stream(
Debug) <<
"There is a conjugate pair on the boundary, resizing sol.numVecs" << std::endl;
828 whichVecs.push_back(whichVecs[sol.numVecs-1]+1);
830 for (
int i=0; i<sol.numVecs; ++i) {
831 printer->stream(
Debug) <<
"whichVecs[" << i <<
"] = " << whichVecs[i] <<
", tmpIndex[" << whichVecs[i] <<
"] = " << tmpIndex[whichVecs[i]] << std::endl;
835 bool keepMore =
false;
836 int numEvecs = sol.numVecs;
837 printer->stream(
Debug) <<
"Number of converged eigenpairs (after) = " << sol.numVecs << std::endl;
838 printer->stream(
Debug) <<
"whichVecs[sol.numVecs-1] > sol.numVecs-1 : " << whichVecs[sol.numVecs-1] <<
" > " << sol.numVecs-1 << std::endl;
839 if (whichVecs[sol.numVecs-1] > (sol.numVecs-1)) {
841 numEvecs = whichVecs[sol.numVecs-1]+1;
842 printer->stream(
Debug) <<
"keepMore = true; numEvecs = " << numEvecs << std::endl;
846 bks_solver->setNumRitzVectors(numEvecs);
847 bks_solver->computeRitzVectors();
853 sol.index = bks_solver->getRitzIndex();
854 sol.Evals = bks_solver->getRitzValues();
855 sol.Evecs = MVT::CloneCopy( *(bks_solver->getRitzVectors()) );
860 sol.Evals.resize(sol.numVecs);
861 sol.index.resize(sol.numVecs);
865 std::vector<Anasazi::Value<ScalarType> > tmpEvals = bks_solver->getRitzValues();
866 for (
int vec_i=0; vec_i<sol.numVecs; ++vec_i) {
867 sol.index[vec_i] = tmpIndex[whichVecs[vec_i]];
868 sol.Evals[vec_i] = tmpEvals[whichVecs[vec_i]];
870 sol.Evecs = MVT::CloneCopy( *(bks_solver->getRitzVectors()), whichVecs );
874 sol.Espace = sol.Evecs;
879 bks_solver->currentStatus(printer->stream(
FinalSummary));
882#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
884 Teuchos::TimeMonitor::summarize( printer->stream(
TimingDetails ) );
888 problem_->setSolution(sol);
889 printer->stream(
Debug) <<
"Returning " << sol.numVecs <<
" eigenpairs to eigenproblem." << std::endl;
892 numIters_ = bks_solver->getNumIters();
894 if (sol.numVecs < nev) {