175 typedef Teuchos::ScalarTraits<ScalarType> SCT;
176 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
177 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
227 const Teuchos::RCP<Teuchos::ParameterList> &pl );
233 virtual Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const {
249 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters()
const;
260 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
261 return Teuchos::tuple(timerSolve_);
291 void setParameters(
const Teuchos::RCP<Teuchos::ParameterList> ¶ms );
332 std::string description()
const;
340 int ARRQR(
int numVecs,
int numOrthVecs,
const Teuchos::SerialDenseMatrix<int,ScalarType>& D);
343 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
346 Teuchos::RCP<OutputManager<ScalarType> > printer_;
347 Teuchos::RCP<std::ostream> outputStream_;
350 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
351 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
352 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
353 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
356 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
359 Teuchos::RCP<Teuchos::ParameterList> params_;
362 static constexpr int maxIters_default_ = 1000;
363 static constexpr int deflatedBlocks_default_ = 2;
364 static constexpr int savedBlocks_default_ = 16;
367 static constexpr int outputFreq_default_ = -1;
368 static constexpr const char * label_default_ =
"Belos";
369 static constexpr const char * orthoType_default_ =
"ICGS";
376 MagnitudeType convtol_;
379 MagnitudeType orthoKappa_;
382 MagnitudeType achievedTol_;
390 int deflatedBlocks_, savedBlocks_, verbosity_, outputStyle_, outputFreq_;
391 std::string orthoType_;
394 Teuchos::RCP<MV> U_, C_, R_;
401 Teuchos::RCP<Teuchos::Time> timerSolve_;
469 if (params_ == Teuchos::null) {
470 params_ = Teuchos::rcp(
new Teuchos::ParameterList(*getValidParameters()) );
473 params->validateParameters(*getValidParameters());
477 if (params->isParameter(
"Maximum Iterations")) {
478 maxIters_ = params->get(
"Maximum Iterations",maxIters_default_);
481 params_->set(
"Maximum Iterations", maxIters_);
482 if (maxIterTest_!=Teuchos::null)
483 maxIterTest_->setMaxIters( maxIters_ );
487 if (params->isParameter(
"Num Saved Blocks")) {
488 savedBlocks_ = params->get(
"Num Saved Blocks",savedBlocks_default_);
489 TEUCHOS_TEST_FOR_EXCEPTION(savedBlocks_ <= 0, std::invalid_argument,
490 "Belos::PCPGSolMgr: \"Num Saved Blocks\" must be strictly positive.");
497 params_->set(
"Num Saved Blocks", savedBlocks_);
499 if (params->isParameter(
"Num Deflated Blocks")) {
500 deflatedBlocks_ = params->get(
"Num Deflated Blocks",deflatedBlocks_default_);
501 TEUCHOS_TEST_FOR_EXCEPTION(deflatedBlocks_ < 0, std::invalid_argument,
502 "Belos::PCPGSolMgr: \"Num Deflated Blocks\" must be positive.");
504 TEUCHOS_TEST_FOR_EXCEPTION(deflatedBlocks_ > savedBlocks_, std::invalid_argument,
505 "Belos::PCPGSolMgr: \"Num Deflated Blocks\" must be <= \"Num Saved Blocks\".");
509 params_->set(
"Num Deflated Blocks",
static_cast<int>(deflatedBlocks_));
513 if (params->isParameter(
"Timer Label")) {
514 std::string tempLabel = params->get(
"Timer Label", label_default_);
517 if (tempLabel != label_) {
519 params_->set(
"Timer Label", label_);
520 std::string solveLabel = label_ +
": PCPGSolMgr total solve time";
521#ifdef BELOS_TEUCHOS_TIME_MONITOR
522 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
524 if (ortho_ != Teuchos::null) {
525 ortho_->setLabel( label_ );
531 if (params->isParameter(
"Verbosity")) {
532 if (Teuchos::isParameterType<int>(*params,
"Verbosity")) {
533 verbosity_ = params->get(
"Verbosity", verbosity_default_);
535 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,
"Verbosity");
539 params_->set(
"Verbosity", verbosity_);
540 if (printer_ != Teuchos::null)
541 printer_->setVerbosity(verbosity_);
545 if (params->isParameter(
"Output Style")) {
546 if (Teuchos::isParameterType<int>(*params,
"Output Style")) {
547 outputStyle_ = params->get(
"Output Style", outputStyle_default_);
549 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,
"Output Style");
553 params_->set(
"Output Style", outputStyle_);
554 outputTest_ = Teuchos::null;
558 if (params->isParameter(
"Output Stream")) {
559 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,
"Output Stream");
562 params_->set(
"Output Stream", outputStream_);
563 if (printer_ != Teuchos::null)
564 printer_->setOStream( outputStream_ );
569 if (params->isParameter(
"Output Frequency")) {
570 outputFreq_ = params->get(
"Output Frequency", outputFreq_default_);
574 params_->set(
"Output Frequency", outputFreq_);
575 if (outputTest_ != Teuchos::null)
576 outputTest_->setOutputFrequency( outputFreq_ );
580 if (printer_ == Teuchos::null) {
585 bool changedOrthoType =
false;
586 if (params->isParameter(
"Orthogonalization")) {
587 std::string tempOrthoType = params->get(
"Orthogonalization",orthoType_default_);
588 if (tempOrthoType != orthoType_) {
589 orthoType_ = tempOrthoType;
590 changedOrthoType =
true;
593 params_->set(
"Orthogonalization", orthoType_);
596 if (params->isParameter(
"Orthogonalization Constant")) {
597 if (params->isType<MagnitudeType> (
"Orthogonalization Constant")) {
598 orthoKappa_ = params->get (
"Orthogonalization Constant",
602 orthoKappa_ = params->get (
"Orthogonalization Constant",
607 params_->set(
"Orthogonalization Constant",orthoKappa_);
608 if (orthoType_==
"DGKS") {
609 if (orthoKappa_ > 0 && ortho_ != Teuchos::null && !changedOrthoType) {
610 Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
616 if (ortho_ == Teuchos::null || changedOrthoType) {
618 Teuchos::RCP<Teuchos::ParameterList> paramsOrtho;
619 if (orthoType_==
"DGKS" && orthoKappa_ > 0) {
620 paramsOrtho->set (
"depTol", orthoKappa_ );
623 ortho_ = factory.
makeMatOrthoManager (orthoType_, Teuchos::null, printer_, label_, paramsOrtho);
631 if (params->isParameter(
"Convergence Tolerance")) {
632 if (params->isType<MagnitudeType> (
"Convergence Tolerance")) {
633 convtol_ = params->get (
"Convergence Tolerance",
641 params_->set(
"Convergence Tolerance", convtol_);
642 if (convTest_ != Teuchos::null)
649 if (maxIterTest_ == Teuchos::null)
652 if (convTest_ == Teuchos::null)
653 convTest_ = Teuchos::rcp(
new StatusTestResNorm_t( convtol_, 1 ) );
655 sTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
663 std::string solverDesc =
" PCPG ";
664 outputTest_->setSolverDesc( solverDesc );
667 if (timerSolve_ == Teuchos::null) {
668 std::string solveLabel = label_ +
": PCPGSolMgr total solve time";
669#ifdef BELOS_TEUCHOS_TIME_MONITOR
670 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
727 if (!isSet_) { setParameters( params_ ); }
729 Teuchos::LAPACK<int,ScalarType> lapack;
730 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
731 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
734 "Belos::PCPGSolMgr::solve(): Linear problem is not a valid object.");
737 "Belos::PCPGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
740 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
741 std::vector<int> currIdx(1);
747 problem_->setLSIndex( currIdx );
750 bool isConverged =
true;
754 Teuchos::ParameterList plist;
755 plist.set(
"Saved Blocks", savedBlocks_);
756 plist.set(
"Block Size", 1);
757 plist.set(
"Keep Diagonal",
true);
758 plist.set(
"Initialize Diagonal",
true);
763 Teuchos::RCP<PCPGIter<ScalarType,MV,OP> > pcpg_iter;
769#ifdef BELOS_TEUCHOS_TIME_MONITOR
770 Teuchos::TimeMonitor slvtimer(*timerSolve_);
772 while ( numRHS2Solve > 0 ) {
775 outputTest_->reset();
778 if (R_ == Teuchos::null)
779 R_ = MVT::Clone( *(problem_->getRHS()), 1 );
781 problem_->computeCurrResVec( &*R_ );
787 if( U_ != Teuchos::null ){
792 Teuchos::RCP<MV> cur_soln_vec = problem_->getCurrLHSVec();
793 std::vector<MagnitudeType> rnorm0(1);
794 MVT::MvNorm( *R_, rnorm0 );
797 std::cout <<
"Solver Manager: dimU_ = " << dimU_ << std::endl;
798 Teuchos::SerialDenseMatrix<int,ScalarType> Z( dimU_, 1 );
800 Teuchos::RCP<const MV> Uactive, Cactive;
801 std::vector<int> active_columns( dimU_ );
802 for (
int i=0; i < dimU_; ++i) active_columns[i] = i;
803 Uactive = MVT::CloneView(*U_, active_columns);
804 Cactive = MVT::CloneView(*C_, active_columns);
807 std::cout <<
" Solver Manager : check duality of seed basis " << std::endl;
808 Teuchos::SerialDenseMatrix<int,ScalarType> H( dimU_, dimU_ );
809 MVT::MvTransMv( one, *Uactive, *Cactive, H );
810 H.print( std::cout );
813 MVT::MvTransMv( one, *Uactive, *R_, Z );
814 Teuchos::RCP<MV> tempU = MVT::Clone( *R_, 1 );
815 MVT::MvTimesMatAddMv( one, *Uactive, Z, zero, *tempU );
816 MVT::MvAddMv( one, *tempU, one, *cur_soln_vec, *cur_soln_vec );
817 MVT::MvTimesMatAddMv( one, *Cactive, Z, zero, *tempU );
818 MVT::MvAddMv( -one, *tempU, one, *R_, *R_ );
819 std::vector<MagnitudeType> rnorm(1);
820 MVT::MvNorm( *R_, rnorm );
821 if( rnorm[0] < rnorm0[0] * .001 ){
822 MVT::MvTransMv( one, *Uactive, *R_, Z );
823 MVT::MvTimesMatAddMv( one, *Uactive, Z, zero, *tempU );
824 MVT::MvAddMv( one, *tempU, one, *cur_soln_vec, *cur_soln_vec );
825 MVT::MvTimesMatAddMv( one, *Cactive, Z, zero, *tempU );
826 MVT::MvAddMv( -one, *tempU, one, *R_, *R_ );
828 Uactive = Teuchos::null;
829 Cactive = Teuchos::null;
830 tempU = Teuchos::null;
841 if( U_ != Teuchos::null ) pcpgState.
U = U_;
842 if( C_ != Teuchos::null ) pcpgState.
C = C_;
843 if( dimU_ > 0 ) pcpgState.
curDim = dimU_;
844 pcpg_iter->initialize(pcpgState);
850 if( !dimU_ ) printer_->stream(
Debug) <<
" No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
851 pcpg_iter->resetNumIters();
853 if( dimU_ > savedBlocks_ )
854 std::cout <<
"Error: dimU_ = " << dimU_ <<
" > savedBlocks_ = " << savedBlocks_ << std::endl;
860 if( debug ) printf(
"********** Calling iterate...\n");
861 pcpg_iter->iterate();
868 if ( convTest_->getStatus() ==
Passed ) {
877 else if ( maxIterTest_->getStatus() ==
Passed ) {
891 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
892 "Belos::PCPGSolMgr::solve(): Invalid return from PCPGIter::iterate().");
895 catch (
const std::exception &e) {
896 printer_->stream(
Errors) <<
"Error! Caught exception in PCPGIter::iterate() at iteration "
897 << pcpg_iter->getNumIters() << std::endl
898 << e.what() << std::endl;
904 Teuchos::RCP<MV> update = pcpg_iter->getCurrentUpdate();
905 problem_->updateSolution( update,
true );
908 problem_->setCurrLS();
916 std::cout <<
"SolverManager: dimU_ " << dimU_ <<
" prevUdim= " << q << std::endl;
918 if( q > deflatedBlocks_ )
919 std::cout <<
"SolverManager: Error deflatedBlocks = " << deflatedBlocks_ << std::endl;
930 rank = ARRQR(dimU_,q, *oldState.
D );
932 std::cout <<
" rank decreased in ARRQR, something to do? " << std::endl;
938 if( dimU_ > deflatedBlocks_ ){
940 if( !deflatedBlocks_ ){
943 dimU_ = deflatedBlocks_;
947 bool Harmonic =
false;
949 Teuchos::RCP<MV> Uorth;
951 std::vector<int> active_cols( dimU_ );
952 for (
int i=0; i < dimU_; ++i) active_cols[i] = i;
955 Uorth = MVT::CloneCopy(*C_, active_cols);
958 Uorth = MVT::CloneCopy(*U_, active_cols);
962 Teuchos::SerialDenseMatrix<int,ScalarType> R(dimU_,dimU_);
963 rank = ortho_->normalize(*Uorth, Teuchos::rcp(&R,
false));
964 Uorth = Teuchos::null;
970 "Belos::PCPGSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
974 Teuchos::SerialDenseMatrix<int,ScalarType> VT;
975 Teuchos::SerialDenseMatrix<int,ScalarType> Ur;
979 if( problem_->isHermitian() ) lrwork = dimU_;
980 std::vector<ScalarType> work(lwork);
981 std::vector<ScalarType> Svec(dimU_);
982 std::vector<ScalarType> rwork(lrwork);
983 lapack.GESVD(
'N',
'O',
984 R.numRows(),R.numCols(),R.values(), R.numRows(),
992 "Belos::PCPGSolMgr::solve(): LAPACK _GESVD failed to compute singular values.");
994 if( work[0] != 67. * dimU_ )
995 std::cout <<
" SVD " << dimU_ <<
" lwork " << work[0] << std::endl;
996 for(
int i=0; i< dimU_; i++)
997 std::cout << i <<
" " << Svec[i] << std::endl;
999 Teuchos::SerialDenseMatrix<int,ScalarType> wholeV( R, Teuchos::TRANS);
1001 int startRow = 0, startCol = 0;
1003 startCol = dimU_ - deflatedBlocks_;
1005 Teuchos::SerialDenseMatrix<int,ScalarType> V(Teuchos::Copy,
1011 std::vector<int> active_columns( dimU_ );
1012 std::vector<int> def_cols( deflatedBlocks_ );
1013 for (
int i=0; i < dimU_; ++i) active_columns[i] = i;
1014 for (
int i=0; i < deflatedBlocks_; ++i) def_cols[i] = i;
1016 Teuchos::RCP<MV> Uactive = MVT::CloneViewNonConst(*U_, def_cols);
1017 Teuchos::RCP<MV> Ucopy = MVT::CloneCopy( *U_, active_columns );
1018 MVT::MvTimesMatAddMv( one, *Ucopy, V, zero, *Uactive );
1019 Ucopy = Teuchos::null;
1020 Uactive = Teuchos::null;
1021 Teuchos::RCP<MV> Cactive = MVT::CloneViewNonConst(*C_, def_cols);
1022 Teuchos::RCP<MV> Ccopy = MVT::CloneCopy( *C_, active_columns );
1023 MVT::MvTimesMatAddMv( one, *Ccopy, V, zero, *Cactive );
1024 Ccopy = Teuchos::null;
1025 Cactive = Teuchos::null;
1026 dimU_ = deflatedBlocks_;
1028 printer_->stream(
Debug) <<
" Generated recycled subspace using RHS index " << currIdx[0] <<
" of dimension " << dimU_ << std::endl << std::endl;
1031 problem_->setCurrLS();
1035 if ( numRHS2Solve > 0 ) {
1039 problem_->setLSIndex( currIdx );
1042 currIdx.resize( numRHS2Solve );
1051#ifdef BELOS_TEUCHOS_TIME_MONITOR
1056 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails) );
1061 using Teuchos::rcp_dynamic_cast;
1064 const std::vector<MagnitudeType>* pTestValues =
1065 rcp_dynamic_cast<conv_test_type>(convTest_)->
getTestValue();
1067 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1068 "Belos::PCPGSolMgr::solve(): The convergence test's getTestValue() "
1069 "method returned NULL. Please report this bug to the Belos developers.");
1071 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1072 "Belos::PCPGSolMgr::solve(): The convergence test's getTestValue() "
1073 "method returned a vector of length zero. Please report this bug to the "
1074 "Belos developers.");
1079 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1083 numIters_ = maxIterTest_->getNumIters();