268 Teuchos::ParameterList &pl ) :
284 inSituRestart_(false),
288#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
289 ,timerSolve_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockKrylovSchurSolMgr::solve()")),
290 timerRestarting_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockKrylovSchurSolMgr restarting"))
293 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument,
"Problem not given to solver manager.");
294 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument,
"Problem not set.");
295 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null, std::invalid_argument,
"Problem does not contain initial vectors to clone from.");
297 const int nev = problem_->getNEV();
300 convtol_ = pl.get(
"Convergence Tolerance",MT::prec());
301 relconvtol_ = pl.get(
"Relative Convergence Tolerance",relconvtol_);
304 maxRestarts_ = pl.get(
"Maximum Restarts",maxRestarts_);
307 blockSize_ = pl.get(
"Block Size",1);
308 TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
309 "Anasazi::BlockKrylovSchurSolMgr: \"Block Size\" must be strictly positive.");
312 xtra_nevBlocks_ = pl.get(
"Extra NEV Blocks",0);
313 if (nev%blockSize_) {
314 nevBlocks_ = nev/blockSize_ + 1;
316 nevBlocks_ = nev/blockSize_;
322 if (pl.isParameter(
"Dynamic Extra NEV")) {
323 if (Teuchos::isParameterType<bool>(pl,
"Dynamic Extra NEV")) {
324 dynXtraNev_ = pl.get(
"Dynamic Extra NEV",dynXtraNev_);
326 dynXtraNev_ = ( Teuchos::getParameter<int>(pl,
"Dynamic Extra NEV") != 0 );
330 numBlocks_ = pl.get(
"Num Blocks",3*nevBlocks_);
331 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= nevBlocks_, std::invalid_argument,
332 "Anasazi::BlockKrylovSchurSolMgr: \"Num Blocks\" must be strictly positive and large enough to compute the requested eigenvalues.");
334 TEUCHOS_TEST_FOR_EXCEPTION(
static_cast<ptrdiff_t
>(numBlocks_)*blockSize_ > MVT::GetGlobalLength(*problem_->getInitVec()),
335 std::invalid_argument,
336 "Anasazi::BlockKrylovSchurSolMgr: Potentially impossible orthogonality requests. Reduce basis size.");
340 stepSize_ = pl.get(
"Step Size", (maxRestarts_+1)*(numBlocks_+1));
342 stepSize_ = pl.get(
"Step Size", numBlocks_+1);
344 TEUCHOS_TEST_FOR_EXCEPTION(stepSize_ < 1, std::invalid_argument,
345 "Anasazi::BlockKrylovSchurSolMgr: \"Step Size\" must be strictly positive.");
348 if (pl.isParameter(
"Sort Manager")) {
349 sort_ = Teuchos::getParameter<Teuchos::RCP<Anasazi::SortManager<MagnitudeType> > >(pl,
"Sort Manager");
352 whch_ = pl.get(
"Which",whch_);
353 TEUCHOS_TEST_FOR_EXCEPTION(whch_ !=
"SM" && whch_ !=
"LM" && whch_ !=
"SR" && whch_ !=
"LR" && whch_ !=
"SI" && whch_ !=
"LI",
354 std::invalid_argument,
"Invalid sorting string.");
359 ortho_ = pl.get(
"Orthogonalization",ortho_);
360 if (ortho_ !=
"DGKS" && ortho_ !=
"SVQB" && ortho_ !=
"ICGS") {
365 ortho_kappa_ = pl.get(
"Orthogonalization Constant",ortho_kappa_);
369 osProc_ = pl.get(
"Output Processor", osProc_);
372 if (pl.isParameter(
"Output Stream")) {
373 osp_ = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,
"Output Stream");
380 if (pl.isParameter(
"Verbosity")) {
381 if (Teuchos::isParameterType<int>(pl,
"Verbosity")) {
382 verbosity_ = pl.get(
"Verbosity", verbosity_);
384 verbosity_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,
"Verbosity");
389 if (pl.isParameter(
"In Situ Restarting")) {
390 if (Teuchos::isParameterType<bool>(pl,
"In Situ Restarting")) {
391 inSituRestart_ = pl.get(
"In Situ Restarting",inSituRestart_);
393 inSituRestart_ = ( Teuchos::getParameter<int>(pl,
"In Situ Restarting") != 0 );
397 printNum_ = pl.get<
int>(
"Print Number of Ritz Values",printNum_);
406 const int nev = problem_->getNEV();
407 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
408 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
410 Teuchos::BLAS<int,ScalarType> blas;
411 Teuchos::LAPACK<int,ScalarType> lapack;
421 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convtest;
422 if (globalTest_ == Teuchos::null) {
426 convtest = globalTest_;
428 Teuchos::RCP<StatusTestWithOrdering<ScalarType,MV,OP> > ordertest
431 Teuchos::Array<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
432 alltests.push_back(ordertest);
434 if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
436 Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > combotest
439 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
440 if ( printer->isVerbosity(
Debug) ) {
449 Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho;
450 if (ortho_==
"SVQB") {
452 }
else if (ortho_==
"DGKS") {
453 if (ortho_kappa_ <= 0) {
458 }
else if (ortho_==
"ICGS") {
461 TEUCHOS_TEST_FOR_EXCEPTION(ortho_!=
"SVQB"&&ortho_!=
"DGKS"&&ortho_!=
"ICGS",std::logic_error,
"Anasazi::BlockKrylovSchurSolMgr::solve(): Invalid orthogonalization type.");
466 Teuchos::ParameterList plist;
467 plist.set(
"Block Size",blockSize_);
468 plist.set(
"Num Blocks",numBlocks_);
469 plist.set(
"Step Size",stepSize_);
470 if (printNum_ == -1) {
471 plist.set(
"Print Number of Ritz Values",nevBlocks_*blockSize_);
474 plist.set(
"Print Number of Ritz Values",printNum_);
479 Teuchos::RCP<BlockKrylovSchur<ScalarType,MV,OP> > bks_solver
482 Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
483 if (probauxvecs != Teuchos::null) {
484 bks_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
494 int maxXtraBlocks = 0;
495 if ( dynXtraNev_ ) maxXtraBlocks = ( ( bks_solver->getMaxSubspaceDim() - nev ) / blockSize_ ) / 2;
497 Teuchos::RCP<MV> workMV;
498 if (maxRestarts_ > 0) {
499 if (inSituRestart_==
true) {
501 workMV = MVT::Clone( *problem_->getInitVec(), 1 );
504 if (problem_->isHermitian()) {
505 workMV = MVT::Clone( *problem_->getInitVec(), (nevBlocks_+maxXtraBlocks)*blockSize_ + blockSize_ );
508 workMV = MVT::Clone( *problem_->getInitVec(), (nevBlocks_+maxXtraBlocks)*blockSize_ + 1 + blockSize_ );
512 workMV = Teuchos::null;
518 problem_->setSolution(sol);
521 int cur_nevBlocks = 0;
525#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
526 Teuchos::TimeMonitor slvtimer(*timerSolve_);
532 bks_solver->iterate();
539 if ( ordertest->getStatus() ==
Passed ) {
557 else if ( (bks_solver->getCurSubspaceDim() == bks_solver->getMaxSubspaceDim()) ||
558 (!problem_->isHermitian() && !conjSplit_ && (bks_solver->getCurSubspaceDim()+1 == bks_solver->getMaxSubspaceDim())) ) {
561 if (!bks_solver->isSchurCurrent()) {
562 bks_solver->computeSchurForm(
true );
565 outputtest->checkStatus( &*bks_solver );
569 if ( numRestarts >= maxRestarts_ || ordertest->getStatus() ==
Passed) {
574#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
575 Teuchos::TimeMonitor restimer(*timerRestarting_);
579 int numConv = ordertest->howMany();
580 cur_nevBlocks = nevBlocks_*blockSize_;
583 int moreNevBlocks = std::min( maxXtraBlocks, std::max( numConv/blockSize_, xtra_nevBlocks_) );
585 cur_nevBlocks += moreNevBlocks * blockSize_;
586 else if ( xtra_nevBlocks_ )
587 cur_nevBlocks += xtra_nevBlocks_ * blockSize_;
595 printer->stream(
Debug) <<
" Performing restart number " << numRestarts <<
" of " << maxRestarts_ << std::endl << std::endl;
596 printer->stream(
Debug) <<
" - Current NEV blocks is " << cur_nevBlocks <<
", the minimum is " << nevBlocks_*blockSize_ << std::endl;
599 ritzValues_ = bks_solver->getRitzValues();
605 int curDim = oldState.
curDim;
608 std::vector<int> ritzIndex = bks_solver->getRitzIndex();
609 if (ritzIndex[cur_nevBlocks-1]==1) {
618 if (problem_->isHermitian() && conjSplit_)
621 <<
" Eigenproblem is Hermitian, complex eigenvalues have been detected, and eigenvalues of interest split a conjugate pair!!!"
623 <<
" Block Krylov-Schur eigensolver cannot guarantee correct behavior in this situation, please turn Hermitian flag off!!!"
630 Teuchos::SerialDenseMatrix<int,ScalarType> Qnev(Teuchos::View, *(oldState.
Q), curDim, cur_nevBlocks);
633 std::vector<int> curind( curDim );
634 for (
int i=0; i<curDim; i++) { curind[i] = i; }
635 Teuchos::RCP<const MV> basistemp = MVT::CloneView( *(oldState.
V), curind );
645 Teuchos::RCP<MV> newF;
646 if (inSituRestart_) {
649 Teuchos::RCP<MV> solverbasis = Teuchos::rcp_const_cast<MV>(oldState.
V);
650 Teuchos::SerialDenseMatrix<int,ScalarType> copyQnev(Teuchos::Copy, Qnev);
653 std::vector<ScalarType> tau(cur_nevBlocks), work(cur_nevBlocks);
655 lapack.GEQRF(curDim,cur_nevBlocks,copyQnev.values(),copyQnev.stride(),&tau[0],&work[0],work.size(),&info);
656 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,
657 "Anasazi::BlockKrylovSchurSolMgr::solve(): error calling GEQRF during restarting.");
659 std::vector<ScalarType> d(cur_nevBlocks);
660 for (
int j=0; j<copyQnev.numCols(); j++) {
661 d[j] = copyQnev(j,j);
663 if (printer->isVerbosity(
Debug)) {
664 Teuchos::SerialDenseMatrix<int,ScalarType> R(Teuchos::Copy,copyQnev,cur_nevBlocks,cur_nevBlocks);
665 for (
int j=0; j<R.numCols(); j++) {
666 R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
667 for (
int i=j+1; i<R.numRows(); i++) {
671 printer->stream(
Debug) <<
"||Triangular factor of Su - I||: " << R.normFrobenius() << std::endl;
677 curind.resize(curDim);
678 for (
int i=0; i<curDim; i++) curind[i] = i;
680 Teuchos::RCP<MV> oldV = MVT::CloneViewNonConst(*solverbasis,curind);
685 curind.resize(cur_nevBlocks);
686 for (
int i=0; i<cur_nevBlocks; i++) { curind[i] = i; }
688 Teuchos::RCP<MV> newV = MVT::CloneViewNonConst( *solverbasis, curind );
689 MVT::MvScale(*newV,d);
692 curind.resize(blockSize_);
693 for (
int i=0; i<blockSize_; i++) { curind[i] = cur_nevBlocks + i; }
694 newF = MVT::CloneViewNonConst( *solverbasis, curind );
698 curind.resize(cur_nevBlocks);
699 for (
int i=0; i<cur_nevBlocks; i++) { curind[i] = i; }
700 Teuchos::RCP<MV> tmp_newV = MVT::CloneViewNonConst(*workMV, curind );
702 MVT::MvTimesMatAddMv( one, *basistemp, Qnev, zero, *tmp_newV );
703 tmp_newV = Teuchos::null;
705 curind.resize(blockSize_);
706 for (
int i=0; i<blockSize_; i++) { curind[i] = cur_nevBlocks + i; }
707 newF = MVT::CloneViewNonConst( *workMV, curind );
711 curind.resize(blockSize_);
712 for (
int i=0; i<blockSize_; i++) { curind[i] = curDim + i; }
713 Teuchos::RCP<const MV> oldF = MVT::CloneView( *(oldState.
V), curind );
714 for (
int i=0; i<blockSize_; i++) { curind[i] = i; }
715 MVT::SetBlock( *oldF, curind, *newF );
716 newF = Teuchos::null;
722 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > newH =
723 Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::Copy, *(oldState.
S), cur_nevBlocks+blockSize_, cur_nevBlocks) );
726 Teuchos::SerialDenseMatrix<int,ScalarType> oldB(Teuchos::View, *(oldState.
H), blockSize_, blockSize_, curDim, curDim-blockSize_);
729 Teuchos::SerialDenseMatrix<int,ScalarType> subQ(Teuchos::View, *(oldState.
Q), blockSize_, cur_nevBlocks, curDim-blockSize_);
732 Teuchos::SerialDenseMatrix<int,ScalarType> newB(Teuchos::View, *newH, blockSize_, cur_nevBlocks, cur_nevBlocks);
735 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, blockSize_, cur_nevBlocks, blockSize_, one,
736 oldB.values(), oldB.stride(), subQ.values(), subQ.stride(), zero, newB.values(), newB.stride() );
742 if (inSituRestart_) {
743 newstate.
V = oldState.
V;
748 newstate.
curDim = cur_nevBlocks;
749 bks_solver->initialize(newstate);
759 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Anasazi::BlockKrylovSchurSolMgr::solve(): Invalid return from bks_solver::iterate().");
764 <<
"Anasazi::BlockKrylovSchurSolMgr::solve() caught unexpected exception from Anasazi::BlockKrylovSchur::iterate() at iteration " << bks_solver->getNumIters() << std::endl
765 << err.what() << std::endl
766 <<
"Anasazi::BlockKrylovSchurSolMgr::solve() returning Unconverged with no solutions." << std::endl;
773 workMV = Teuchos::null;
776 ritzValues_ = bks_solver->getRitzValues();
778 sol.
numVecs = ordertest->howMany();
779 printer->stream(
Debug) <<
"ordertest->howMany() : " << sol.
numVecs << std::endl;
780 std::vector<int> whichVecs = ordertest->whichVecs();
786 std::vector<int> tmpIndex = bks_solver->getRitzIndex();
787 for (
int i=0; i<(int)ritzValues_.size(); ++i) {
788 printer->stream(
Debug) << ritzValues_[i].realpart <<
" + i " << ritzValues_[i].imagpart <<
", Index = " << tmpIndex[i] << std::endl;
790 printer->stream(
Debug) <<
"Number of converged eigenpairs (before) = " << sol.
numVecs << std::endl;
791 for (
int i=0; i<sol.
numVecs; ++i) {
792 printer->stream(
Debug) <<
"whichVecs[" << i <<
"] = " << whichVecs[i] <<
", tmpIndex[" << whichVecs[i] <<
"] = " << tmpIndex[whichVecs[i]] << std::endl;
794 if (tmpIndex[whichVecs[sol.
numVecs-1]]==1) {
795 printer->stream(
Debug) <<
"There is a conjugate pair on the boundary, resizing sol.numVecs" << std::endl;
796 whichVecs.push_back(whichVecs[sol.
numVecs-1]+1);
798 for (
int i=0; i<sol.
numVecs; ++i) {
799 printer->stream(
Debug) <<
"whichVecs[" << i <<
"] = " << whichVecs[i] <<
", tmpIndex[" << whichVecs[i] <<
"] = " << tmpIndex[whichVecs[i]] << std::endl;
803 bool keepMore =
false;
805 printer->stream(
Debug) <<
"Number of converged eigenpairs (after) = " << sol.
numVecs << std::endl;
806 printer->stream(
Debug) <<
"whichVecs[sol.numVecs-1] > sol.numVecs-1 : " << whichVecs[sol.
numVecs-1] <<
" > " << sol.
numVecs-1 << std::endl;
809 numEvecs = whichVecs[sol.
numVecs-1]+1;
810 printer->stream(
Debug) <<
"keepMore = true; numEvecs = " << numEvecs << std::endl;
814 bks_solver->setNumRitzVectors(numEvecs);
815 bks_solver->computeRitzVectors();
821 sol.
index = bks_solver->getRitzIndex();
822 sol.
Evals = bks_solver->getRitzValues();
823 sol.
Evecs = MVT::CloneCopy( *(bks_solver->getRitzVectors()) );
833 std::vector<Anasazi::Value<ScalarType> > tmpEvals = bks_solver->getRitzValues();
834 for (
int vec_i=0; vec_i<sol.
numVecs; ++vec_i) {
835 sol.
index[vec_i] = tmpIndex[whichVecs[vec_i]];
836 sol.
Evals[vec_i] = tmpEvals[whichVecs[vec_i]];
838 sol.
Evecs = MVT::CloneCopy( *(bks_solver->getRitzVectors()), whichVecs );
847 bks_solver->currentStatus(printer->stream(
FinalSummary));
850#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
852 Teuchos::TimeMonitor::summarize( printer->stream(
TimingDetails ) );
856 problem_->setSolution(sol);
857 printer->stream(
Debug) <<
"Returning " << sol.
numVecs <<
" eigenpairs to eigenproblem." << std::endl;
860 numIters_ = bks_solver->getNumIters();