42 #ifndef BELOS_PCPG_SOLMGR_HPP 43 #define BELOS_PCPG_SOLMGR_HPP 64 #include "Teuchos_BLAS.hpp" 65 #include "Teuchos_LAPACK.hpp" 66 #ifdef BELOS_TEUCHOS_TIME_MONITOR 67 # include "Teuchos_TimeMonitor.hpp" 69 #if defined(HAVE_TEUCHOSCORE_CXX11) 70 # include <type_traits> 71 #endif // defined(HAVE_TEUCHOSCORE_CXX11) 72 #include "Teuchos_TypeTraits.hpp" 151 template<
class ScalarType,
class MV,
class OP,
152 const bool supportsScalarType =
154 ! Teuchos::ScalarTraits<ScalarType>::isComplex>
157 Belos::Details::LapackSupportsScalar<ScalarType>::value &&
158 ! Teuchos::ScalarTraits<ScalarType>::isComplex>
160 static const bool scalarTypeIsSupported =
162 ! Teuchos::ScalarTraits<ScalarType>::isComplex;
171 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
177 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
182 template<
class ScalarType,
class MV,
class OP>
188 typedef Teuchos::ScalarTraits<ScalarType> SCT;
189 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
190 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
240 const Teuchos::RCP<Teuchos::ParameterList> &pl );
246 virtual Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const {
262 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters()
const;
273 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
274 return Teuchos::tuple(timerSolve_);
304 void setParameters(
const Teuchos::RCP<Teuchos::ParameterList> ¶ms );
345 std::string description()
const;
353 int ARRQR(
int numVecs,
int numOrthVecs,
const Teuchos::SerialDenseMatrix<int,ScalarType>& D);
356 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
359 Teuchos::RCP<OutputManager<ScalarType> > printer_;
360 Teuchos::RCP<std::ostream> outputStream_;
363 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
364 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
365 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
366 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
369 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
372 Teuchos::RCP<Teuchos::ParameterList> params_;
375 static constexpr
int maxIters_default_ = 1000;
376 static constexpr
int deflatedBlocks_default_ = 2;
377 static constexpr
int savedBlocks_default_ = 16;
380 static constexpr
int outputFreq_default_ = -1;
381 static constexpr
const char * label_default_ =
"Belos";
382 static constexpr
const char * orthoType_default_ =
"DGKS";
383 static constexpr std::ostream * outputStream_default_ = &std::cout;
390 MagnitudeType convtol_;
393 MagnitudeType orthoKappa_;
396 MagnitudeType achievedTol_;
404 int deflatedBlocks_, savedBlocks_, verbosity_, outputStyle_, outputFreq_;
405 std::string orthoType_;
408 Teuchos::RCP<MV> U_, C_, R_;
415 Teuchos::RCP<Teuchos::Time> timerSolve_;
423 template<
class ScalarType,
class MV,
class OP>
425 outputStream_(Teuchos::rcp(outputStream_default_,false)),
428 achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
430 maxIters_(maxIters_default_),
431 deflatedBlocks_(deflatedBlocks_default_),
432 savedBlocks_(savedBlocks_default_),
433 verbosity_(verbosity_default_),
434 outputStyle_(outputStyle_default_),
435 outputFreq_(outputFreq_default_),
436 orthoType_(orthoType_default_),
438 label_(label_default_),
444 template<
class ScalarType,
class MV,
class OP>
447 const Teuchos::RCP<Teuchos::ParameterList> &pl ) :
449 outputStream_(Teuchos::rcp(outputStream_default_,false)),
453 achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
455 maxIters_(maxIters_default_),
456 deflatedBlocks_(deflatedBlocks_default_),
457 savedBlocks_(savedBlocks_default_),
458 verbosity_(verbosity_default_),
459 outputStyle_(outputStyle_default_),
460 outputFreq_(outputFreq_default_),
461 orthoType_(orthoType_default_),
463 label_(label_default_),
466 TEUCHOS_TEST_FOR_EXCEPTION(
467 problem_.is_null (), std::invalid_argument,
468 "Belos::PCPGSolMgr two-argument constructor: " 469 "'problem' is null. You must supply a non-null Belos::LinearProblem " 470 "instance when calling this constructor.");
472 if (! pl.is_null ()) {
479 template<
class ScalarType,
class MV,
class OP>
483 if (params_ == Teuchos::null) {
484 params_ = Teuchos::rcp(
new Teuchos::ParameterList(*getValidParameters()) );
487 params->validateParameters(*getValidParameters());
491 if (params->isParameter(
"Maximum Iterations")) {
492 maxIters_ = params->get(
"Maximum Iterations",maxIters_default_);
495 params_->set(
"Maximum Iterations", maxIters_);
496 if (maxIterTest_!=Teuchos::null)
497 maxIterTest_->setMaxIters( maxIters_ );
501 if (params->isParameter(
"Num Saved Blocks")) {
502 savedBlocks_ = params->get(
"Num Saved Blocks",savedBlocks_default_);
503 TEUCHOS_TEST_FOR_EXCEPTION(savedBlocks_ <= 0, std::invalid_argument,
504 "Belos::PCPGSolMgr: \"Num Saved Blocks\" must be strictly positive.");
511 params_->set(
"Num Saved Blocks", savedBlocks_);
513 if (params->isParameter(
"Num Deflated Blocks")) {
514 deflatedBlocks_ = params->get(
"Num Deflated Blocks",deflatedBlocks_default_);
515 TEUCHOS_TEST_FOR_EXCEPTION(deflatedBlocks_ < 0, std::invalid_argument,
516 "Belos::PCPGSolMgr: \"Num Deflated Blocks\" must be positive.");
518 TEUCHOS_TEST_FOR_EXCEPTION(deflatedBlocks_ > savedBlocks_, std::invalid_argument,
519 "Belos::PCPGSolMgr: \"Num Deflated Blocks\" must be <= \"Num Saved Blocks\".");
523 params_->set(
"Num Deflated Blocks", static_cast<int>(deflatedBlocks_));
527 if (params->isParameter(
"Timer Label")) {
528 std::string tempLabel = params->get(
"Timer Label", label_default_);
531 if (tempLabel != label_) {
533 params_->set(
"Timer Label", label_);
534 std::string solveLabel = label_ +
": PCPGSolMgr total solve time";
535 #ifdef BELOS_TEUCHOS_TIME_MONITOR 536 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
538 if (ortho_ != Teuchos::null) {
539 ortho_->setLabel( label_ );
545 if (params->isParameter(
"Orthogonalization")) {
546 std::string tempOrthoType = params->get(
"Orthogonalization",orthoType_default_);
547 TEUCHOS_TEST_FOR_EXCEPTION( tempOrthoType !=
"DGKS" && tempOrthoType !=
"ICGS" && tempOrthoType !=
"IMGS",
548 std::invalid_argument,
549 "Belos::PCPGSolMgr: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
550 if (tempOrthoType != orthoType_) {
551 orthoType_ = tempOrthoType;
552 params_->set(
"Orthogonalization", orthoType_);
554 if (orthoType_==
"DGKS") {
555 if (orthoKappa_ <= 0) {
563 else if (orthoType_==
"ICGS") {
566 else if (orthoType_==
"IMGS") {
573 if (params->isParameter(
"Orthogonalization Constant")) {
574 if (params->isType<MagnitudeType> (
"Orthogonalization Constant")) {
575 orthoKappa_ = params->get (
"Orthogonalization Constant",
579 orthoKappa_ = params->get (
"Orthogonalization Constant",
584 params_->set(
"Orthogonalization Constant",orthoKappa_);
585 if (orthoType_==
"DGKS") {
586 if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
593 if (params->isParameter(
"Verbosity")) {
594 if (Teuchos::isParameterType<int>(*params,
"Verbosity")) {
595 verbosity_ = params->get(
"Verbosity", verbosity_default_);
597 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,
"Verbosity");
601 params_->set(
"Verbosity", verbosity_);
602 if (printer_ != Teuchos::null)
603 printer_->setVerbosity(verbosity_);
607 if (params->isParameter(
"Output Style")) {
608 if (Teuchos::isParameterType<int>(*params,
"Output Style")) {
609 outputStyle_ = params->get(
"Output Style", outputStyle_default_);
611 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,
"Output Style");
615 params_->set(
"Output Style", outputStyle_);
616 outputTest_ = Teuchos::null;
620 if (params->isParameter(
"Output Stream")) {
621 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,
"Output Stream");
624 params_->set(
"Output Stream", outputStream_);
625 if (printer_ != Teuchos::null)
626 printer_->setOStream( outputStream_ );
631 if (params->isParameter(
"Output Frequency")) {
632 outputFreq_ = params->get(
"Output Frequency", outputFreq_default_);
636 params_->set(
"Output Frequency", outputFreq_);
637 if (outputTest_ != Teuchos::null)
638 outputTest_->setOutputFrequency( outputFreq_ );
642 if (printer_ == Teuchos::null) {
651 if (params->isParameter(
"Convergence Tolerance")) {
652 if (params->isType<MagnitudeType> (
"Convergence Tolerance")) {
653 convtol_ = params->get (
"Convergence Tolerance",
661 params_->set(
"Convergence Tolerance", convtol_);
662 if (convTest_ != Teuchos::null)
669 if (maxIterTest_ == Teuchos::null)
672 if (convTest_ == Teuchos::null)
673 convTest_ = Teuchos::rcp(
new StatusTestResNorm_t( convtol_, 1 ) );
675 sTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
683 std::string solverDesc =
" PCPG ";
684 outputTest_->setSolverDesc( solverDesc );
688 if (ortho_ == Teuchos::null) {
689 params_->set(
"Orthogonalization", orthoType_);
690 if (orthoType_==
"DGKS") {
691 if (orthoKappa_ <= 0) {
699 else if (orthoType_==
"ICGS") {
702 else if (orthoType_==
"IMGS") {
706 TEUCHOS_TEST_FOR_EXCEPTION(orthoType_!=
"ICGS"&&orthoType_!=
"DGKS"&&orthoType_!=
"IMGS",std::logic_error,
707 "Belos::PCPGSolMgr(): Invalid orthogonalization type.");
712 if (timerSolve_ == Teuchos::null) {
713 std::string solveLabel = label_ +
": PCPGSolMgr total solve time";
714 #ifdef BELOS_TEUCHOS_TIME_MONITOR 715 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
724 template<
class ScalarType,
class MV,
class OP>
725 Teuchos::RCP<const Teuchos::ParameterList>
728 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
729 if (is_null(validPL)) {
730 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
733 "The relative residual tolerance that needs to be achieved by the\n" 734 "iterative solver in order for the linear system to be declared converged.");
735 pl->set(
"Maximum Iterations", static_cast<int>(maxIters_default_),
736 "The maximum number of iterations allowed for each\n" 737 "set of RHS solved.");
738 pl->set(
"Num Deflated Blocks", static_cast<int>(deflatedBlocks_default_),
739 "The maximum number of vectors in the seed subspace." );
740 pl->set(
"Num Saved Blocks", static_cast<int>(savedBlocks_default_),
741 "The maximum number of vectors saved from old Krylov subspaces." );
742 pl->set(
"Verbosity", static_cast<int>(verbosity_default_),
743 "What type(s) of solver information should be outputted\n" 744 "to the output stream.");
745 pl->set(
"Output Style", static_cast<int>(outputStyle_default_),
746 "What style is used for the solver information outputted\n" 747 "to the output stream.");
748 pl->set(
"Output Frequency", static_cast<int>(outputFreq_default_),
749 "How often convergence information should be outputted\n" 750 "to the output stream.");
751 pl->set(
"Output Stream", Teuchos::rcp(outputStream_default_,
false),
752 "A reference-counted pointer to the output stream where all\n" 753 "solver output is sent.");
754 pl->set(
"Timer Label", static_cast<const char *>(label_default_),
755 "The string to use as a prefix for the timer labels.");
756 pl->set(
"Orthogonalization", static_cast<const char *>(orthoType_default_),
757 "The type of orthogonalization to use: DGKS, ICGS, IMGS");
759 "The constant used by DGKS orthogonalization to determine\n" 760 "whether another step of classical Gram-Schmidt is necessary.");
768 template<
class ScalarType,
class MV,
class OP>
772 if (!isSet_) { setParameters( params_ ); }
774 Teuchos::BLAS<int,ScalarType> blas;
775 Teuchos::LAPACK<int,ScalarType> lapack;
776 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
777 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
780 "Belos::PCPGSolMgr::solve(): Linear problem is not a valid object.");
783 "Belos::PCPGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
786 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
787 std::vector<int> currIdx(1);
793 problem_->setLSIndex( currIdx );
796 bool isConverged =
true;
800 Teuchos::ParameterList plist;
801 plist.set(
"Saved Blocks", savedBlocks_);
802 plist.set(
"Block Size", 1);
803 plist.set(
"Keep Diagonal",
true);
804 plist.set(
"Initialize Diagonal",
true);
809 Teuchos::RCP<PCPGIter<ScalarType,MV,OP> > pcpg_iter;
815 #ifdef BELOS_TEUCHOS_TIME_MONITOR 816 Teuchos::TimeMonitor slvtimer(*timerSolve_);
818 while ( numRHS2Solve > 0 ) {
821 outputTest_->reset();
824 if (R_ == Teuchos::null)
825 R_ = MVT::Clone( *(problem_->getRHS()), 1 );
827 problem_->computeCurrResVec( &*R_ );
833 if( U_ != Teuchos::null ){
838 Teuchos::RCP<MV> cur_soln_vec = problem_->getCurrLHSVec();
839 std::vector<MagnitudeType> rnorm0(1);
840 MVT::MvNorm( *R_, rnorm0 );
843 std::cout <<
"Solver Manager: dimU_ = " << dimU_ << std::endl;
844 Teuchos::SerialDenseMatrix<int,ScalarType> Z( dimU_, 1 );
846 Teuchos::RCP<const MV> Uactive, Cactive;
847 std::vector<int> active_columns( dimU_ );
848 for (
int i=0; i < dimU_; ++i) active_columns[i] = i;
849 Uactive = MVT::CloneView(*U_, active_columns);
850 Cactive = MVT::CloneView(*C_, active_columns);
853 std::cout <<
" Solver Manager : check duality of seed basis " << std::endl;
854 Teuchos::SerialDenseMatrix<int,ScalarType> H( dimU_, dimU_ );
855 MVT::MvTransMv( one, *Uactive, *Cactive, H );
856 H.print( std::cout );
859 MVT::MvTransMv( one, *Uactive, *R_, Z );
860 Teuchos::RCP<MV> tempU = MVT::Clone( *R_, 1 );
861 MVT::MvTimesMatAddMv( one, *Uactive, Z, zero, *tempU );
862 MVT::MvAddMv( one, *tempU, one, *cur_soln_vec, *cur_soln_vec );
863 MVT::MvTimesMatAddMv( one, *Cactive, Z, zero, *tempU );
864 MVT::MvAddMv( -one, *tempU, one, *R_, *R_ );
865 std::vector<MagnitudeType> rnorm(1);
866 MVT::MvNorm( *R_, rnorm );
867 if( rnorm[0] < rnorm0[0] * .001 ){
868 MVT::MvTransMv( one, *Uactive, *R_, Z );
869 MVT::MvTimesMatAddMv( one, *Uactive, Z, zero, *tempU );
870 MVT::MvAddMv( one, *tempU, one, *cur_soln_vec, *cur_soln_vec );
871 MVT::MvTimesMatAddMv( one, *Cactive, Z, zero, *tempU );
872 MVT::MvAddMv( -one, *tempU, one, *R_, *R_ );
874 Uactive = Teuchos::null;
875 Cactive = Teuchos::null;
876 tempU = Teuchos::null;
887 if( U_ != Teuchos::null ) pcpgState.
U = U_;
888 if( C_ != Teuchos::null ) pcpgState.
C = C_;
889 if( dimU_ > 0 ) pcpgState.
curDim = dimU_;
890 pcpg_iter->initialize(pcpgState);
896 if( !dimU_ ) printer_->stream(
Debug) <<
" No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
897 pcpg_iter->resetNumIters();
899 if( dimU_ > savedBlocks_ )
900 std::cout <<
"Error: dimU_ = " << dimU_ <<
" > savedBlocks_ = " << savedBlocks_ << std::endl;
906 if( debug ) printf(
"********** Calling iterate...\n");
907 pcpg_iter->iterate();
914 if ( convTest_->getStatus() ==
Passed ) {
923 else if ( maxIterTest_->getStatus() ==
Passed ) {
937 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
938 "Belos::PCPGSolMgr::solve(): Invalid return from PCPGIter::iterate().");
944 sTest_->checkStatus( &*pcpg_iter );
945 if (convTest_->getStatus() !=
Passed)
949 catch (
const std::exception &e) {
950 printer_->stream(
Errors) <<
"Error! Caught exception in PCPGIter::iterate() at iteration " 951 << pcpg_iter->getNumIters() << std::endl
952 << e.what() << std::endl;
958 Teuchos::RCP<MV> update = pcpg_iter->getCurrentUpdate();
959 problem_->updateSolution( update,
true );
962 problem_->setCurrLS();
970 std::cout <<
"SolverManager: dimU_ " << dimU_ <<
" prevUdim= " << q << std::endl;
972 if( q > deflatedBlocks_ )
973 std::cout <<
"SolverManager: Error deflatedBlocks = " << deflatedBlocks_ << std::endl;
984 rank = ARRQR(dimU_,q, *oldState.
D );
986 std::cout <<
" rank decreased in ARRQR, something to do? " << std::endl;
992 if( dimU_ > deflatedBlocks_ ){
994 if( !deflatedBlocks_ ){
997 dimU_ = deflatedBlocks_;
1001 bool Harmonic =
false;
1003 Teuchos::RCP<MV> Uorth;
1005 std::vector<int> active_cols( dimU_ );
1006 for (
int i=0; i < dimU_; ++i) active_cols[i] = i;
1009 Uorth = MVT::CloneCopy(*C_, active_cols);
1012 Uorth = MVT::CloneCopy(*U_, active_cols);
1016 Teuchos::SerialDenseMatrix<int,ScalarType> R(dimU_,dimU_);
1017 rank = ortho_->normalize(*Uorth, Teuchos::rcp(&R,
false));
1018 Uorth = Teuchos::null;
1024 "Belos::PCPGSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
1028 Teuchos::SerialDenseMatrix<int,ScalarType> VT;
1029 Teuchos::SerialDenseMatrix<int,ScalarType> Ur;
1030 int lwork = 5*dimU_;
1033 if( problem_->isHermitian() ) lrwork = dimU_;
1034 std::vector<ScalarType> work(lwork);
1035 std::vector<ScalarType> Svec(dimU_);
1036 std::vector<ScalarType> rwork(lrwork);
1037 lapack.GESVD(
'N',
'O',
1038 R.numRows(),R.numCols(),R.values(), R.numRows(),
1046 "Belos::PCPGSolMgr::solve(): LAPACK _GESVD failed to compute singular values.");
1048 if( work[0] != 67. * dimU_ )
1049 std::cout <<
" SVD " << dimU_ <<
" lwork " << work[0] << std::endl;
1050 for(
int i=0; i< dimU_; i++)
1051 std::cout << i <<
" " << Svec[i] << std::endl;
1053 Teuchos::SerialDenseMatrix<int,ScalarType> wholeV( R,
Teuchos::TRANS);
1055 int startRow = 0, startCol = 0;
1057 startCol = dimU_ - deflatedBlocks_;
1059 Teuchos::SerialDenseMatrix<int,ScalarType> V(Teuchos::Copy,
1065 std::vector<int> active_columns( dimU_ );
1066 std::vector<int> def_cols( deflatedBlocks_ );
1067 for (
int i=0; i < dimU_; ++i) active_columns[i] = i;
1068 for (
int i=0; i < deflatedBlocks_; ++i) def_cols[i] = i;
1070 Teuchos::RCP<MV> Uactive = MVT::CloneViewNonConst(*U_, def_cols);
1071 Teuchos::RCP<MV> Ucopy = MVT::CloneCopy( *U_, active_columns );
1072 MVT::MvTimesMatAddMv( one, *Ucopy, V, zero, *Uactive );
1073 Ucopy = Teuchos::null;
1074 Uactive = Teuchos::null;
1075 Teuchos::RCP<MV> Cactive = MVT::CloneViewNonConst(*C_, def_cols);
1076 Teuchos::RCP<MV> Ccopy = MVT::CloneCopy( *C_, active_columns );
1077 MVT::MvTimesMatAddMv( one, *Ccopy, V, zero, *Cactive );
1078 Ccopy = Teuchos::null;
1079 Cactive = Teuchos::null;
1080 dimU_ = deflatedBlocks_;
1082 printer_->stream(
Debug) <<
" Generated recycled subspace using RHS index " << currIdx[0] <<
" of dimension " << dimU_ << std::endl << std::endl;
1085 problem_->setCurrLS();
1089 if ( numRHS2Solve > 0 ) {
1093 problem_->setLSIndex( currIdx );
1096 currIdx.resize( numRHS2Solve );
1105 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1110 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails) );
1115 using Teuchos::rcp_dynamic_cast;
1118 const std::vector<MagnitudeType>* pTestValues =
1119 rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
1121 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1122 "Belos::PCPGSolMgr::solve(): The convergence test's getTestValue() " 1123 "method returned NULL. Please report this bug to the Belos developers.");
1125 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1126 "Belos::PCPGSolMgr::solve(): The convergence test's getTestValue() " 1127 "method returned a vector of length zero. Please report this bug to the " 1128 "Belos developers.");
1133 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1137 numIters_ = maxIterTest_->getNumIters();
1148 template<
class ScalarType,
class MV,
class OP>
1152 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1153 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1156 Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
1157 Teuchos::SerialDenseMatrix<int,ScalarType> gamma( 1, 1 );
1158 Teuchos::SerialDenseMatrix<int,ScalarType> anorm( 1, 1 );
1159 std::vector<int> curind(1);
1160 std::vector<int> ipiv(p - q);
1161 std::vector<ScalarType> Pivots(p);
1163 ScalarType rteps = 1.5e-8;
1166 for( i = q ; i < p ; i++ ){
1169 RCP<MV> P = MVT::CloneViewNonConst(*U_,curind);
1170 RCP<MV> AP = MVT::CloneViewNonConst(*C_,curind);
1171 anorm(0,0) = one / Teuchos::ScalarTraits<ScalarType>::squareroot( D(i-q,i-q) ) ;
1172 MVT::MvAddMv( anorm(0,0), *P, zero, *AP, *P );
1173 MVT::MvAddMv( zero, *P, anorm(0,0), *AP, *AP );
1177 for( i = q ; i < p ; i++ ){
1178 if( q < i && i < p-1 ){
1181 for( j = i+1 ; j < p ; j++ ){
1182 const int k = ipiv[j-q];
1183 if( Pivots[k] > Pivots[l] ){
1190 ipiv[imax-q] = ipiv[i-q];
1196 if( Pivots[k] > 1.5625e-2 ){
1197 anorm(0,0) = Pivots[k];
1201 RCP<const MV> P = MVT::CloneView(*U_,curind);
1202 RCP<const MV> AP = MVT::CloneView(*C_,curind);
1203 MVT::MvTransMv( one, *P, *AP, anorm );
1204 anorm(0,0) = Teuchos::ScalarTraits<ScalarType>::squareroot( anorm(0,0) ) ;
1206 if( rteps <= anorm(0,0) && anorm(0,0) < 9.765625e-4){
1214 std::cout <<
"ARRQR: Bad case not implemented" << std::endl;
1216 if( anorm(0,0) < rteps ){
1217 std::cout <<
"ARRQR : deficient case not implemented " << std::endl;
1225 RCP<MV> P = MVT::CloneViewNonConst(*U_,curind);
1226 RCP<MV> AP = MVT::CloneViewNonConst(*C_,curind);
1227 MVT::MvAddMv( anorm(0,0), *P, zero, *AP, *P );
1228 MVT::MvAddMv( zero, *P, anorm(0,0), *AP, *AP );
1232 P = MVT::CloneViewNonConst(*U_,curind);
1233 AP = MVT::CloneViewNonConst(*C_,curind);
1234 for( j = i+1 ; j < p ; j++ ){
1237 RCP<MV> Q = MVT::CloneViewNonConst(*U_,curind);
1238 MVT::MvTransMv( one, *Q, *AP, alpha);
1239 MVT::MvAddMv( -alpha(0,0), *P, one, *Q, *Q );
1241 RCP<MV> AQ = MVT::CloneViewNonConst(*C_,curind);
1242 MVT::MvAddMv( -alpha(0,0), *AP, one, *AQ, *AQ );
1244 gamma(0,0) = ( Pivots[l] - alpha(0,0))*( Pivots[l] + alpha(0,0));
1245 if( gamma(0,0) > 0){
1246 Pivots[l] = Teuchos::ScalarTraits<ScalarType>::squareroot( gamma(0,0) );
1257 template<
class ScalarType,
class MV,
class OP>
1260 std::ostringstream oss;
1261 oss <<
"Belos::PCPGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<
">";
1263 oss <<
"Ortho Type='"<<orthoType_;
bool isLOADetected() const
Return whether a loss of accuracy was detected by this solver during the most current solve...
Collection of types and exceptions used within the Belos solvers.
int prevUdim
Number of block columns in matrices C and U before current iteration.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
PCPGSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
Class which manages the output and verbosity of the Belos solvers.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
Belos concrete class to iterate Preconditioned Conjugate Projected Gradients.
static constexpr double orthoKappa
DGKS orthogonalization constant.
PCPGSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate orthonormal...
Teuchos::RCP< MV > R
The current residual.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get current linear problem being solved for in this object.
Teuchos::RCP< MV > C
C = AU, U spans recycled subspace.
Base class for Belos::SolverManager subclasses which normally can only compile with real ScalarType t...
PCPGSolMgrRecyclingFailure is thrown when any problem occurs in using/creating the recycling subspace...
A factory class for generating StatusTestOutput objects.
int curDim
The current dimension of the reduction.
An implementation of StatusTestResNorm using a family of residual norms.
PCPGSolMgrOrthoFailure(const std::string &what_arg)
Structure to contain pointers to PCPGIter state variables.
PCPGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i.e.
Belos::StatusTest class for specifying a maximum number of iterations.
A factory class for generating StatusTestOutput objects.
int getNumIters() const
Get the iteration count for the most recent call to solve().
Iterated Modified Gram-Schmidt (IMGS) implementation of the Belos::OrthoManager class.
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > D
The current Hessenberg matrix.
Classical Gram-Schmidt (with DGKS correction) implementation of the Belos::OrthoManager class...
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
Pure virtual base class which describes the basic interface for a solver manager. ...
static constexpr double convTol
Default convergence tolerance.
void reset(const ResetType type)
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
virtual Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const
clone for Inverted Injection (DII)
A linear system to solve, and its associated information.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem)
Set the linear problem that needs to be solved.
Class which describes the linear problem to be solved by the iterative solver.
PCPGSolMgrLAPACKFailure is thrown when a nonzero value is retuned from an LAPACK call.
PCPGSolMgrLAPACKFailure(const std::string &what_arg)
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
ReturnType
Whether the Belos solve converged for all linear systems.
MagnitudeType achievedTol() const
Tolerance achieved by the last solve() invocation.
Iterated Classical Gram-Schmidt (ICGS) implementation of the Belos::OrthoManager class.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
Teuchos::RCP< MV > U
The recycled subspace.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
int setTolerance(MagnitudeType tolerance)
Set the value of the tolerance.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > create(const Teuchos::RCP< OutputManager< ScalarType > > &printer, Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test, int mod, int printStates)
Create the StatusTestOutput object specified by the outputStyle.
PCPGIterOrthoFailure is thrown when the PCPGIter object is unable to compute independent direction ve...
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
PCPG iterative linear solver.
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
A class for extending the status testing capabilities of Belos via logical combinations.
virtual ~PCPGSolMgr()
Destructor.
Class which defines basic traits for the operator type.
Parent class to all Belos exceptions.
Default parameters common to most Belos solvers.
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
PCPGSolMgrRecyclingFailure(const std::string &what_arg)
This class implements the PCPG iteration, where a single-std::vector Krylov subspace is constructed...
Belos header file which uses auto-configuration information to include necessary C++ headers...
PCPGSolMgrLinearProblemFailure(const std::string &what_arg)