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>
162 ! Teuchos::ScalarTraits<ScalarType>::isComplex;
171 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
177 template<
class ScalarType,
class MV,
class OP>
183 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
184 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
MagnitudeType;
185 typedef Teuchos::ScalarTraits<MagnitudeType>
MT;
235 const Teuchos::RCP<Teuchos::ParameterList> &pl );
252 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters()
const;
263 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
264 return Teuchos::tuple(timerSolve_);
294 void setParameters(
const Teuchos::RCP<Teuchos::ParameterList> ¶ms );
335 std::string description()
const;
343 int ARRQR(
int numVecs,
int numOrthVecs,
const Teuchos::SerialDenseMatrix<int,ScalarType>& D);
346 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >
problem_;
353 Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
sTest_;
355 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> >
convTest_;
359 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> >
ortho_;
396 int deflatedBlocks_, savedBlocks_,
verbosity_, outputStyle_, outputFreq_;
400 Teuchos::RCP<MV>
U_, C_, R_;
415 template<
class ScalarType,
class MV,
class OP>
419 template<
class ScalarType,
class MV,
class OP>
423 template<
class ScalarType,
class MV,
class OP>
426 template<
class ScalarType,
class MV,
class OP>
429 template<
class ScalarType,
class MV,
class OP>
432 template<
class ScalarType,
class MV,
class OP>
435 template<
class ScalarType,
class MV,
class OP>
438 template<
class ScalarType,
class MV,
class OP>
441 template<
class ScalarType,
class MV,
class OP>
444 template<
class ScalarType,
class MV,
class OP>
447 template<
class ScalarType,
class MV,
class OP>
452 template<
class ScalarType,
class MV,
class OP>
454 outputStream_(outputStream_default_),
455 convtol_(convtol_default_),
456 orthoKappa_(orthoKappa_default_),
459 maxIters_(maxIters_default_),
460 deflatedBlocks_(deflatedBlocks_default_),
461 savedBlocks_(savedBlocks_default_),
462 verbosity_(verbosity_default_),
463 outputStyle_(outputStyle_default_),
464 outputFreq_(outputFreq_default_),
465 orthoType_(orthoType_default_),
467 label_(label_default_),
473 template<
class ScalarType,
class MV,
class OP>
476 const Teuchos::RCP<Teuchos::ParameterList> &pl ) :
478 outputStream_(outputStream_default_),
479 convtol_(convtol_default_),
480 orthoKappa_(orthoKappa_default_),
483 maxIters_(maxIters_default_),
484 deflatedBlocks_(deflatedBlocks_default_),
485 savedBlocks_(savedBlocks_default_),
486 verbosity_(verbosity_default_),
487 outputStyle_(outputStyle_default_),
488 outputFreq_(outputFreq_default_),
489 orthoType_(orthoType_default_),
491 label_(label_default_),
494 TEUCHOS_TEST_FOR_EXCEPTION(
495 problem_.is_null (), std::invalid_argument,
496 "Belos::PCPGSolMgr two-argument constructor: " 497 "'problem' is null. You must supply a non-null Belos::LinearProblem " 498 "instance when calling this constructor.");
500 if (! pl.is_null ()) {
507 template<
class ScalarType,
class MV,
class OP>
511 if (params_ == Teuchos::null) {
512 params_ = Teuchos::rcp(
new Teuchos::ParameterList(*getValidParameters()) );
515 params->validateParameters(*getValidParameters());
519 if (params->isParameter(
"Maximum Iterations")) {
520 maxIters_ = params->get(
"Maximum Iterations",maxIters_default_);
523 params_->set(
"Maximum Iterations", maxIters_);
524 if (maxIterTest_!=Teuchos::null)
525 maxIterTest_->setMaxIters( maxIters_ );
529 if (params->isParameter(
"Num Saved Blocks")) {
530 savedBlocks_ = params->get(
"Num Saved Blocks",savedBlocks_default_);
531 TEUCHOS_TEST_FOR_EXCEPTION(savedBlocks_ <= 0, std::invalid_argument,
532 "Belos::PCPGSolMgr: \"Num Saved Blocks\" must be strictly positive.");
539 params_->set(
"Num Saved Blocks", savedBlocks_);
541 if (params->isParameter(
"Num Deflated Blocks")) {
542 deflatedBlocks_ = params->get(
"Num Deflated Blocks",deflatedBlocks_default_);
543 TEUCHOS_TEST_FOR_EXCEPTION(deflatedBlocks_ < 0, std::invalid_argument,
544 "Belos::PCPGSolMgr: \"Num Deflated Blocks\" must be positive.");
546 TEUCHOS_TEST_FOR_EXCEPTION(deflatedBlocks_ > savedBlocks_, std::invalid_argument,
547 "Belos::PCPGSolMgr: \"Num Deflated Blocks\" must be <= \"Num Saved Blocks\".");
550 params_->set(
"Num Deflated Blocks", deflatedBlocks_);
554 if (params->isParameter(
"Timer Label")) {
555 std::string tempLabel = params->get(
"Timer Label", label_default_);
558 if (tempLabel != label_) {
560 params_->set(
"Timer Label", label_);
561 std::string solveLabel = label_ +
": PCPGSolMgr total solve time";
562 #ifdef BELOS_TEUCHOS_TIME_MONITOR 563 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
565 if (ortho_ != Teuchos::null) {
566 ortho_->setLabel( label_ );
572 if (params->isParameter(
"Orthogonalization")) {
573 std::string tempOrthoType = params->get(
"Orthogonalization",orthoType_default_);
574 TEUCHOS_TEST_FOR_EXCEPTION( tempOrthoType !=
"DGKS" && tempOrthoType !=
"ICGS" && tempOrthoType !=
"IMGS",
575 std::invalid_argument,
576 "Belos::PCPGSolMgr: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
577 if (tempOrthoType != orthoType_) {
578 orthoType_ = tempOrthoType;
580 if (orthoType_==
"DGKS") {
581 if (orthoKappa_ <= 0) {
589 else if (orthoType_==
"ICGS") {
592 else if (orthoType_==
"IMGS") {
599 if (params->isParameter(
"Orthogonalization Constant")) {
600 orthoKappa_ = params->get(
"Orthogonalization Constant",orthoKappa_default_);
603 params_->set(
"Orthogonalization Constant",orthoKappa_);
604 if (orthoType_==
"DGKS") {
605 if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
612 if (params->isParameter(
"Verbosity")) {
613 if (Teuchos::isParameterType<int>(*params,
"Verbosity")) {
614 verbosity_ = params->get(
"Verbosity", verbosity_default_);
616 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,
"Verbosity");
620 params_->set(
"Verbosity", verbosity_);
621 if (printer_ != Teuchos::null)
622 printer_->setVerbosity(verbosity_);
626 if (params->isParameter(
"Output Style")) {
627 if (Teuchos::isParameterType<int>(*params,
"Output Style")) {
628 outputStyle_ = params->get(
"Output Style", outputStyle_default_);
630 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,
"Output Style");
634 params_->set(
"Output Style", outputStyle_);
635 outputTest_ = Teuchos::null;
639 if (params->isParameter(
"Output Stream")) {
640 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,
"Output Stream");
643 params_->set(
"Output Stream", outputStream_);
644 if (printer_ != Teuchos::null)
645 printer_->setOStream( outputStream_ );
650 if (params->isParameter(
"Output Frequency")) {
651 outputFreq_ = params->get(
"Output Frequency", outputFreq_default_);
655 params_->set(
"Output Frequency", outputFreq_);
656 if (outputTest_ != Teuchos::null)
657 outputTest_->setOutputFrequency( outputFreq_ );
661 if (printer_ == Teuchos::null) {
670 if (params->isParameter(
"Convergence Tolerance")) {
671 convtol_ = params->get(
"Convergence Tolerance",convtol_default_);
674 params_->set(
"Convergence Tolerance", convtol_);
675 if (convTest_ != Teuchos::null)
682 if (maxIterTest_ == Teuchos::null)
685 if (convTest_ == Teuchos::null)
686 convTest_ = Teuchos::rcp(
new StatusTestResNorm_t( convtol_, 1 ) );
688 sTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
696 std::string solverDesc =
" PCPG ";
697 outputTest_->setSolverDesc( solverDesc );
701 if (ortho_ == Teuchos::null) {
702 if (orthoType_==
"DGKS") {
703 if (orthoKappa_ <= 0) {
711 else if (orthoType_==
"ICGS") {
714 else if (orthoType_==
"IMGS") {
718 TEUCHOS_TEST_FOR_EXCEPTION(orthoType_!=
"ICGS"&&orthoType_!=
"DGKS"&&orthoType_!=
"IMGS",std::logic_error,
719 "Belos::PCPGSolMgr(): Invalid orthogonalization type.");
724 if (timerSolve_ == Teuchos::null) {
725 std::string solveLabel = label_ +
": PCPGSolMgr total solve time";
726 #ifdef BELOS_TEUCHOS_TIME_MONITOR 727 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
736 template<
class ScalarType,
class MV,
class OP>
737 Teuchos::RCP<const Teuchos::ParameterList>
740 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
741 if (is_null(validPL)) {
742 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
744 pl->set(
"Convergence Tolerance", convtol_default_,
745 "The relative residual tolerance that needs to be achieved by the\n" 746 "iterative solver in order for the linear system to be declared converged.");
747 pl->set(
"Maximum Iterations", maxIters_default_,
748 "The maximum number of iterations allowed for each\n" 749 "set of RHS solved.");
750 pl->set(
"Num Deflated Blocks", deflatedBlocks_default_,
751 "The maximum number of vectors in the seed subspace." );
752 pl->set(
"Num Saved Blocks", savedBlocks_default_,
753 "The maximum number of vectors saved from old Krylov subspaces." );
754 pl->set(
"Verbosity", verbosity_default_,
755 "What type(s) of solver information should be outputted\n" 756 "to the output stream.");
757 pl->set(
"Output Style", outputStyle_default_,
758 "What style is used for the solver information outputted\n" 759 "to the output stream.");
760 pl->set(
"Output Frequency", outputFreq_default_,
761 "How often convergence information should be outputted\n" 762 "to the output stream.");
763 pl->set(
"Output Stream", outputStream_default_,
764 "A reference-counted pointer to the output stream where all\n" 765 "solver output is sent.");
766 pl->set(
"Timer Label", label_default_,
767 "The string to use as a prefix for the timer labels.");
769 pl->set(
"Orthogonalization", orthoType_default_,
770 "The type of orthogonalization to use: DGKS, ICGS, IMGS");
771 pl->set(
"Orthogonalization Constant",orthoKappa_default_,
772 "The constant used by DGKS orthogonalization to determine\n" 773 "whether another step of classical Gram-Schmidt is necessary.");
781 template<
class ScalarType,
class MV,
class OP>
785 if (!isSet_) { setParameters( params_ ); }
787 Teuchos::BLAS<int,ScalarType> blas;
788 Teuchos::LAPACK<int,ScalarType> lapack;
789 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
790 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
793 "Belos::PCPGSolMgr::solve(): Linear problem is not a valid object.");
796 "Belos::PCPGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
799 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
800 std::vector<int> currIdx(1);
806 problem_->setLSIndex( currIdx );
809 bool isConverged =
true;
813 Teuchos::ParameterList plist;
814 plist.set(
"Saved Blocks", savedBlocks_);
815 plist.set(
"Block Size", 1);
816 plist.set(
"Keep Diagonal",
true);
817 plist.set(
"Initialize Diagonal",
true);
822 Teuchos::RCP<PCPGIter<ScalarType,MV,OP> > pcpg_iter;
828 #ifdef BELOS_TEUCHOS_TIME_MONITOR 829 Teuchos::TimeMonitor slvtimer(*timerSolve_);
831 while ( numRHS2Solve > 0 ) {
834 outputTest_->reset();
837 if (R_ == Teuchos::null)
838 R_ = MVT::Clone( *(problem_->getRHS()), 1 );
840 problem_->computeCurrResVec( &*R_ );
846 if( U_ != Teuchos::null ){
851 Teuchos::RCP<MV> cur_soln_vec = problem_->getCurrLHSVec();
852 std::vector<MagnitudeType> rnorm0(1);
853 MVT::MvNorm( *R_, rnorm0 );
856 std::cout <<
"Solver Manager: dimU_ = " << dimU_ << std::endl;
857 Teuchos::SerialDenseMatrix<int,ScalarType> Z( dimU_, 1 );
859 Teuchos::RCP<const MV> Uactive, Cactive;
860 std::vector<int> active_columns( dimU_ );
861 for (
int i=0; i < dimU_; ++i) active_columns[i] = i;
862 Uactive = MVT::CloneView(*U_, active_columns);
863 Cactive = MVT::CloneView(*C_, active_columns);
866 std::cout <<
" Solver Manager : check duality of seed basis " << std::endl;
867 Teuchos::SerialDenseMatrix<int,ScalarType> H( dimU_, dimU_ );
868 MVT::MvTransMv( one, *Uactive, *Cactive, H );
869 H.print( std::cout );
872 MVT::MvTransMv( one, *Uactive, *R_, Z );
873 Teuchos::RCP<MV> tempU = MVT::Clone( *R_, 1 );
874 MVT::MvTimesMatAddMv( one, *Uactive, Z, zero, *tempU );
875 MVT::MvAddMv( one, *tempU, one, *cur_soln_vec, *cur_soln_vec );
876 MVT::MvTimesMatAddMv( one, *Cactive, Z, zero, *tempU );
877 MVT::MvAddMv( -one, *tempU, one, *R_, *R_ );
878 std::vector<MagnitudeType> rnorm(1);
879 MVT::MvNorm( *R_, rnorm );
880 if( rnorm[0] < rnorm0[0] * .001 ){
881 MVT::MvTransMv( one, *Uactive, *R_, Z );
882 MVT::MvTimesMatAddMv( one, *Uactive, Z, zero, *tempU );
883 MVT::MvAddMv( one, *tempU, one, *cur_soln_vec, *cur_soln_vec );
884 MVT::MvTimesMatAddMv( one, *Cactive, Z, zero, *tempU );
885 MVT::MvAddMv( -one, *tempU, one, *R_, *R_ );
887 Uactive = Teuchos::null;
888 Cactive = Teuchos::null;
889 tempU = Teuchos::null;
900 if( U_ != Teuchos::null ) pcpgState.
U = U_;
901 if( C_ != Teuchos::null ) pcpgState.
C = C_;
902 if( dimU_ > 0 ) pcpgState.
curDim = dimU_;
903 pcpg_iter->initialize(pcpgState);
909 if( !dimU_ ) printer_->stream(
Debug) <<
" No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
910 pcpg_iter->resetNumIters();
912 if( dimU_ > savedBlocks_ )
913 std::cout <<
"Error: dimU_ = " << dimU_ <<
" > savedBlocks_ = " << savedBlocks_ << std::endl;
919 if( debug ) printf(
"********** Calling iterate...\n");
920 pcpg_iter->iterate();
927 if ( convTest_->getStatus() ==
Passed ) {
936 else if ( maxIterTest_->getStatus() ==
Passed ) {
950 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
951 "Belos::PCPGSolMgr::solve(): Invalid return from PCPGIter::iterate().");
957 sTest_->checkStatus( &*pcpg_iter );
958 if (convTest_->getStatus() !=
Passed)
962 catch (
const std::exception &e) {
963 printer_->stream(
Errors) <<
"Error! Caught exception in PCPGIter::iterate() at iteration " 964 << pcpg_iter->getNumIters() << std::endl
965 << e.what() << std::endl;
971 Teuchos::RCP<MV> update = pcpg_iter->getCurrentUpdate();
972 problem_->updateSolution( update,
true );
975 problem_->setCurrLS();
983 std::cout <<
"SolverManager: dimU_ " << dimU_ <<
" prevUdim= " << q << std::endl;
985 if( q > deflatedBlocks_ )
986 std::cout <<
"SolverManager: Error deflatedBlocks = " << deflatedBlocks_ << std::endl;
997 rank = ARRQR(dimU_,q, *oldState.
D );
999 std::cout <<
" rank decreased in ARRQR, something to do? " << std::endl;
1005 if( dimU_ > deflatedBlocks_ ){
1007 if( !deflatedBlocks_ ){
1010 dimU_ = deflatedBlocks_;
1014 bool Harmonic =
false;
1016 Teuchos::RCP<MV> Uorth;
1018 std::vector<int> active_cols( dimU_ );
1019 for (
int i=0; i < dimU_; ++i) active_cols[i] = i;
1022 Uorth = MVT::CloneCopy(*C_, active_cols);
1025 Uorth = MVT::CloneCopy(*U_, active_cols);
1029 Teuchos::SerialDenseMatrix<int,ScalarType> R(dimU_,dimU_);
1030 rank = ortho_->normalize(*Uorth, Teuchos::rcp(&R,
false));
1031 Uorth = Teuchos::null;
1037 "Belos::PCPGSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
1041 Teuchos::SerialDenseMatrix<int,ScalarType> VT;
1042 Teuchos::SerialDenseMatrix<int,ScalarType> Ur;
1043 int lwork = 5*dimU_;
1046 if( problem_->isHermitian() ) lrwork = dimU_;
1047 std::vector<ScalarType> work(lwork);
1048 std::vector<ScalarType> Svec(dimU_);
1049 std::vector<ScalarType> rwork(lrwork);
1050 lapack.GESVD(
'N',
'O',
1051 R.numRows(),R.numCols(),R.values(), R.numRows(),
1059 "Belos::PCPGSolMgr::solve(): LAPACK _GESVD failed to compute singular values.");
1061 if( work[0] != 67. * dimU_ )
1062 std::cout <<
" SVD " << dimU_ <<
" lwork " << work[0] << std::endl;
1063 for(
int i=0; i< dimU_; i++)
1064 std::cout << i <<
" " << Svec[i] << std::endl;
1066 Teuchos::SerialDenseMatrix<int,ScalarType> wholeV( R,
Teuchos::TRANS);
1068 int startRow = 0, startCol = 0;
1070 startCol = dimU_ - deflatedBlocks_;
1072 Teuchos::SerialDenseMatrix<int,ScalarType> V(Teuchos::Copy,
1078 std::vector<int> active_columns( dimU_ );
1079 std::vector<int> def_cols( deflatedBlocks_ );
1080 for (
int i=0; i < dimU_; ++i) active_columns[i] = i;
1081 for (
int i=0; i < deflatedBlocks_; ++i) def_cols[i] = i;
1083 Teuchos::RCP<MV> Uactive = MVT::CloneViewNonConst(*U_, def_cols);
1084 Teuchos::RCP<MV> Ucopy = MVT::CloneCopy( *U_, active_columns );
1085 MVT::MvTimesMatAddMv( one, *Ucopy, V, zero, *Uactive );
1086 Ucopy = Teuchos::null;
1087 Uactive = Teuchos::null;
1088 Teuchos::RCP<MV> Cactive = MVT::CloneViewNonConst(*C_, def_cols);
1089 Teuchos::RCP<MV> Ccopy = MVT::CloneCopy( *C_, active_columns );
1090 MVT::MvTimesMatAddMv( one, *Ccopy, V, zero, *Cactive );
1091 Ccopy = Teuchos::null;
1092 Cactive = Teuchos::null;
1093 dimU_ = deflatedBlocks_;
1095 printer_->stream(
Debug) <<
" Generated recycled subspace using RHS index " << currIdx[0] <<
" of dimension " << dimU_ << std::endl << std::endl;
1098 problem_->setCurrLS();
1102 if ( numRHS2Solve > 0 ) {
1106 problem_->setLSIndex( currIdx );
1109 currIdx.resize( numRHS2Solve );
1118 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1123 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails) );
1128 using Teuchos::rcp_dynamic_cast;
1131 const std::vector<MagnitudeType>* pTestValues =
1132 rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
1134 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1135 "Belos::PCPGSolMgr::solve(): The convergence test's getTestValue() " 1136 "method returned NULL. Please report this bug to the Belos developers.");
1138 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1139 "Belos::PCPGSolMgr::solve(): The convergence test's getTestValue() " 1140 "method returned a vector of length zero. Please report this bug to the " 1141 "Belos developers.");
1146 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1150 numIters_ = maxIterTest_->getNumIters();
1161 template<
class ScalarType,
class MV,
class OP>
1165 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1166 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1169 Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
1170 Teuchos::SerialDenseMatrix<int,ScalarType> gamma( 1, 1 );
1171 Teuchos::SerialDenseMatrix<int,ScalarType> anorm( 1, 1 );
1172 std::vector<int> curind(1);
1173 std::vector<int> ipiv(p - q);
1174 std::vector<ScalarType> Pivots(p);
1176 ScalarType rteps = 1.5e-8;
1179 for( i = q ; i < p ; i++ ){
1182 RCP<MV> P = MVT::CloneViewNonConst(*U_,curind);
1183 RCP<MV> AP = MVT::CloneViewNonConst(*C_,curind);
1184 anorm(0,0) = one / Teuchos::ScalarTraits<ScalarType>::squareroot( D(i-q,i-q) ) ;
1185 MVT::MvAddMv( anorm(0,0), *P, zero, *AP, *P );
1186 MVT::MvAddMv( zero, *P, anorm(0,0), *AP, *AP );
1190 for( i = q ; i < p ; i++ ){
1191 if( q < i && i < p-1 ){
1194 for( j = i+1 ; j < p ; j++ ){
1195 const int k = ipiv[j-q];
1196 if( Pivots[k] > Pivots[l] ){
1203 ipiv[imax-q] = ipiv[i-q];
1209 if( Pivots[k] > 1.5625e-2 ){
1210 anorm(0,0) = Pivots[k];
1214 RCP<const MV> P = MVT::CloneView(*U_,curind);
1215 RCP<const MV> AP = MVT::CloneView(*C_,curind);
1216 MVT::MvTransMv( one, *P, *AP, anorm );
1217 anorm(0,0) = Teuchos::ScalarTraits<ScalarType>::squareroot( anorm(0,0) ) ;
1219 if( rteps <= anorm(0,0) && anorm(0,0) < 9.765625e-4){
1227 std::cout <<
"ARRQR: Bad case not implemented" << std::endl;
1229 if( anorm(0,0) < rteps ){
1230 std::cout <<
"ARRQR : deficient case not implemented " << std::endl;
1238 RCP<MV> P = MVT::CloneViewNonConst(*U_,curind);
1239 RCP<MV> AP = MVT::CloneViewNonConst(*C_,curind);
1240 MVT::MvAddMv( anorm(0,0), *P, zero, *AP, *P );
1241 MVT::MvAddMv( zero, *P, anorm(0,0), *AP, *AP );
1245 P = MVT::CloneViewNonConst(*U_,curind);
1246 AP = MVT::CloneViewNonConst(*C_,curind);
1247 for( j = i+1 ; j < p ; j++ ){
1250 RCP<MV> Q = MVT::CloneViewNonConst(*U_,curind);
1251 MVT::MvTransMv( one, *Q, *AP, alpha);
1252 MVT::MvAddMv( -alpha(0,0), *P, one, *Q, *Q );
1254 RCP<MV> AQ = MVT::CloneViewNonConst(*C_,curind);
1255 MVT::MvAddMv( -alpha(0,0), *AP, one, *AQ, *AQ );
1257 gamma(0,0) = ( Pivots[l] - alpha(0,0))*( Pivots[l] + alpha(0,0));
1258 if( gamma(0,0) > 0){
1259 Pivots[l] = Teuchos::ScalarTraits<ScalarType>::squareroot( gamma(0,0) );
1270 template<
class ScalarType,
class MV,
class OP>
1273 std::ostringstream oss;
1274 oss <<
"Belos::PCPGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<
">";
1276 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.
static const std::string orthoType_default_
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
static const MagnitudeType orthoKappa_default_
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.
static const bool scalarTypeIsSupported
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
Belos concrete class to iterate Preconditioned Conjugate Projected Gradients.
int numIters_
Number of iterations taken by the last solve() invocation.
PCPGSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate orthonormal...
Teuchos::RCP< MV > R
The current residual.
static const int verbosity_default_
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)
int maxIters_
Maximum iteration count (read from parameter list).
MagnitudeType orthoKappa_
Orthogonalization parameter (read from parameter list).
Structure to contain pointers to PCPGIter state variables.
Teuchos::RCP< Teuchos::Time > timerSolve_
PCPGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i.e.
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
Belos::StatusTest class for specifying a maximum number of iterations.
Teuchos::RCP< std::ostream > outputStream_
static const MagnitudeType convtol_default_
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.
OperatorTraits< ScalarType, MV, OP > OPT
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. ...
MagnitudeType convtol_
Convergence tolerance (read from parameter list).
void reset(const ResetType type)
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
static const int savedBlocks_default_
MagnitudeType achievedTol_
Tolerance achieved by the last solve() invocation.
A linear system to solve, and its associated information.
MultiVecTraits< ScalarType, MV > MVT
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem)
Set the linear problem that needs to be solved.
static const int outputFreq_default_
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)
static const int maxIters_default_
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
static const Teuchos::RCP< std::ostream > outputStream_default_
ReturnType
Whether the Belos solve converged for all linear systems.
static const int deflatedBlocks_default_
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< OutputManager< ScalarType > > printer_
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::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
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...
Teuchos::ScalarTraits< ScalarType > SCT
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
PCPG iterative linear solver.
static const std::string label_default_
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.
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > convTest_
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
virtual ~PCPGSolMgr()
Destructor.
Class which defines basic traits for the operator type.
Parent class to all Belos exceptions.
Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > ortho_
Details::SolverManagerRequiresRealLapack< ScalarType, MV, OP, scalarTypeIsSupported > base_type
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...
static const int outputStyle_default_
Teuchos::RCP< Teuchos::ParameterList > params_
PCPGSolMgrLinearProblemFailure(const std::string &what_arg)
Teuchos::ScalarTraits< MagnitudeType > MT