42 #ifndef BELOS_RCG_SOLMGR_HPP
43 #define BELOS_RCG_SOLMGR_HPP
63 #ifdef BELOS_TEUCHOS_TIME_MONITOR
136 template<
class ScalarType,
class MV,
class OP,
137 const bool supportsScalarType =
142 Belos::Details::LapackSupportsScalar<ScalarType>::value &&
143 ! Teuchos::ScalarTraits<ScalarType>::isComplex>
170 template<
class ScalarType,
class MV,
class OP>
244 return Teuchos::tuple(timerSolve_);
318 std::string description()
const override;
334 void sort(std::vector<ScalarType>& dlist,
int n, std::vector<int>& iperm);
337 void initializeStateStorage();
356 static constexpr
int maxIters_default_ = 1000;
357 static constexpr
int blockSize_default_ = 1;
358 static constexpr
int numBlocks_default_ = 25;
359 static constexpr
int recycleBlocks_default_ = 3;
360 static constexpr
bool showMaxResNormOnly_default_ =
false;
363 static constexpr
int outputFreq_default_ = -1;
364 static constexpr
const char * label_default_ =
"Belos";
462 template<
class ScalarType,
class MV,
class OP>
471 template<
class ScalarType,
class MV,
class OP>
489 template<
class ScalarType,
class MV,
class OP>
492 outputStream_ = Teuchos::rcpFromRef(std::cout);
494 maxIters_ = maxIters_default_;
495 numBlocks_ = numBlocks_default_;
496 recycleBlocks_ = recycleBlocks_default_;
497 verbosity_ = verbosity_default_;
498 outputStyle_= outputStyle_default_;
499 outputFreq_= outputFreq_default_;
500 showMaxResNormOnly_ = showMaxResNormOnly_default_;
501 label_ = label_default_;
540 template<
class ScalarType,
class MV,
class OP>
553 maxIters_ = params->
get(
"Maximum Iterations",maxIters_default_);
556 params_->set(
"Maximum Iterations", maxIters_);
558 maxIterTest_->setMaxIters( maxIters_ );
563 numBlocks_ = params->
get(
"Num Blocks",numBlocks_default_);
565 "Belos::RCGSolMgr: \"Num Blocks\" must be strictly positive.");
568 params_->set(
"Num Blocks", numBlocks_);
573 recycleBlocks_ = params->
get(
"Num Recycled Blocks",recycleBlocks_default_);
575 "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be strictly positive.");
578 "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be less than \"Num Blocks\".");
581 params_->set(
"Num Recycled Blocks", recycleBlocks_);
586 std::string tempLabel = params->
get(
"Timer Label", label_default_);
589 if (tempLabel != label_) {
591 params_->set(
"Timer Label", label_);
592 std::string solveLabel = label_ +
": RCGSolMgr total solve time";
593 #ifdef BELOS_TEUCHOS_TIME_MONITOR
601 if (Teuchos::isParameterType<int>(*params,
"Verbosity")) {
602 verbosity_ = params->
get(
"Verbosity", verbosity_default_);
604 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,
"Verbosity");
608 params_->set(
"Verbosity", verbosity_);
610 printer_->setVerbosity(verbosity_);
615 if (Teuchos::isParameterType<int>(*params,
"Output Style")) {
616 outputStyle_ = params->
get(
"Output Style", outputStyle_default_);
618 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,
"Output Style");
622 params_->set(
"Output Style", outputStyle_);
628 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,
"Output Stream");
631 params_->set(
"Output Stream", outputStream_);
633 printer_->setOStream( outputStream_ );
639 outputFreq_ = params->
get(
"Output Frequency", outputFreq_default_);
643 params_->set(
"Output Frequency", outputFreq_);
645 outputTest_->setOutputFrequency( outputFreq_ );
658 if (params->
isParameter(
"Convergence Tolerance")) {
660 convtol_ = params->
get (
"Convergence Tolerance",
668 params_->set(
"Convergence Tolerance", convtol_);
670 convTest_->setTolerance( convtol_ );
673 if (params->
isParameter(
"Show Maximum Residual Norm Only")) {
674 showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,
"Show Maximum Residual Norm Only");
677 params_->set(
"Show Maximum Residual Norm Only", showMaxResNormOnly_);
679 convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
690 convTest_ =
Teuchos::rcp(
new StatusTestResNorm_t( convtol_, 1 ) );
693 sTest_ =
Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
703 std::string solverDesc =
" Recycling CG ";
704 outputTest_->setSolverDesc( solverDesc );
709 std::string solveLabel = label_ +
": RCGSolMgr total solve time";
710 #ifdef BELOS_TEUCHOS_TIME_MONITOR
720 template<
class ScalarType,
class MV,
class OP>
730 "The relative residual tolerance that needs to be achieved by the\n"
731 "iterative solver in order for the linear system to be declared converged.");
732 pl->
set(
"Maximum Iterations", static_cast<int>(maxIters_default_),
733 "The maximum number of block iterations allowed for each\n"
734 "set of RHS solved.");
735 pl->
set(
"Block Size", static_cast<int>(blockSize_default_),
736 "Block Size Parameter -- currently must be 1 for RCG");
737 pl->
set(
"Num Blocks", static_cast<int>(numBlocks_default_),
738 "The length of a cycle (and this max number of search vectors kept)\n");
739 pl->
set(
"Num Recycled Blocks", static_cast<int>(recycleBlocks_default_),
740 "The number of vectors in the recycle subspace.");
741 pl->
set(
"Verbosity", static_cast<int>(verbosity_default_),
742 "What type(s) of solver information should be outputted\n"
743 "to the output stream.");
744 pl->
set(
"Output Style", static_cast<int>(outputStyle_default_),
745 "What style is used for the solver information outputted\n"
746 "to the output stream.");
747 pl->
set(
"Output Frequency", static_cast<int>(outputFreq_default_),
748 "How often convergence information should be outputted\n"
749 "to the output stream.");
750 pl->
set(
"Output Stream", Teuchos::rcpFromRef(std::cout),
751 "A reference-counted pointer to the output stream where all\n"
752 "solver output is sent.");
753 pl->
set(
"Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
754 "When convergence information is printed, only show the maximum\n"
755 "relative residual norm when the block size is greater than one.");
756 pl->
set(
"Timer Label", static_cast<const char *>(label_default_),
757 "The string to use as a prefix for the timer labels.");
764 template<
class ScalarType,
class MV,
class OP>
777 "Belos::RCGSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
781 P_ = MVT::Clone( *rhsMV, numBlocks_+2 );
785 if (MVT::GetNumberVecs(*P_) < numBlocks_+2) {
787 P_ = MVT::Clone( *tmp, numBlocks_+2 );
793 Ap_ = MVT::Clone( *rhsMV, 1 );
797 r_ = MVT::Clone( *rhsMV, 1 );
801 z_ = MVT::Clone( *rhsMV, 1 );
805 U_ = MVT::Clone( *rhsMV, recycleBlocks_ );
809 if (MVT::GetNumberVecs(*U_) < recycleBlocks_) {
811 U_ = MVT::Clone( *tmp, recycleBlocks_ );
817 AU_ = MVT::Clone( *rhsMV, recycleBlocks_ );
821 if (MVT::GetNumberVecs(*AU_) < recycleBlocks_) {
823 AU_ = MVT::Clone( *tmp, recycleBlocks_ );
829 U1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
833 if (MVT::GetNumberVecs(*U1_) < recycleBlocks_) {
835 U1_ = MVT::Clone( *tmp, recycleBlocks_ );
843 if ( (Alpha_->numRows() != numBlocks_) || (Alpha_->numCols() != 1) )
844 Alpha_->reshape( numBlocks_, 1 );
851 if ( (Beta_->numRows() != (numBlocks_+1)) || (Beta_->numCols() != 1) )
852 Beta_->reshape( numBlocks_ + 1, 1 );
859 if ( (D_->numRows() != numBlocks_) || (D_->numCols() != 1) )
860 D_->reshape( numBlocks_, 1 );
867 if ( (Delta_->numRows() != recycleBlocks_) || (Delta_->numCols() != (numBlocks_ + 1)) )
868 Delta_->reshape( recycleBlocks_, numBlocks_ + 1 );
875 if ( (UTAU_->numRows() != recycleBlocks_) || (UTAU_->numCols() != recycleBlocks_) )
876 UTAU_->reshape( recycleBlocks_, recycleBlocks_ );
883 if ( (LUUTAU_->numRows() != recycleBlocks_) || (LUUTAU_->numCols() != recycleBlocks_) )
884 LUUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
889 ipiv_ =
Teuchos::rcp(
new std::vector<int>(recycleBlocks_) );
891 if ( (
int)ipiv_->size() != recycleBlocks_ )
892 ipiv_->resize(recycleBlocks_);
899 if ( (AUTAU_->numRows() != recycleBlocks_) || (AUTAU_->numCols() != recycleBlocks_) )
900 AUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
907 if ( (rTz_old_->numRows() != 1) || (rTz_old_->numCols() != 1) )
908 rTz_old_->reshape( 1, 1 );
915 if ( (F_->numRows() != (numBlocks_+recycleBlocks_)) || (F_->numCols() != numBlocks_+recycleBlocks_) )
916 F_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
923 if ( (G_->numRows() != (numBlocks_+recycleBlocks_)) || (G_->numCols() != numBlocks_+recycleBlocks_) )
924 G_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
931 if ( (Y_->numRows() != (numBlocks_+recycleBlocks_)) || (Y_->numCols() != numBlocks_+recycleBlocks_) )
932 Y_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
939 if ( (L2_->numRows() != (numBlocks_+1)) || (L2_->numCols() != numBlocks_) )
940 L2_->reshape( numBlocks_+1, numBlocks_ );
947 if ( (DeltaL2_->numRows() != (recycleBlocks_)) || (DeltaL2_->numCols() != (numBlocks_) ) )
948 DeltaL2_->reshape( recycleBlocks_, numBlocks_ );
955 if ( (AU1TUDeltaL2_->numRows() != (recycleBlocks_)) || (AU1TUDeltaL2_->numCols() != (numBlocks_) ) )
956 AU1TUDeltaL2_->reshape( recycleBlocks_, numBlocks_ );
963 if ( (AU1TAU1_->numRows() != (recycleBlocks_)) || (AU1TAU1_->numCols() != (recycleBlocks_) ) )
964 AU1TAU1_->reshape( recycleBlocks_, recycleBlocks_ );
971 if ( (GY_->numRows() != (numBlocks_ + recycleBlocks_)) || (GY_->numCols() != (recycleBlocks_) ) )
972 GY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
979 if ( (AU1TU1_->numRows() != (recycleBlocks_)) || (AU1TU1_->numCols() != (recycleBlocks_) ) )
980 AU1TU1_->reshape( recycleBlocks_, recycleBlocks_ );
987 if ( (FY_->numRows() != (numBlocks_ + recycleBlocks_)) || (FY_->numCols() != (recycleBlocks_) ) )
988 FY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
995 if ( (AU1TAP_->numRows() != (recycleBlocks_)) || (AU1TAP_->numCols() != (numBlocks_) ) )
996 AU1TAP_->reshape( recycleBlocks_, numBlocks_ );
1003 if ( (APTAP_->numRows() != (numBlocks_)) || (APTAP_->numCols() != (numBlocks_) ) )
1004 APTAP_->reshape( numBlocks_, numBlocks_ );
1009 U1Y1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
1013 if (MVT::GetNumberVecs(*U1Y1_) < recycleBlocks_) {
1015 U1Y1_ = MVT::Clone( *tmp, recycleBlocks_ );
1021 PY2_ = MVT::Clone( *rhsMV, recycleBlocks_ );
1025 if (MVT::GetNumberVecs(*PY2_) < recycleBlocks_) {
1027 PY2_ = MVT::Clone( *tmp, recycleBlocks_ );
1035 if ( (AUTAP_->numRows() != (recycleBlocks_)) || (AUTAP_->numCols() != (numBlocks_) ) )
1036 AUTAP_->reshape( recycleBlocks_, numBlocks_ );
1043 if ( (AU1TU_->numRows() != (recycleBlocks_)) || (AU1TU_->numCols() != (recycleBlocks_) ) )
1044 AU1TU_->reshape( recycleBlocks_, recycleBlocks_ );
1051 template<
class ScalarType,
class MV,
class OP>
1055 std::vector<int> index(recycleBlocks_);
1066 setParameters(Teuchos::parameterList(*getValidParameters()));
1070 "Belos::RCGSolMgr::solve(): Linear problem is not a valid object.");
1072 "Belos::RCGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
1075 "Belos::RCGSolMgr::solve(): RCG does not support split preconditioning, only set left or right preconditioner.");
1080 precObj = Teuchos::rcp_const_cast<OP>(problem_->getLeftPrec());
1083 precObj = Teuchos::rcp_const_cast<OP>(problem_->getRightPrec());
1087 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1088 std::vector<int> currIdx(1);
1092 problem_->setLSIndex( currIdx );
1095 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1096 if (numBlocks_ > dim) {
1097 numBlocks_ = Teuchos::asSafe<int>(dim);
1098 params_->set(
"Num Blocks", numBlocks_);
1100 "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
1101 " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
1105 initializeStateStorage();
1109 plist.
set(
"Num Blocks",numBlocks_);
1110 plist.
set(
"Recycled Blocks",recycleBlocks_);
1113 outputTest_->reset();
1116 bool isConverged =
true;
1120 index.resize(recycleBlocks_);
1121 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1123 index.resize(recycleBlocks_);
1124 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1127 problem_->applyOp( *Utmp, *AUtmp );
1130 MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
1134 index.resize(recycleBlocks_);
1135 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1136 index.resize(recycleBlocks_);
1137 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1139 OPT::Apply( *precObj, *AUtmp, *PCAU );
1140 MVT::MvTransMv( one, *AUtmp, *PCAU, AUTAUtmp );
1142 MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
1154 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1158 while ( numRHS2Solve > 0 ) {
1161 if (printer_->isVerbosity(
Debug ) ) {
1162 if (existU_) printer_->print(
Debug,
"Using recycle space generated from previous call to solve()." );
1163 else printer_->print(
Debug,
"No recycle space exists." );
1167 rcg_iter->resetNumIters();
1170 rcg_iter->setSize( recycleBlocks_, numBlocks_ );
1173 outputTest_->resetNumCalls();
1182 problem_->computeCurrResVec( &*r_ );
1188 index.resize(recycleBlocks_);
1189 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1191 MVT::MvTransMv( one, *Utmp, *r_, Utr );
1194 LUUTAUtmp.
assign(UTAUtmp);
1198 "Belos::RCGSolMgr::solve(): LAPACK GESV failed to compute a solution.");
1201 MVT::MvTimesMatAddMv( one, *Utmp, Utr, one, *problem_->getCurrLHSVec() );
1204 index.resize(recycleBlocks_);
1205 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1207 MVT::MvTimesMatAddMv( -one, *AUtmp, Utr, one, *r_ );
1211 OPT::Apply( *precObj, *r_, *z_ );
1217 MVT::MvTransMv( one, *r_, *z_, *rTz_old_ );
1222 index.resize(recycleBlocks_);
1223 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1225 MVT::MvTransMv( one, *AUtmp, *z_, mu );
1231 "Belos::RCGSolMgr::solve(): LAPACK GETRS failed to compute a solution.");
1236 MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
1237 MVT::MvTimesMatAddMv( -one, *U_, mu, one, *Ptmp );
1243 MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
1250 index.resize( numBlocks_+1 );
1251 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1252 newstate.
P = MVT::CloneViewNonConst( *P_, index );
1253 index.resize( recycleBlocks_ );
1254 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1255 newstate.
U = MVT::CloneViewNonConst( *U_, index );
1256 index.resize( recycleBlocks_ );
1257 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1258 newstate.
AU = MVT::CloneViewNonConst( *AU_, index );
1269 newstate.
existU = existU_;
1270 newstate.
ipiv = ipiv_;
1273 rcg_iter->initialize(newstate);
1279 rcg_iter->iterate();
1286 if ( convTest_->getStatus() ==
Passed ) {
1295 else if ( maxIterTest_->getStatus() ==
Passed ) {
1297 isConverged =
false;
1305 else if ( rcg_iter->getCurSubspaceDim() == rcg_iter->getMaxSubspaceDim() ) {
1310 if (recycleBlocks_ > 0) {
1321 for (
int ii=0;ii<numBlocks_;ii++) {
1322 Gtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
1324 Gtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1325 Gtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1327 Ftmp(ii,ii) = Dtmp(ii,0);
1332 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1335 index.resize( numBlocks_ );
1336 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
1338 index.resize( recycleBlocks_ );
1339 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1341 MVT::MvTimesMatAddMv( one, *Ptmp, Ytmp, zero, *U1tmp );
1361 ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
1362 for (
int ii=0; ii<recycleBlocks_; ++ii) {
1363 AU1TAPtmp(ii,0) = Ytmp(numBlocks_-1,ii) * alphatmp;
1371 lastBeta = numBlocks_-1;
1380 AU1TAPtmp.
scale(Dtmp(0,0));
1386 for (
int ii=0; ii<numBlocks_; ii++) {
1387 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
1389 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1390 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1405 for(
int ii=0;ii<numBlocks_;ii++) {
1406 F22(ii,ii) = Dtmp(ii,0);
1419 for (
int ii=0;ii<recycleBlocks_;++ii)
1420 for (
int jj=0;jj<numBlocks_;++jj)
1421 G21(jj,ii) = G12(ii,jj);
1426 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1429 index.resize( numBlocks_ );
1430 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
1432 index.resize( recycleBlocks_ );
1433 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1436 index.resize( recycleBlocks_ );
1437 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1439 index.resize( recycleBlocks_ );
1440 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1443 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1444 MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
1445 MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
1457 ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
1458 for (
int ii=0; ii<recycleBlocks_; ++ii) {
1459 AU1TAPtmp(ii,0) = Ytmp(numBlocks_+recycleBlocks_-1,ii) * alphatmp;
1469 lastp = numBlocks_+1;
1470 lastBeta = numBlocks_;
1481 for (
int ii=0; ii<numBlocks_; ii++) {
1482 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
1484 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1485 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1491 for(
int ii=0;ii<numBlocks_;ii++) {
1492 L2tmp(ii,ii) = 1./Alphatmp(ii,0);
1493 L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
1514 for(
int ii=0;ii<numBlocks_;ii++) {
1515 F22(ii,ii) = Dtmp(ii,0);
1528 for (
int ii=0;ii<recycleBlocks_;++ii)
1529 for (
int jj=0;jj<numBlocks_;++jj)
1530 G21(jj,ii) = G12(ii,jj);
1535 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1538 index.resize( recycleBlocks_ );
1539 for (
int ii=0; ii<(recycleBlocks_); ++ii) { index[ii] = ii; }
1541 index.resize( numBlocks_ );
1542 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
1544 index.resize( recycleBlocks_ );
1545 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1548 index.resize( recycleBlocks_ );
1549 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1552 index.resize( recycleBlocks_ );
1553 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1555 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1556 MVT::MvTimesMatAddMv( one, *Utmp, Y1, zero, *UY1tmp );
1557 MVT::MvAddMv(one,*UY1tmp, one, *PY2tmp, *U1tmp);
1575 AU1TUtmp.
assign(UTAUtmp);
1578 dold = Dtmp(numBlocks_-1,0);
1585 lastBeta = numBlocks_-1;
1592 for (
int ii=0; ii<numBlocks_; ii++) {
1593 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
1595 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1596 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1601 for(
int ii=0;ii<numBlocks_;ii++) {
1602 L2tmp(ii,ii) = 1./Alphatmp(ii,0);
1603 L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
1619 ScalarType val = dold * (-Betatmp(0,0)/Alphatmp(0,0));
1620 for(
int ii=0;ii<recycleBlocks_;ii++) {
1621 AU1TAPtmp(ii,0) += Y2(numBlocks_-1,ii)*val;
1627 AU1TUtmp.
assign(Y1TAU1TU);
1637 for(
int ii=0;ii<numBlocks_;ii++) {
1638 F22(ii,ii) = Dtmp(ii,0);
1651 for (
int ii=0;ii<recycleBlocks_;++ii)
1652 for (
int jj=0;jj<numBlocks_;++jj)
1653 G21(jj,ii) = G12(ii,jj);
1658 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1661 index.resize( numBlocks_ );
1662 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
1664 index.resize( recycleBlocks_ );
1665 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1667 index.resize( recycleBlocks_ );
1668 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1670 index.resize( recycleBlocks_ );
1671 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1673 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1674 MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
1675 MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
1690 dold = Dtmp(numBlocks_-1,0);
1693 lastp = numBlocks_+1;
1694 lastBeta = numBlocks_;
1704 index[0] = lastp-1; index[1] = lastp;
1706 index[0] = 0; index[1] = 1;
1707 MVT::SetBlock(*Ptmp2,index,*P_);
1710 (*Beta_)(0,0) = (*Beta_)(lastBeta,0);
1721 index.resize( numBlocks_+1 );
1722 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii+1; }
1723 newstate.
P = MVT::CloneViewNonConst( *P_, index );
1734 rcg_iter->initialize(newstate);
1747 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
1748 "Belos::RCGSolMgr::solve(): Invalid return from RCGIter::iterate().");
1753 achievedTol_ = MT::one();
1755 MVT::MvInit( *X, SCT::zero() );
1756 printer_->stream(
Warnings) <<
"Belos::RCGSolMgr::solve(): Warning! NaN has been detected!"
1760 catch (
const std::exception &e) {
1761 printer_->stream(
Errors) <<
"Error! Caught std::exception in RCGIter::iterate() at iteration "
1762 << rcg_iter->getNumIters() << std::endl
1763 << e.what() << std::endl;
1769 problem_->setCurrLS();
1773 if ( numRHS2Solve > 0 ) {
1776 problem_->setLSIndex( currIdx );
1779 currIdx.resize( numRHS2Solve );
1785 index.resize(recycleBlocks_);
1786 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1787 MVT::SetBlock(*U1_,index,*U_);
1790 if (numRHS2Solve > 0) {
1807 index.resize(recycleBlocks_);
1808 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1810 index.resize(recycleBlocks_);
1811 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1814 problem_->applyOp( *Utmp, *AUtmp );
1817 MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
1821 index.resize(recycleBlocks_);
1822 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1823 index.resize(recycleBlocks_);
1824 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1826 OPT::Apply( *precObj, *AUtmp, *LeftPCAU );
1827 MVT::MvTransMv( one, *AUtmp, *LeftPCAU, AUTAUtmp );
1829 MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
1842 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1851 numIters_ = maxIterTest_->getNumIters();
1855 using Teuchos::rcp_dynamic_cast;
1858 const std::vector<MagnitudeType>* pTestValues =
1859 rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
1861 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1862 "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() "
1863 "method returned NULL. Please report this bug to the Belos developers.");
1865 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1866 "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() "
1867 "method returned a vector of length zero. Please report this bug to the "
1868 "Belos developers.");
1873 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1883 template<
class ScalarType,
class MV,
class OP>
1894 std::vector<MagnitudeType> w(n);
1897 std::vector<int> iperm(n);
1904 std::vector<ScalarType> work(1);
1912 lapack.
SYGV(itype, jobz, uplo, n, G2.
values(), G2.
stride(), F2.
values(), F2.
stride(), &w[0], &work[0], lwork, &info);
1914 "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to query optimal work size.");
1915 lwork = (int)work[0];
1917 lapack.
SYGV(itype, jobz, uplo, n, G2.
values(), G2.
stride(), F2.
values(), F2.
stride(), &w[0], &work[0], lwork, &info);
1919 "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to compute eigensolutions.");
1923 this->sort(w,n,iperm);
1926 for(
int i=0; i<recycleBlocks_; i++ ) {
1927 for(
int j=0; j<n; j++ ) {
1928 Y(j,i) = G2(j,iperm[i]);
1935 template<
class ScalarType,
class MV,
class OP>
1938 int l, r, j, i, flag;
1967 if (dlist[j] > dlist[j - 1]) j = j + 1;
1969 if (dlist[j - 1] > dK) {
1970 dlist[i - 1] = dlist[j - 1];
1971 iperm[i - 1] = iperm[j - 1];
1984 dlist[r] = dlist[0];
1985 iperm[r] = iperm[0];
2000 template<
class ScalarType,
class MV,
class OP>
2003 std::ostringstream oss;
ScalarType * values() const
Teuchos::RCP< MV > r
The current residual.
Collection of types and exceptions used within the Belos solvers.
Teuchos::RCP< std::ostream > outputStream_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > rTz_old
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
Teuchos::RCP< MV > P
The current Krylov basis.
Teuchos::ScalarTraits< MagnitudeType > MT
Teuchos::RCP< std::vector< int > > ipiv
Data from LU factorization of U^T A U.
Class which manages the output and verbosity of the Belos solvers.
Teuchos::RCP< Teuchos::Time > timerSolve_
bool is_null(const boost::shared_ptr< T > &p)
static const bool scalarTypeIsSupported
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
MagnitudeType achievedTol_
Tolerance achieved by the last solve() invocation.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > D_
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
int getNumIters() const override
Get the iteration count for the most recent call to solve().
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > AUTAU_
T & get(ParameterList &l, const std::string &name)
virtual ~RCGSolMgr()
Destructor.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > GY_
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
bool is_null(const std::shared_ptr< T > &p)
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > L2_
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
Teuchos::RCP< MV > U
The recycled subspace and its image.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > LUUTAU
The LU factorization of the matrix U^T A U.
Base class for Belos::SolverManager subclasses which normally can only compile with real ScalarType t...
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
A factory class for generating StatusTestOutput objects.
Implementation of the RCG (Recycling Conjugate Gradient) iterative linear solver. ...
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Beta
This class implements the RCG iteration, where a single-std::vector Krylov subspace is constructed...
int numIters_
Number of iterations taken by the last solve() invocation.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > AUTAP_
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
An implementation of StatusTestResNorm using a family of residual norms.
int scale(const ScalarType alpha)
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > UTAU_
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > rTz_old_
static const double convTol
Default convergence tolerance.
Belos::StatusTest class for specifying a maximum number of iterations.
Teuchos::RCP< Teuchos::ParameterList > params_
static std::string name()
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > AU1TU1_
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve...
A factory class for generating StatusTestOutput objects.
RCGSolMgrLinearProblemFailure(const std::string &what_arg)
RCGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i.e.
MagnitudeType convtol_
Convergence tolerance (read from parameter list).
Traits class which defines basic operations on multivectors.
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
Belos::StatusTest for logically combining several status tests.
bool isParameter(const std::string &name) const
Structure to contain pointers to RCGIter state variables.
Belos concrete class for performing the RCG iteration.
MultiVecTraits< ScalarType, MV > MVT
int maxIters_
Maximum iteration count (read from parameter list).
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
bool existU
Flag to indicate the recycle space should be used.
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Return a reference to the linear problem being solved by this solver manager.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Pure virtual base class which describes the basic interface for a solver manager. ...
Teuchos::RCP< MV > z
The current preconditioned residual.
Teuchos::RCP< MV > Ap
A times current search vector.
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Alpha_
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
void SYGV(const OrdinalType &itype, const char &JOBZ, const char &UPLO, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, ScalarType *B, const OrdinalType &ldb, ScalarType *W, ScalarType *WORK, const OrdinalType &lwork, OrdinalType *info) const
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
OperatorTraits< ScalarType, MV, OP > OPT
ReturnType
Whether the Belos solve converged for all linear systems.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > APTAP_
Teuchos::RCP< std::vector< int > > ipiv_
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > convTest_
Teuchos::RCP< OutputManager< ScalarType > > printer_
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Alpha
Coefficients arising in RCG iteration.
OrdinalType numCols() const
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.
void GESV(const OrdinalType &n, const OrdinalType &nrhs, ScalarType *A, const OrdinalType &lda, OrdinalType *IPIV, ScalarType *B, const OrdinalType &ldb, OrdinalType *info) const
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Beta_
void GETRS(const char &TRANS, const OrdinalType &n, const OrdinalType &nrhs, const ScalarType *A, const OrdinalType &lda, const OrdinalType *IPIV, ScalarType *B, const OrdinalType &ldb, OrdinalType *info) const
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Y_
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
int curDim
The current dimension of the reduction.
RCGSolMgrLAPACKFailure(const std::string &what_arg)
bool isType(const std::string &name) const
RCGSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > LUUTAU_
A class for extending the status testing capabilities of Belos via logical combinations.
Details::SolverManagerRequiresRealLapack< ScalarType, MV, OP, scalarTypeIsSupported > base_type
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.
Class which defines basic traits for the operator type.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > D
Teuchos::ScalarTraits< ScalarType > SCT
Parent class to all Belos exceptions.
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
Belos header file which uses auto-configuration information to include necessary C++ headers...
RCGSolMgrLAPACKFailure is thrown when a nonzero value is retuned from an LAPACK call.
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Delta
Solutions to local least-squares problems.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
OrdinalType stride() const
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Delta_
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...