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>
145 static const bool scalarTypeIsSupported =
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";
371 MagnitudeType convtol_;
377 MagnitudeType achievedTol_;
385 int numBlocks_, recycleBlocks_;
386 bool showMaxResNormOnly_;
387 int verbosity_, outputStyle_, outputFreq_;
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_;
512 Alpha_ = Teuchos::null;
513 Beta_ = Teuchos::null;
515 Delta_ = Teuchos::null;
516 UTAU_ = Teuchos::null;
517 LUUTAU_ = Teuchos::null;
518 ipiv_ = Teuchos::null;
519 AUTAU_ = Teuchos::null;
520 rTz_old_ = Teuchos::null;
525 DeltaL2_ = Teuchos::null;
526 AU1TUDeltaL2_ = Teuchos::null;
527 AU1TAU1_ = Teuchos::null;
528 AU1TU1_ = Teuchos::null;
529 AU1TAP_ = Teuchos::null;
532 APTAP_ = Teuchos::null;
533 U1Y1_ = Teuchos::null;
534 PY2_ = Teuchos::null;
535 AUTAP_ = Teuchos::null;
536 AU1TU_ = Teuchos::null;
540 template<
class ScalarType,
class MV,
class OP>
544 if (params_ == Teuchos::null) {
553 maxIters_ = params->
get(
"Maximum Iterations",maxIters_default_);
556 params_->set(
"Maximum Iterations", maxIters_);
557 if (maxIterTest_!=Teuchos::null)
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_);
609 if (printer_ != Teuchos::null)
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_);
623 outputTest_ = Teuchos::null;
628 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,
"Output Stream");
631 params_->set(
"Output Stream", outputStream_);
632 if (printer_ != Teuchos::null)
633 printer_->setOStream( outputStream_ );
639 outputFreq_ = params->
get(
"Output Frequency", outputFreq_default_);
643 params_->set(
"Output Frequency", outputFreq_);
644 if (outputTest_ != Teuchos::null)
645 outputTest_->setOutputFrequency( outputFreq_ );
649 if (printer_ == Teuchos::null) {
658 if (params->
isParameter(
"Convergence Tolerance")) {
659 if (params->
isType<MagnitudeType> (
"Convergence Tolerance")) {
660 convtol_ = params->
get (
"Convergence Tolerance",
668 params_->set(
"Convergence Tolerance", convtol_);
669 if (convTest_ != Teuchos::null)
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_);
678 if (convTest_ != Teuchos::null)
679 convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
685 if (maxIterTest_ == Teuchos::null)
689 if (convTest_ == Teuchos::null)
690 convTest_ =
Teuchos::rcp(
new StatusTestResNorm_t( convtol_, 1 ) );
692 if (sTest_ == Teuchos::null)
693 sTest_ =
Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
695 if (outputTest_ == Teuchos::null) {
703 std::string solverDesc =
" Recycling CG ";
704 outputTest_->setSolverDesc( solverDesc );
708 if (timerSolve_ == Teuchos::null) {
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>
769 if (rhsMV == Teuchos::null) {
777 "Belos::RCGSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
780 if (P_ == Teuchos::null) {
781 P_ = MVT::Clone( *rhsMV, numBlocks_+2 );
785 if (MVT::GetNumberVecs(*P_) < numBlocks_+2) {
787 P_ = MVT::Clone( *tmp, numBlocks_+2 );
792 if (Ap_ == Teuchos::null)
793 Ap_ = MVT::Clone( *rhsMV, 1 );
796 if (r_ == Teuchos::null)
797 r_ = MVT::Clone( *rhsMV, 1 );
800 if (z_ == Teuchos::null)
801 z_ = MVT::Clone( *rhsMV, 1 );
804 if (U_ == Teuchos::null) {
805 U_ = MVT::Clone( *rhsMV, recycleBlocks_ );
809 if (MVT::GetNumberVecs(*U_) < recycleBlocks_) {
811 U_ = MVT::Clone( *tmp, recycleBlocks_ );
816 if (AU_ == Teuchos::null) {
817 AU_ = MVT::Clone( *rhsMV, recycleBlocks_ );
821 if (MVT::GetNumberVecs(*AU_) < recycleBlocks_) {
823 AU_ = MVT::Clone( *tmp, recycleBlocks_ );
828 if (U1_ == Teuchos::null) {
829 U1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
833 if (MVT::GetNumberVecs(*U1_) < recycleBlocks_) {
835 U1_ = MVT::Clone( *tmp, recycleBlocks_ );
840 if (Alpha_ == Teuchos::null)
843 if ( (Alpha_->numRows() != numBlocks_) || (Alpha_->numCols() != 1) )
844 Alpha_->reshape( numBlocks_, 1 );
848 if (Beta_ == Teuchos::null)
851 if ( (Beta_->numRows() != (numBlocks_+1)) || (Beta_->numCols() != 1) )
852 Beta_->reshape( numBlocks_ + 1, 1 );
856 if (D_ == Teuchos::null)
859 if ( (D_->numRows() != numBlocks_) || (D_->numCols() != 1) )
860 D_->reshape( numBlocks_, 1 );
864 if (Delta_ == Teuchos::null)
867 if ( (Delta_->numRows() != recycleBlocks_) || (Delta_->numCols() != (numBlocks_ + 1)) )
868 Delta_->reshape( recycleBlocks_, numBlocks_ + 1 );
872 if (UTAU_ == Teuchos::null)
875 if ( (UTAU_->numRows() != recycleBlocks_) || (UTAU_->numCols() != recycleBlocks_) )
876 UTAU_->reshape( recycleBlocks_, recycleBlocks_ );
880 if (LUUTAU_ == Teuchos::null)
883 if ( (LUUTAU_->numRows() != recycleBlocks_) || (LUUTAU_->numCols() != recycleBlocks_) )
884 LUUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
888 if (ipiv_ == Teuchos::null)
889 ipiv_ =
Teuchos::rcp(
new std::vector<int>(recycleBlocks_) );
891 if ( (
int)ipiv_->size() != recycleBlocks_ )
892 ipiv_->resize(recycleBlocks_);
896 if (AUTAU_ == Teuchos::null)
899 if ( (AUTAU_->numRows() != recycleBlocks_) || (AUTAU_->numCols() != recycleBlocks_) )
900 AUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
904 if (rTz_old_ == Teuchos::null)
907 if ( (rTz_old_->numRows() != 1) || (rTz_old_->numCols() != 1) )
908 rTz_old_->reshape( 1, 1 );
912 if (F_ == Teuchos::null)
915 if ( (F_->numRows() != (numBlocks_+recycleBlocks_)) || (F_->numCols() != numBlocks_+recycleBlocks_) )
916 F_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
920 if (G_ == Teuchos::null)
923 if ( (G_->numRows() != (numBlocks_+recycleBlocks_)) || (G_->numCols() != numBlocks_+recycleBlocks_) )
924 G_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
928 if (Y_ == Teuchos::null)
931 if ( (Y_->numRows() != (numBlocks_+recycleBlocks_)) || (Y_->numCols() != numBlocks_+recycleBlocks_) )
932 Y_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
936 if (L2_ == Teuchos::null)
939 if ( (L2_->numRows() != (numBlocks_+1)) || (L2_->numCols() != numBlocks_) )
940 L2_->reshape( numBlocks_+1, numBlocks_ );
944 if (DeltaL2_ == Teuchos::null)
947 if ( (DeltaL2_->numRows() != (recycleBlocks_)) || (DeltaL2_->numCols() != (numBlocks_) ) )
948 DeltaL2_->reshape( recycleBlocks_, numBlocks_ );
952 if (AU1TUDeltaL2_ == Teuchos::null)
955 if ( (AU1TUDeltaL2_->numRows() != (recycleBlocks_)) || (AU1TUDeltaL2_->numCols() != (numBlocks_) ) )
956 AU1TUDeltaL2_->reshape( recycleBlocks_, numBlocks_ );
960 if (AU1TAU1_ == Teuchos::null)
963 if ( (AU1TAU1_->numRows() != (recycleBlocks_)) || (AU1TAU1_->numCols() != (recycleBlocks_) ) )
964 AU1TAU1_->reshape( recycleBlocks_, recycleBlocks_ );
968 if (GY_ == Teuchos::null)
971 if ( (GY_->numRows() != (numBlocks_ + recycleBlocks_)) || (GY_->numCols() != (recycleBlocks_) ) )
972 GY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
976 if (AU1TU1_ == Teuchos::null)
979 if ( (AU1TU1_->numRows() != (recycleBlocks_)) || (AU1TU1_->numCols() != (recycleBlocks_) ) )
980 AU1TU1_->reshape( recycleBlocks_, recycleBlocks_ );
984 if (FY_ == Teuchos::null)
987 if ( (FY_->numRows() != (numBlocks_ + recycleBlocks_)) || (FY_->numCols() != (recycleBlocks_) ) )
988 FY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
992 if (AU1TAP_ == Teuchos::null)
995 if ( (AU1TAP_->numRows() != (recycleBlocks_)) || (AU1TAP_->numCols() != (numBlocks_) ) )
996 AU1TAP_->reshape( recycleBlocks_, numBlocks_ );
1000 if (APTAP_ == Teuchos::null)
1003 if ( (APTAP_->numRows() != (numBlocks_)) || (APTAP_->numCols() != (numBlocks_) ) )
1004 APTAP_->reshape( numBlocks_, numBlocks_ );
1008 if (U1Y1_ == Teuchos::null) {
1009 U1Y1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
1013 if (MVT::GetNumberVecs(*U1Y1_) < recycleBlocks_) {
1015 U1Y1_ = MVT::Clone( *tmp, recycleBlocks_ );
1020 if (PY2_ == Teuchos::null) {
1021 PY2_ = MVT::Clone( *rhsMV, recycleBlocks_ );
1025 if (MVT::GetNumberVecs(*PY2_) < recycleBlocks_) {
1027 PY2_ = MVT::Clone( *tmp, recycleBlocks_ );
1032 if (AUTAP_ == Teuchos::null)
1035 if ( (AUTAP_->numRows() != (recycleBlocks_)) || (AUTAP_->numCols() != (numBlocks_) ) )
1036 AUTAP_->reshape( recycleBlocks_, numBlocks_ );
1040 if (AU1TU_ == Teuchos::null)
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.");
1079 if (problem_->getLeftPrec() != Teuchos::null) {
1080 precObj = Teuchos::rcp_const_cast<OP>(problem_->getLeftPrec());
1082 else if (problem_->getRightPrec() != Teuchos::null) {
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 );
1133 if ( precObj != Teuchos::null ) {
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_ );
1210 if ( precObj != Teuchos::null ) {
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);
1720 newstate.
P = Teuchos::null;
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 );
1725 newstate.
Beta = Teuchos::null;
1728 newstate.
Delta = Teuchos::null;
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) {
1792 newstate.
P = Teuchos::null;
1793 newstate.
Ap = Teuchos::null;
1794 newstate.
r = Teuchos::null;
1795 newstate.
z = Teuchos::null;
1796 newstate.
U = Teuchos::null;
1797 newstate.
AU = Teuchos::null;
1798 newstate.
Alpha = Teuchos::null;
1799 newstate.
Beta = Teuchos::null;
1800 newstate.
D = Teuchos::null;
1801 newstate.
Delta = Teuchos::null;
1802 newstate.
LUUTAU = Teuchos::null;
1803 newstate.
ipiv = Teuchos::null;
1804 newstate.
rTz_old = Teuchos::null;
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 );
1820 if ( precObj != Teuchos::null ) {
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>
1936 void RCGSolMgr<ScalarType,MV,OP,true>::sort(std::vector<ScalarType>& dlist,
int n, std::vector<int>& iperm)
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< 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::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.
bool is_null(const boost::shared_ptr< T > &p)
int getNumIters() const override
Get the iteration count for the most recent call to solve().
T & get(const std::string &name, T def_value)
virtual ~RCGSolMgr()
Destructor.
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)
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...
An implementation of StatusTestResNorm using a family of residual norms.
int scale(const ScalarType alpha)
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
static const double convTol
Default convergence tolerance.
Belos::StatusTest class for specifying a maximum number of iterations.
static std::string name()
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.
Traits class which defines basic operations on multivectors.
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.
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)
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...
ReturnType
Whether the Belos solve converged for all linear systems.
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
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
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
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
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)
A class for extending the status testing capabilities of Belos via logical combinations.
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
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
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...