42 #ifndef BELOS_RCG_SOLMGR_HPP 43 #define BELOS_RCG_SOLMGR_HPP 67 #ifdef BELOS_TEUCHOS_TIME_MONITOR 150 template<
class ScalarType,
class MV,
class OP,
151 const bool supportsScalarType =
156 Belos::Details::LapackSupportsScalar<ScalarType>::value &&
157 ! Teuchos::ScalarTraits<ScalarType>::isComplex>
179 template<
class ScalarType,
class MV,
class OP>
248 return Teuchos::tuple(timerSolve_);
322 std::string description()
const;
338 void sort(std::vector<ScalarType>& dlist,
int n, std::vector<int>& iperm);
341 void initializeStateStorage();
468 template<
class ScalarType,
class MV,
class OP>
472 template<
class ScalarType,
class MV,
class OP>
475 template<
class ScalarType,
class MV,
class OP>
478 template<
class ScalarType,
class MV,
class OP>
481 template<
class ScalarType,
class MV,
class OP>
484 template<
class ScalarType,
class MV,
class OP>
487 template<
class ScalarType,
class MV,
class OP>
490 template<
class ScalarType,
class MV,
class OP>
493 template<
class ScalarType,
class MV,
class OP>
496 template<
class ScalarType,
class MV,
class OP>
499 template<
class ScalarType,
class MV,
class OP>
503 template<
class ScalarType,
class MV,
class OP>
512 template<
class ScalarType,
class MV,
class OP>
530 template<
class ScalarType,
class MV,
class OP>
533 outputStream_ = outputStream_default_;
534 convtol_ = convtol_default_;
535 maxIters_ = maxIters_default_;
536 numBlocks_ = numBlocks_default_;
537 recycleBlocks_ = recycleBlocks_default_;
538 verbosity_ = verbosity_default_;
539 outputStyle_= outputStyle_default_;
540 outputFreq_= outputFreq_default_;
541 showMaxResNormOnly_ = showMaxResNormOnly_default_;
542 label_ = label_default_;
553 Alpha_ = Teuchos::null;
554 Beta_ = Teuchos::null;
556 Delta_ = Teuchos::null;
557 UTAU_ = Teuchos::null;
558 LUUTAU_ = Teuchos::null;
559 ipiv_ = Teuchos::null;
560 AUTAU_ = Teuchos::null;
561 rTz_old_ = Teuchos::null;
566 DeltaL2_ = Teuchos::null;
567 AU1TUDeltaL2_ = Teuchos::null;
568 AU1TAU1_ = Teuchos::null;
569 AU1TU1_ = Teuchos::null;
570 AU1TAP_ = Teuchos::null;
573 APTAP_ = Teuchos::null;
574 U1Y1_ = Teuchos::null;
575 PY2_ = Teuchos::null;
576 AUTAP_ = Teuchos::null;
577 AU1TU_ = Teuchos::null;
581 template<
class ScalarType,
class MV,
class OP>
585 if (params_ == Teuchos::null) {
594 maxIters_ = params->
get(
"Maximum Iterations",maxIters_default_);
597 params_->set(
"Maximum Iterations", maxIters_);
598 if (maxIterTest_!=Teuchos::null)
599 maxIterTest_->setMaxIters( maxIters_ );
604 numBlocks_ = params->
get(
"Num Blocks",numBlocks_default_);
606 "Belos::RCGSolMgr: \"Num Blocks\" must be strictly positive.");
609 params_->set(
"Num Blocks", numBlocks_);
614 recycleBlocks_ = params->
get(
"Num Recycled Blocks",recycleBlocks_default_);
616 "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be strictly positive.");
619 "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be less than \"Num Blocks\".");
622 params_->set(
"Num Recycled Blocks", recycleBlocks_);
627 std::string tempLabel = params->
get(
"Timer Label", label_default_);
630 if (tempLabel != label_) {
632 params_->set(
"Timer Label", label_);
633 std::string solveLabel = label_ +
": RCGSolMgr total solve time";
634 #ifdef BELOS_TEUCHOS_TIME_MONITOR 642 if (Teuchos::isParameterType<int>(*params,
"Verbosity")) {
643 verbosity_ = params->
get(
"Verbosity", verbosity_default_);
645 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,
"Verbosity");
649 params_->set(
"Verbosity", verbosity_);
650 if (printer_ != Teuchos::null)
651 printer_->setVerbosity(verbosity_);
656 if (Teuchos::isParameterType<int>(*params,
"Output Style")) {
657 outputStyle_ = params->
get(
"Output Style", outputStyle_default_);
659 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,
"Output Style");
663 params_->set(
"Output Style", outputStyle_);
664 outputTest_ = Teuchos::null;
669 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,
"Output Stream");
672 params_->set(
"Output Stream", outputStream_);
673 if (printer_ != Teuchos::null)
674 printer_->setOStream( outputStream_ );
680 outputFreq_ = params->
get(
"Output Frequency", outputFreq_default_);
684 params_->set(
"Output Frequency", outputFreq_);
685 if (outputTest_ != Teuchos::null)
686 outputTest_->setOutputFrequency( outputFreq_ );
690 if (printer_ == Teuchos::null) {
699 if (params->
isParameter(
"Convergence Tolerance")) {
700 convtol_ = params->
get(
"Convergence Tolerance",convtol_default_);
703 params_->set(
"Convergence Tolerance", convtol_);
704 if (convTest_ != Teuchos::null)
705 convTest_->setTolerance( convtol_ );
708 if (params->
isParameter(
"Show Maximum Residual Norm Only")) {
709 showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,
"Show Maximum Residual Norm Only");
712 params_->set(
"Show Maximum Residual Norm Only", showMaxResNormOnly_);
713 if (convTest_ != Teuchos::null)
714 convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
720 if (maxIterTest_ == Teuchos::null)
724 if (convTest_ == Teuchos::null)
725 convTest_ =
Teuchos::rcp(
new StatusTestResNorm_t( convtol_, 1 ) );
727 if (sTest_ == Teuchos::null)
728 sTest_ =
Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
730 if (outputTest_ == Teuchos::null) {
738 std::string solverDesc =
" Recycling CG ";
739 outputTest_->setSolverDesc( solverDesc );
743 if (timerSolve_ == Teuchos::null) {
744 std::string solveLabel = label_ +
": RCGSolMgr total solve time";
745 #ifdef BELOS_TEUCHOS_TIME_MONITOR 755 template<
class ScalarType,
class MV,
class OP>
764 pl->
set(
"Convergence Tolerance", convtol_default_,
765 "The relative residual tolerance that needs to be achieved by the\n" 766 "iterative solver in order for the linear system to be declared converged.");
767 pl->
set(
"Maximum Iterations", maxIters_default_,
768 "The maximum number of block iterations allowed for each\n" 769 "set of RHS solved.");
770 pl->
set(
"Block Size", blockSize_default_,
771 "Block Size Parameter -- currently must be 1 for RCG");
772 pl->
set(
"Num Blocks", numBlocks_default_,
773 "The length of a cycle (and this max number of search vectors kept)\n");
774 pl->
set(
"Num Recycled Blocks", recycleBlocks_default_,
775 "The number of vectors in the recycle subspace.");
776 pl->
set(
"Verbosity", verbosity_default_,
777 "What type(s) of solver information should be outputted\n" 778 "to the output stream.");
779 pl->
set(
"Output Style", outputStyle_default_,
780 "What style is used for the solver information outputted\n" 781 "to the output stream.");
782 pl->
set(
"Output Frequency", outputFreq_default_,
783 "How often convergence information should be outputted\n" 784 "to the output stream.");
785 pl->
set(
"Output Stream", outputStream_default_,
786 "A reference-counted pointer to the output stream where all\n" 787 "solver output is sent.");
788 pl->
set(
"Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
789 "When convergence information is printed, only show the maximum\n" 790 "relative residual norm when the block size is greater than one.");
791 pl->
set(
"Timer Label", label_default_,
792 "The string to use as a prefix for the timer labels.");
800 template<
class ScalarType,
class MV,
class OP>
805 if (rhsMV == Teuchos::null) {
813 "Belos::RCGSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
816 if (P_ == Teuchos::null) {
817 P_ = MVT::Clone( *rhsMV, numBlocks_+2 );
821 if (MVT::GetNumberVecs(*P_) < numBlocks_+2) {
823 P_ = MVT::Clone( *tmp, numBlocks_+2 );
828 if (Ap_ == Teuchos::null)
829 Ap_ = MVT::Clone( *rhsMV, 1 );
832 if (r_ == Teuchos::null)
833 r_ = MVT::Clone( *rhsMV, 1 );
836 if (z_ == Teuchos::null)
837 z_ = MVT::Clone( *rhsMV, 1 );
840 if (U_ == Teuchos::null) {
841 U_ = MVT::Clone( *rhsMV, recycleBlocks_ );
845 if (MVT::GetNumberVecs(*U_) < recycleBlocks_) {
847 U_ = MVT::Clone( *tmp, recycleBlocks_ );
852 if (AU_ == Teuchos::null) {
853 AU_ = MVT::Clone( *rhsMV, recycleBlocks_ );
857 if (MVT::GetNumberVecs(*AU_) < recycleBlocks_) {
859 AU_ = MVT::Clone( *tmp, recycleBlocks_ );
864 if (U1_ == Teuchos::null) {
865 U1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
869 if (MVT::GetNumberVecs(*U1_) < recycleBlocks_) {
871 U1_ = MVT::Clone( *tmp, recycleBlocks_ );
876 if (Alpha_ == Teuchos::null)
879 if ( (Alpha_->numRows() != numBlocks_) || (Alpha_->numCols() != 1) )
880 Alpha_->reshape( numBlocks_, 1 );
884 if (Beta_ == Teuchos::null)
887 if ( (Beta_->numRows() != (numBlocks_+1)) || (Beta_->numCols() != 1) )
888 Beta_->reshape( numBlocks_ + 1, 1 );
892 if (D_ == Teuchos::null)
895 if ( (D_->numRows() != numBlocks_) || (D_->numCols() != 1) )
896 D_->reshape( numBlocks_, 1 );
900 if (Delta_ == Teuchos::null)
903 if ( (Delta_->numRows() != recycleBlocks_) || (Delta_->numCols() != (numBlocks_ + 1)) )
904 Delta_->reshape( recycleBlocks_, numBlocks_ + 1 );
908 if (UTAU_ == Teuchos::null)
911 if ( (UTAU_->numRows() != recycleBlocks_) || (UTAU_->numCols() != recycleBlocks_) )
912 UTAU_->reshape( recycleBlocks_, recycleBlocks_ );
916 if (LUUTAU_ == Teuchos::null)
919 if ( (LUUTAU_->numRows() != recycleBlocks_) || (LUUTAU_->numCols() != recycleBlocks_) )
920 LUUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
924 if (ipiv_ == Teuchos::null)
925 ipiv_ =
Teuchos::rcp(
new std::vector<int>(recycleBlocks_) );
927 if ( (
int)ipiv_->size() != recycleBlocks_ )
928 ipiv_->resize(recycleBlocks_);
932 if (AUTAU_ == Teuchos::null)
935 if ( (AUTAU_->numRows() != recycleBlocks_) || (AUTAU_->numCols() != recycleBlocks_) )
936 AUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
940 if (rTz_old_ == Teuchos::null)
943 if ( (rTz_old_->numRows() != 1) || (rTz_old_->numCols() != 1) )
944 rTz_old_->reshape( 1, 1 );
948 if (F_ == Teuchos::null)
951 if ( (F_->numRows() != (numBlocks_+recycleBlocks_)) || (F_->numCols() != numBlocks_+recycleBlocks_) )
952 F_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
956 if (G_ == Teuchos::null)
959 if ( (G_->numRows() != (numBlocks_+recycleBlocks_)) || (G_->numCols() != numBlocks_+recycleBlocks_) )
960 G_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
964 if (Y_ == Teuchos::null)
967 if ( (Y_->numRows() != (numBlocks_+recycleBlocks_)) || (Y_->numCols() != numBlocks_+recycleBlocks_) )
968 Y_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
972 if (L2_ == Teuchos::null)
975 if ( (L2_->numRows() != (numBlocks_+1)) || (L2_->numCols() != numBlocks_) )
976 L2_->reshape( numBlocks_+1, numBlocks_ );
980 if (DeltaL2_ == Teuchos::null)
983 if ( (DeltaL2_->numRows() != (recycleBlocks_)) || (DeltaL2_->numCols() != (numBlocks_) ) )
984 DeltaL2_->reshape( recycleBlocks_, numBlocks_ );
988 if (AU1TUDeltaL2_ == Teuchos::null)
991 if ( (AU1TUDeltaL2_->numRows() != (recycleBlocks_)) || (AU1TUDeltaL2_->numCols() != (numBlocks_) ) )
992 AU1TUDeltaL2_->reshape( recycleBlocks_, numBlocks_ );
996 if (AU1TAU1_ == Teuchos::null)
999 if ( (AU1TAU1_->numRows() != (recycleBlocks_)) || (AU1TAU1_->numCols() != (recycleBlocks_) ) )
1000 AU1TAU1_->reshape( recycleBlocks_, recycleBlocks_ );
1004 if (GY_ == Teuchos::null)
1007 if ( (GY_->numRows() != (numBlocks_ + recycleBlocks_)) || (GY_->numCols() != (recycleBlocks_) ) )
1008 GY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
1012 if (AU1TU1_ == Teuchos::null)
1015 if ( (AU1TU1_->numRows() != (recycleBlocks_)) || (AU1TU1_->numCols() != (recycleBlocks_) ) )
1016 AU1TU1_->reshape( recycleBlocks_, recycleBlocks_ );
1020 if (FY_ == Teuchos::null)
1023 if ( (FY_->numRows() != (numBlocks_ + recycleBlocks_)) || (FY_->numCols() != (recycleBlocks_) ) )
1024 FY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
1028 if (AU1TAP_ == Teuchos::null)
1031 if ( (AU1TAP_->numRows() != (recycleBlocks_)) || (AU1TAP_->numCols() != (numBlocks_) ) )
1032 AU1TAP_->reshape( recycleBlocks_, numBlocks_ );
1036 if (APTAP_ == Teuchos::null)
1039 if ( (APTAP_->numRows() != (numBlocks_)) || (APTAP_->numCols() != (numBlocks_) ) )
1040 APTAP_->reshape( numBlocks_, numBlocks_ );
1044 if (U1Y1_ == Teuchos::null) {
1045 U1Y1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
1049 if (MVT::GetNumberVecs(*U1Y1_) < recycleBlocks_) {
1051 U1Y1_ = MVT::Clone( *tmp, recycleBlocks_ );
1056 if (PY2_ == Teuchos::null) {
1057 PY2_ = MVT::Clone( *rhsMV, recycleBlocks_ );
1061 if (MVT::GetNumberVecs(*PY2_) < recycleBlocks_) {
1063 PY2_ = MVT::Clone( *tmp, recycleBlocks_ );
1068 if (AUTAP_ == Teuchos::null)
1071 if ( (AUTAP_->numRows() != (recycleBlocks_)) || (AUTAP_->numCols() != (numBlocks_) ) )
1072 AUTAP_->reshape( recycleBlocks_, numBlocks_ );
1076 if (AU1TU_ == Teuchos::null)
1079 if ( (AU1TU_->numRows() != (recycleBlocks_)) || (AU1TU_->numCols() != (recycleBlocks_) ) )
1080 AU1TU_->reshape( recycleBlocks_, recycleBlocks_ );
1087 template<
class ScalarType,
class MV,
class OP>
1091 std::vector<int> index(recycleBlocks_);
1102 setParameters(Teuchos::parameterList(*getValidParameters()));
1106 "Belos::RCGSolMgr::solve(): Linear problem is not a valid object.");
1108 "Belos::RCGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
1111 "Belos::RCGSolMgr::solve(): RCG does not support split preconditioning, only set left or right preconditioner.");
1115 if (problem_->getLeftPrec() != Teuchos::null) {
1116 precObj = Teuchos::rcp_const_cast<OP>(problem_->getLeftPrec());
1118 else if (problem_->getRightPrec() != Teuchos::null) {
1119 precObj = Teuchos::rcp_const_cast<OP>(problem_->getRightPrec());
1123 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1124 std::vector<int> currIdx(1);
1128 problem_->setLSIndex( currIdx );
1131 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1132 if (numBlocks_ > dim) {
1133 numBlocks_ = Teuchos::asSafe<int>(dim);
1134 params_->set(
"Num Blocks", numBlocks_);
1136 "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
1137 " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
1141 initializeStateStorage();
1145 plist.
set(
"Num Blocks",numBlocks_);
1146 plist.
set(
"Recycled Blocks",recycleBlocks_);
1149 outputTest_->reset();
1152 bool isConverged =
true;
1156 index.resize(recycleBlocks_);
1157 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1159 index.resize(recycleBlocks_);
1160 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1163 problem_->applyOp( *Utmp, *AUtmp );
1166 MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
1169 if ( precObj != Teuchos::null ) {
1170 index.resize(recycleBlocks_);
1171 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1172 index.resize(recycleBlocks_);
1173 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1175 OPT::Apply( *precObj, *AUtmp, *PCAU );
1176 MVT::MvTransMv( one, *AUtmp, *PCAU, AUTAUtmp );
1178 MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
1190 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1194 while ( numRHS2Solve > 0 ) {
1197 if (printer_->isVerbosity(
Debug ) ) {
1198 if (existU_) printer_->print(
Debug,
"Using recycle space generated from previous call to solve()." );
1199 else printer_->print(
Debug,
"No recycle space exists." );
1203 rcg_iter->resetNumIters();
1206 rcg_iter->setSize( recycleBlocks_, numBlocks_ );
1209 outputTest_->resetNumCalls();
1218 problem_->computeCurrResVec( &*r_ );
1224 index.resize(recycleBlocks_);
1225 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1227 MVT::MvTransMv( one, *Utmp, *r_, Utr );
1230 LUUTAUtmp.
assign(UTAUtmp);
1234 "Belos::RCGSolMgr::solve(): LAPACK GESV failed to compute a solution.");
1237 MVT::MvTimesMatAddMv( one, *Utmp, Utr, one, *problem_->getCurrLHSVec() );
1240 index.resize(recycleBlocks_);
1241 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1243 MVT::MvTimesMatAddMv( -one, *AUtmp, Utr, one, *r_ );
1246 if ( precObj != Teuchos::null ) {
1247 OPT::Apply( *precObj, *r_, *z_ );
1253 MVT::MvTransMv( one, *r_, *z_, *rTz_old_ );
1258 index.resize(recycleBlocks_);
1259 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1261 MVT::MvTransMv( one, *AUtmp, *z_, mu );
1267 "Belos::RCGSolMgr::solve(): LAPACK GETRS failed to compute a solution.");
1272 MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
1273 MVT::MvTimesMatAddMv( -one, *U_, mu, one, *Ptmp );
1279 MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
1286 index.resize( numBlocks_+1 );
1287 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1288 newstate.
P = MVT::CloneViewNonConst( *P_, index );
1289 index.resize( recycleBlocks_ );
1290 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1291 newstate.
U = MVT::CloneViewNonConst( *U_, index );
1292 index.resize( recycleBlocks_ );
1293 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1294 newstate.
AU = MVT::CloneViewNonConst( *AU_, index );
1305 newstate.
existU = existU_;
1306 newstate.
ipiv = ipiv_;
1309 rcg_iter->initialize(newstate);
1315 rcg_iter->iterate();
1322 if ( convTest_->getStatus() ==
Passed ) {
1331 else if ( maxIterTest_->getStatus() ==
Passed ) {
1333 isConverged =
false;
1341 else if ( rcg_iter->getCurSubspaceDim() == rcg_iter->getMaxSubspaceDim() ) {
1346 if (recycleBlocks_ > 0) {
1357 for (
int ii=0;ii<numBlocks_;ii++) {
1358 Gtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
1360 Gtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1361 Gtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1363 Ftmp(ii,ii) = Dtmp(ii,0);
1368 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1371 index.resize( numBlocks_ );
1372 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
1374 index.resize( recycleBlocks_ );
1375 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1377 MVT::MvTimesMatAddMv( one, *Ptmp, Ytmp, zero, *U1tmp );
1397 ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
1398 for (
int ii=0; ii<recycleBlocks_; ++ii) {
1399 AU1TAPtmp(ii,0) = Ytmp(numBlocks_-1,ii) * alphatmp;
1407 lastBeta = numBlocks_-1;
1416 AU1TAPtmp.
scale(Dtmp(0,0));
1422 for (
int ii=0; ii<numBlocks_; ii++) {
1423 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
1425 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1426 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1441 for(
int ii=0;ii<numBlocks_;ii++) {
1442 F22(ii,ii) = Dtmp(ii,0);
1455 for (
int ii=0;ii<recycleBlocks_;++ii)
1456 for (
int jj=0;jj<numBlocks_;++jj)
1457 G21(jj,ii) = G12(ii,jj);
1462 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1465 index.resize( numBlocks_ );
1466 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
1468 index.resize( recycleBlocks_ );
1469 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1472 index.resize( recycleBlocks_ );
1473 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1475 index.resize( recycleBlocks_ );
1476 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1479 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1480 MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
1481 MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
1493 ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
1494 for (
int ii=0; ii<recycleBlocks_; ++ii) {
1495 AU1TAPtmp(ii,0) = Ytmp(numBlocks_+recycleBlocks_-1,ii) * alphatmp;
1505 lastp = numBlocks_+1;
1506 lastBeta = numBlocks_;
1517 for (
int ii=0; ii<numBlocks_; ii++) {
1518 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
1520 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1521 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1527 for(
int ii=0;ii<numBlocks_;ii++) {
1528 L2tmp(ii,ii) = 1./Alphatmp(ii,0);
1529 L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
1550 for(
int ii=0;ii<numBlocks_;ii++) {
1551 F22(ii,ii) = Dtmp(ii,0);
1564 for (
int ii=0;ii<recycleBlocks_;++ii)
1565 for (
int jj=0;jj<numBlocks_;++jj)
1566 G21(jj,ii) = G12(ii,jj);
1571 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1574 index.resize( recycleBlocks_ );
1575 for (
int ii=0; ii<(recycleBlocks_); ++ii) { index[ii] = ii; }
1577 index.resize( numBlocks_ );
1578 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
1580 index.resize( recycleBlocks_ );
1581 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1584 index.resize( recycleBlocks_ );
1585 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1588 index.resize( recycleBlocks_ );
1589 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1591 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1592 MVT::MvTimesMatAddMv( one, *Utmp, Y1, zero, *UY1tmp );
1593 MVT::MvAddMv(one,*UY1tmp, one, *PY2tmp, *U1tmp);
1611 AU1TUtmp.
assign(UTAUtmp);
1614 dold = Dtmp(numBlocks_-1,0);
1621 lastBeta = numBlocks_-1;
1628 for (
int ii=0; ii<numBlocks_; ii++) {
1629 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
1631 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1632 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1637 for(
int ii=0;ii<numBlocks_;ii++) {
1638 L2tmp(ii,ii) = 1./Alphatmp(ii,0);
1639 L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
1655 ScalarType val = dold * (-Betatmp(0,0)/Alphatmp(0,0));
1656 for(
int ii=0;ii<recycleBlocks_;ii++) {
1657 AU1TAPtmp(ii,0) += Y2(numBlocks_-1,ii)*val;
1663 AU1TUtmp.
assign(Y1TAU1TU);
1673 for(
int ii=0;ii<numBlocks_;ii++) {
1674 F22(ii,ii) = Dtmp(ii,0);
1687 for (
int ii=0;ii<recycleBlocks_;++ii)
1688 for (
int jj=0;jj<numBlocks_;++jj)
1689 G21(jj,ii) = G12(ii,jj);
1694 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1697 index.resize( numBlocks_ );
1698 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
1700 index.resize( recycleBlocks_ );
1701 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1703 index.resize( recycleBlocks_ );
1704 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1706 index.resize( recycleBlocks_ );
1707 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1709 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1710 MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
1711 MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
1726 dold = Dtmp(numBlocks_-1,0);
1729 lastp = numBlocks_+1;
1730 lastBeta = numBlocks_;
1740 index[0] = lastp-1; index[1] = lastp;
1742 index[0] = 0; index[1] = 1;
1743 MVT::SetBlock(*Ptmp2,index,*P_);
1746 (*Beta_)(0,0) = (*Beta_)(lastBeta,0);
1756 newstate.
P = Teuchos::null;
1757 index.resize( numBlocks_+1 );
1758 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii+1; }
1759 newstate.
P = MVT::CloneViewNonConst( *P_, index );
1761 newstate.
Beta = Teuchos::null;
1764 newstate.
Delta = Teuchos::null;
1770 rcg_iter->initialize(newstate);
1784 "Belos::RCGSolMgr::solve(): Invalid return from RCGIter::iterate().");
1787 catch (
const std::exception &e) {
1788 printer_->stream(
Errors) <<
"Error! Caught std::exception in RCGIter::iterate() at iteration " 1789 << rcg_iter->getNumIters() << std::endl
1790 << e.what() << std::endl;
1796 problem_->setCurrLS();
1800 if ( numRHS2Solve > 0 ) {
1803 problem_->setLSIndex( currIdx );
1806 currIdx.resize( numRHS2Solve );
1812 index.resize(recycleBlocks_);
1813 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1814 MVT::SetBlock(*U1_,index,*U_);
1817 if (numRHS2Solve > 0) {
1819 newstate.
P = Teuchos::null;
1820 newstate.
Ap = Teuchos::null;
1821 newstate.
r = Teuchos::null;
1822 newstate.
z = Teuchos::null;
1823 newstate.
U = Teuchos::null;
1824 newstate.
AU = Teuchos::null;
1825 newstate.
Alpha = Teuchos::null;
1826 newstate.
Beta = Teuchos::null;
1827 newstate.
D = Teuchos::null;
1828 newstate.
Delta = Teuchos::null;
1829 newstate.
LUUTAU = Teuchos::null;
1830 newstate.
ipiv = Teuchos::null;
1831 newstate.
rTz_old = Teuchos::null;
1834 index.resize(recycleBlocks_);
1835 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1837 index.resize(recycleBlocks_);
1838 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1841 problem_->applyOp( *Utmp, *AUtmp );
1844 MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
1847 if ( precObj != Teuchos::null ) {
1848 index.resize(recycleBlocks_);
1849 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1850 index.resize(recycleBlocks_);
1851 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1853 OPT::Apply( *precObj, *AUtmp, *LeftPCAU );
1854 MVT::MvTransMv( one, *AUtmp, *LeftPCAU, AUTAUtmp );
1856 MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
1869 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1878 numIters_ = maxIterTest_->getNumIters();
1882 using Teuchos::rcp_dynamic_cast;
1885 const std::vector<MagnitudeType>* pTestValues =
1886 rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
1889 "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() " 1890 "method returned NULL. Please report this bug to the Belos developers.");
1893 "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() " 1894 "method returned a vector of length zero. Please report this bug to the " 1895 "Belos developers.");
1900 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1910 template<
class ScalarType,
class MV,
class OP>
1921 std::vector<MagnitudeType> w(
n);
1924 std::vector<int> iperm(
n);
1931 std::vector<ScalarType> work(1);
1939 lapack.
SYGV(itype, jobz, uplo,
n, G2.
values(), G2.
stride(), F2.
values(), F2.
stride(), &w[0], &work[0], lwork, &info);
1941 "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to query optimal work size.");
1942 lwork = (int)work[0];
1944 lapack.
SYGV(itype, jobz, uplo,
n, G2.
values(), G2.
stride(), F2.
values(), F2.
stride(), &w[0], &work[0], lwork, &info);
1946 "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to compute eigensolutions.");
1950 this->sort(w,
n,iperm);
1953 for(
int i=0; i<recycleBlocks_; i++ ) {
1954 for(
int j=0; j<
n; j++ ) {
1955 Y(j,i) = G2(j,iperm[i]);
1962 template<
class ScalarType,
class MV,
class OP>
1965 int l, r, j, i, flag;
1994 if (dlist[j] > dlist[j - 1]) j = j + 1;
1996 if (dlist[j - 1] > dK) {
1997 dlist[i - 1] = dlist[j - 1];
1998 iperm[i - 1] = iperm[j - 1];
2011 dlist[r] = dlist[0];
2012 iperm[r] = iperm[0];
2027 template<
class ScalarType,
class MV,
class OP>
2030 std::ostringstream oss;
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
void GESV(const OrdinalType n, const OrdinalType nrhs, ScalarType *A, const OrdinalType lda, OrdinalType *IPIV, ScalarType *B, const OrdinalType ldb, OrdinalType *info) const
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.
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
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem)
Set the linear problem that needs to be solved.
Teuchos::RCP< Teuchos::Time > timerSolve_
ScalarType * values() const
bool is_null(const boost::shared_ptr< T > &p)
RCGSolMgrRecyclingFailure is thrown when any problem occurs in using/creating the recycling subspace...
static const bool scalarTypeIsSupported
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
MagnitudeType achievedTol_
Tolerance achieved by the last solve() invocation.
static const int blockSize_default_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > D_
static const int outputStyle_default_
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
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)
static const bool showMaxResNormOnly_default_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > UTAU_
int getNumIters() const
Get the iteration count for the most recent call to solve().
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > rTz_old_
Belos::StatusTest class for specifying a maximum number of iterations.
static const int numBlocks_default_
Teuchos::RCP< Teuchos::ParameterList > params_
static std::string name()
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > AU1TU1_
A factory class for generating StatusTestOutput objects.
Iterated Modified Gram-Schmidt (IMGS) implementation of the Belos::OrthoManager class.
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.
static const int recycleBlocks_default_
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
Structure to contain pointers to RCGIter state variables.
Belos concrete class for performing the RCG iteration.
MultiVecTraits< ScalarType, MV > MVT
Classical Gram-Schmidt (with DGKS correction) implementation of the Belos::OrthoManager class...
int maxIters_
Maximum iteration count (read from parameter list).
static const int outputFreq_default_
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
A Belos::StatusTest class for specifying a maximum number of iterations.
static const MagnitudeType convtol_default_
ResetType
How to reset the solver.
bool existU
Flag to indicate the recycle space should be used.
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.
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
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())
static const Teuchos::RCP< std::ostream > outputStream_default_
MagnitudeType achievedTol() const
Tolerance achieved by the last solve() invocation.
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
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_
bool isLOADetected() const
Return whether a loss of accuracy was detected by this solver during the most current solve...
Iterated Classical Gram-Schmidt (ICGS) implementation of the Belos::OrthoManager class.
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > convTest_
Teuchos::RCP< OutputManager< ScalarType > > printer_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Alpha
Coefficients arising in RCG iteration.
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.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Beta_
static const std::string label_default_
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)
void reset(const ResetType type)
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
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.
bool isParameter(const std::string &name) const
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
Details::SolverManagerRequiresRealLapack< ScalarType, MV, OP, scalarTypeIsSupported > base_type
Class which defines basic traits for the operator type.
OrdinalType stride() const
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > D
Teuchos::ScalarTraits< ScalarType > SCT
Parent class to all Belos exceptions.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Return a reference to the linear problem being solved by this solver manager.
static const int verbosity_default_
OrdinalType numCols() const
RCGSolMgrRecyclingFailure(const std::string &what_arg)
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.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Delta_
static const int maxIters_default_