45 #ifndef BELOS_BLOCK_GCRODR_SOLMGR_HPP 46 #define BELOS_BLOCK_GCRODR_SOLMGR_HPP 62 #include "Teuchos_LAPACK.hpp" 63 #include "Teuchos_as.hpp" 64 #ifdef BELOS_TEUCHOS_TIME_MONITOR 65 #include "Teuchos_TimeMonitor.hpp" 66 #endif // BELOS_TEUCHOS_TIME_MONITOR 125 template<
class ScalarType,
class MV,
class OP>
131 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
132 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
MagnitudeType;
133 typedef Teuchos::ScalarTraits<MagnitudeType>
MT;
134 typedef Teuchos::ScalarTraits<MagnitudeType>
SMT;
136 typedef Teuchos::SerialDenseMatrix<int,ScalarType>
SDM;
137 typedef Teuchos::SerialDenseVector<int,ScalarType>
SDV;
176 const Teuchos::RCP<Teuchos::ParameterList> &pl);
223 TEUCHOS_TEST_FOR_EXCEPTION(problem.is_null(), std::invalid_argument,
224 "Belos::BlockGCRODRSolMgr::setProblem: The input LinearProblem cannot be null.");
227 if (! problem->isProblemSet()) {
228 const bool success = problem->setProblem();
229 TEUCHOS_TEST_FOR_EXCEPTION(success, std::runtime_error,
230 "Belos::BlockGCRODRSolMgr::setProblem: Calling the input LinearProblem's setProblem() method failed. This likely means that the " 231 "LinearProblem has a missing (null) matrix A, solution vector X, or right-hand side vector B. Please set these items in the LinearProblem and try again.");
238 void setParameters(
const Teuchos::RCP<Teuchos::ParameterList> ¶ms );
317 const Teuchos::RCP<const MV>& VV,
321 void sort (std::vector<MagnitudeType>& dlist,
int n, std::vector<int>& iperm);
327 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >
problem_;
334 Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
sTest_;
348 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> >
ortho_;
400 Teuchos::RCP<SDM >
G_;
401 Teuchos::RCP<SDM >
H_;
402 Teuchos::RCP<SDM >
B_;
407 Teuchos::RCP<SDM >
F_;
426 template<
class ScalarType,
class MV,
class OP>
429 template<
class ScalarType,
class MV,
class OP>
436 template<
class ScalarType,
class MV,
class OP>
442 template<
class ScalarType,
class MV,
class OP>
445 const Teuchos::RCP<Teuchos::ParameterList> &pl ) {
450 TEUCHOS_TEST_FOR_EXCEPTION(problem.is_null(), std::invalid_argument,
451 "Belos::BlockGCRODR constructor: The solver manager's constructor needs " 452 "the linear problem argument 'problem' to be nonnull.");
462 template<
class ScalarType,
class MV,
class OP>
464 adaptiveBlockSize_ = adaptiveBlockSize_default_;
465 recycleMethod_ = recycleMethod_default_;
467 loaDetected_ =
false;
468 builtRecycleSpace_ =
false;
539 template<
class ScalarType,
class MV,
class OP>
541 std::ostringstream oss;
542 oss <<
"Belos::BlockGCRODRSolMgr<" << SCT::name() <<
", ...>";
544 oss <<
"Ortho Type='"<<orthoType_ ;
545 oss <<
", Num Blocks=" <<numBlocks_;
546 oss <<
", Block Size=" <<blockSize_;
547 oss <<
", Num Recycle Blocks=" << recycledBlocks_;
548 oss <<
", Max Restarts=" << maxRestarts_;
553 template<
class ScalarType,
class MV,
class OP>
554 Teuchos::RCP<const Teuchos::ParameterList>
556 using Teuchos::ParameterList;
557 using Teuchos::parameterList;
559 using Teuchos::rcpFromRef;
561 if (defaultParams_.is_null()) {
562 RCP<ParameterList> pl = parameterList ();
564 const MagnitudeType convTol = SMT::squareroot (SCT::magnitude (SCT::eps()));
565 const int maxRestarts = 1000;
566 const int maxIters = 5000;
567 const int blockSize = 2;
568 const int numBlocks = 100;
569 const int numRecycledBlocks = 25;
572 const int outputFreq = 1;
573 RCP<std::ostream> outputStream = rcpFromRef (std::cout);
574 const std::string impResScale (
"Norm of Preconditioned Initial Residual");
575 const std::string expResScale (
"Norm of Initial Residual");
576 const std::string timerLabel (
"Belos");
577 const std::string orthoType (
"ICGS");
578 RCP<const ParameterList> orthoParams = orthoFactory_.getDefaultParameters (orthoType);
588 pl->set (
"Convergence Tolerance", convTol,
589 "The tolerance that the solver needs to achieve in order for " 590 "the linear system(s) to be declared converged. The meaning " 591 "of this tolerance depends on the convergence test details.");
592 pl->set(
"Maximum Restarts", maxRestarts,
593 "The maximum number of restart cycles (not counting the first) " 594 "allowed for each set of right-hand sides solved.");
595 pl->set(
"Maximum Iterations", maxIters,
596 "The maximum number of iterations allowed for each set of " 597 "right-hand sides solved.");
598 pl->set(
"Block Size", blockSize,
599 "How many linear systems to solve at once.");
600 pl->set(
"Num Blocks", numBlocks,
601 "The maximum number of blocks allowed in the Krylov subspace " 602 "for each set of right-hand sides solved. This includes the " 603 "initial block. Total Krylov basis storage, not counting the " 604 "recycled subspace, is \"Num Blocks\" times \"Block Size\".");
605 pl->set(
"Num Recycled Blocks", numRecycledBlocks,
606 "The maximum number of vectors in the recycled subspace." );
607 pl->set(
"Verbosity", verbosity,
608 "What type(s) of solver information should be written " 609 "to the output stream.");
610 pl->set(
"Output Style", outputStyle,
611 "What style is used for the solver information to write " 612 "to the output stream.");
613 pl->set(
"Output Frequency", outputFreq,
614 "How often convergence information should be written " 615 "to the output stream.");
616 pl->set(
"Output Stream", outputStream,
617 "A reference-counted pointer to the output stream where all " 618 "solver output is sent.");
619 pl->set(
"Implicit Residual Scaling", impResScale,
620 "The type of scaling used in the implicit residual convergence test.");
621 pl->set(
"Explicit Residual Scaling", expResScale,
622 "The type of scaling used in the explicit residual convergence test.");
623 pl->set(
"Timer Label", timerLabel,
624 "The string to use as a prefix for the timer labels.");
625 pl->set(
"Orthogonalization", orthoType,
626 "The orthogonalization method to use. Valid options: " +
627 orthoFactory_.validNamesString());
628 pl->set (
"Orthogonalization Parameters", *orthoParams,
629 "Sublist of parameters specific to the orthogonalization method to use.");
630 pl->set(
"Orthogonalization Constant", orthoKappa,
631 "When using DGKS orthogonalization: the \"depTol\" constant, used " 632 "to determine whether another step of classical Gram-Schmidt is " 633 "necessary. Otherwise ignored. Nonpositive values are ignored.");
636 return defaultParams_;
639 template<
class ScalarType,
class MV,
class OP>
643 using Teuchos::isParameterType;
644 using Teuchos::getParameter;
646 using Teuchos::ParameterList;
647 using Teuchos::parameterList;
650 using Teuchos::rcp_dynamic_cast;
651 using Teuchos::rcpFromRef;
652 using Teuchos::Exceptions::InvalidParameter;
653 using Teuchos::Exceptions::InvalidParameterName;
654 using Teuchos::Exceptions::InvalidParameterType;
657 RCP<const ParameterList> defaultParams = getValidParameters();
659 if (params.is_null()) {
664 params_ = parameterList (*defaultParams);
676 params->validateParametersAndSetDefaults (*defaultParams);
692 const int maxRestarts = params_->get<
int> (
"Maximum Restarts");
693 TEUCHOS_TEST_FOR_EXCEPTION(maxRestarts <= 0, std::invalid_argument,
694 "Belos::BlockGCRODRSolMgr: The \"Maximum Restarts\" parameter " 695 "must be nonnegative, but you specified a negative value of " 696 << maxRestarts <<
".");
698 const int maxIters = params_->get<
int> (
"Maximum Iterations");
699 TEUCHOS_TEST_FOR_EXCEPTION(maxIters <= 0, std::invalid_argument,
700 "Belos::BlockGCRODRSolMgr: The \"Maximum Iterations\" parameter " 701 "must be positive, but you specified a nonpositive value of " 704 const int numBlocks = params_->get<
int> (
"Num Blocks");
705 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument,
706 "Belos::BlockGCRODRSolMgr: The \"Num Blocks\" parameter must be " 707 "positive, but you specified a nonpositive value of " << numBlocks
710 const int blockSize = params_->get<
int> (
"Block Size");
711 TEUCHOS_TEST_FOR_EXCEPTION(blockSize <= 0, std::invalid_argument,
712 "Belos::BlockGCRODRSolMgr: The \"Block Size\" parameter must be " 713 "positive, but you specified a nonpositive value of " << blockSize
716 const int recycledBlocks = params_->get<
int> (
"Num Recycled Blocks");
717 TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks <= 0, std::invalid_argument,
718 "Belos::BlockGCRODRSolMgr: The \"Num Recycled Blocks\" parameter must " 719 "be positive, but you specified a nonpositive value of " 720 << recycledBlocks <<
".");
721 TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks >= numBlocks,
722 std::invalid_argument,
"Belos::BlockGCRODRSolMgr: The \"Num Recycled " 723 "Blocks\" parameter must be less than the \"Num Blocks\" parameter, " 724 "but you specified \"Num Recycled Blocks\" = " << recycledBlocks
725 <<
" and \"Num Blocks\" = " << numBlocks <<
".");
729 maxRestarts_ = maxRestarts;
730 maxIters_ = maxIters;
731 numBlocks_ = numBlocks;
732 blockSize_ = blockSize;
733 recycledBlocks_ = recycledBlocks;
740 std::string tempLabel = params_->get<std::string> (
"Timer Label");
741 const bool labelChanged = (tempLabel != label_);
743 #ifdef BELOS_TEUCHOS_TIME_MONITOR 744 std::string solveLabel = label_ +
": BlockGCRODRSolMgr total solve time";
745 if (timerSolve_.is_null()) {
747 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
748 }
else if (labelChanged) {
754 Teuchos::TimeMonitor::clearCounter (solveLabel);
755 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
757 #endif // BELOS_TEUCHOS_TIME_MONITOR 763 if (params_->isParameter (
"Verbosity")) {
764 if (isParameterType<int> (*params_,
"Verbosity")) {
765 verbosity_ = params_->get<
int> (
"Verbosity");
768 verbosity_ = (int) getParameter<MsgType> (*params_,
"Verbosity");
775 if (params_->isParameter (
"Output Style")) {
776 if (isParameterType<int> (*params_,
"Output Style")) {
777 outputStyle_ = params_->get<
int> (
"Output Style");
780 outputStyle_ = (int) getParameter<OutputType> (*params_,
"Output Style");
808 if (params_->isParameter (
"Output Stream")) {
810 outputStream_ = getParameter<RCP<std::ostream> > (*params_,
"Output Stream");
812 catch (InvalidParameter&) {
813 outputStream_ = rcpFromRef (std::cout);
820 if (outputStream_.is_null()) {
821 outputStream_ = rcp (
new Teuchos::oblackholestream);
825 outputFreq_ = params_->get<
int> (
"Output Frequency");
828 if (! outputTest_.is_null()) {
829 outputTest_->setOutputFrequency (outputFreq_);
837 if (printer_.is_null()) {
840 printer_->setVerbosity (verbosity_);
841 printer_->setOStream (outputStream_);
852 if (params_->isParameter (
"Orthogonalization")) {
853 const std::string& tempOrthoType =
854 params_->get<std::string> (
"Orthogonalization");
856 if (! orthoFactory_.isValidName (tempOrthoType)) {
857 std::ostringstream os;
858 os <<
"Belos::BlockGCRODRSolMgr: Invalid orthogonalization name \"" 859 << tempOrthoType <<
"\". The following are valid options " 860 <<
"for the \"Orthogonalization\" name parameter: ";
861 orthoFactory_.printValidNames (os);
862 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument, os.str());
864 if (tempOrthoType != orthoType_) {
866 orthoType_ = tempOrthoType;
883 RCP<ParameterList> orthoParams = sublist (params,
"Orthogonalization Parameters",
true);
884 TEUCHOS_TEST_FOR_EXCEPTION(orthoParams.is_null(), std::logic_error,
885 "Failed to get orthogonalization parameters. " 886 "Please report this bug to the Belos developers.");
912 ortho_ = orthoFactory_.makeMatOrthoManager (orthoType_, null, printer_,
913 label_, orthoParams);
925 if (params_->isParameter (
"Orthogonalization Constant")) {
928 if (orthoKappa > 0) {
929 orthoKappa_ = orthoKappa;
932 if (orthoType_ ==
"DGKS" && ! ortho_.is_null()) {
938 rcp_dynamic_cast<ortho_man_type>(ortho_)->setDepTol (orthoKappa_);
948 convTol_ = params_->get<
MagnitudeType> (
"Convergence Tolerance");
949 if (! impConvTest_.is_null()) {
952 if (! expConvTest_.is_null()) {
953 expConvTest_->setTolerance (convTol_);
957 if (params_->isParameter (
"Implicit Residual Scaling")) {
958 std::string tempImpResScale =
959 getParameter<std::string> (*params_,
"Implicit Residual Scaling");
962 if (impResScale_ != tempImpResScale) {
964 impResScale_ = tempImpResScale;
966 if (! impConvTest_.is_null()) {
979 if (params_->isParameter(
"Explicit Residual Scaling")) {
980 std::string tempExpResScale =
981 getParameter<std::string> (*params_,
"Explicit Residual Scaling");
984 if (expResScale_ != tempExpResScale) {
986 expResScale_ = tempExpResScale;
988 if (! expConvTest_.is_null()) {
1006 if (maxIterTest_.is_null()) {
1009 maxIterTest_->setMaxIters (maxIters_);
1014 if (impConvTest_.is_null()) {
1015 impConvTest_ = rcp (
new StatusTestResNorm_t (convTol_));
1021 if (expConvTest_.is_null()) {
1022 expConvTest_ = rcp (
new StatusTestResNorm_t (convTol_));
1023 expConvTest_->defineResForm (StatusTestResNorm_t::Explicit,
Belos::TwoNorm);
1029 if (convTest_.is_null()) {
1030 convTest_ = rcp (
new StatusTestCombo_t (StatusTestCombo_t::SEQ,
1038 sTest_ = rcp (
new StatusTestCombo_t (StatusTestCombo_t::OR,
1039 maxIterTest_, convTest_));
1043 outputTest_ = stoFactory.
create (printer_, sTest_, outputFreq_,
1047 std::string solverDesc =
"Block GCRODR ";
1048 outputTest_->setSolverDesc (solverDesc);
1055 template<
class ScalarType,
class MV,
class OP>
1060 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1063 Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
1066 int KrylSpaDim = (numBlocks_ - 1) * blockSize_;
1069 int augSpaDim = KrylSpaDim + recycledBlocks_ + 1;
1072 int augSpaImgDim = KrylSpaDim + blockSize_ + recycledBlocks_+1;
1096 if (U_ == Teuchos::null) {
1097 U_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1101 if (MVT::GetNumberVecs(*U_) < recycledBlocks_+1) {
1102 Teuchos::RCP<const MV> tmp = U_;
1103 U_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1108 if (C_ == Teuchos::null) {
1109 C_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1113 if (MVT::GetNumberVecs(*C_) < recycledBlocks_+1) {
1114 Teuchos::RCP<const MV> tmp = C_;
1115 C_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1120 if (U1_ == Teuchos::null) {
1121 U1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1125 if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1126 Teuchos::RCP<const MV> tmp = U1_;
1127 U1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1132 if (C1_ == Teuchos::null) {
1133 C1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1137 if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1138 Teuchos::RCP<const MV> tmp = C1_;
1139 C1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1144 if (R_ == Teuchos::null){
1145 R_ = MVT::Clone( *rhsMV, blockSize_ );
1151 if (G_ == Teuchos::null){
1152 G_ = Teuchos::rcp(
new SDM( augSpaImgDim, augSpaDim ) );
1155 if ( (G_->numRows() != augSpaImgDim) || (G_->numCols() != augSpaDim) )
1157 G_->reshape( augSpaImgDim, augSpaDim );
1159 G_->putScalar(zero);
1163 if (H_ == Teuchos::null){
1164 H_ = Teuchos::rcp (
new SDM ( Teuchos::View, *G_, KrylSpaDim + blockSize_, KrylSpaDim, recycledBlocks_+1 ,recycledBlocks_+1 ) );
1168 if (F_ == Teuchos::null){
1169 F_ = Teuchos::rcp(
new SDM( recycledBlocks_+1, recycledBlocks_+1 ) );
1172 if ( (F_->numRows() != recycledBlocks_+1) || (F_->numCols() != recycledBlocks_+1) ){
1173 F_->reshape( recycledBlocks_+1, recycledBlocks_+1 );
1176 F_->putScalar(zero);
1179 if (PP_ == Teuchos::null){
1180 PP_ = Teuchos::rcp(
new SDM( augSpaImgDim, recycledBlocks_+1 ) );
1183 if ( (PP_->numRows() != augSpaImgDim) || (PP_->numCols() != recycledBlocks_+1) ){
1184 PP_->reshape( augSpaImgDim, recycledBlocks_+1 );
1189 if (HP_ == Teuchos::null)
1190 HP_ = Teuchos::rcp(
new SDM( augSpaImgDim, augSpaDim ) );
1192 if ( (HP_->numRows() != augSpaImgDim) || (HP_->numCols() != augSpaDim) ){
1193 HP_->reshape( augSpaImgDim, augSpaDim );
1198 tau_.resize(recycledBlocks_+1);
1201 work_.resize(recycledBlocks_+1);
1204 ipiv_.resize(recycledBlocks_+1);
1208 template<
class ScalarType,
class MV,
class OP>
1211 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1212 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1214 int p = block_gmres_iter->getState().curDim;
1215 std::vector<int> index(keff);
1219 SDM HH(Teuchos::Copy, *H_, p+blockSize_, p);
1220 if(recycledBlocks_ >= p + blockSize_){
1224 for (
int ii=0; ii < p; ++ii) index[ii] = ii;
1225 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1226 MVT::SetBlock(*V_, index, *Utmp);
1231 Teuchos::RCP<SDM > PPtmp = Teuchos::rcp (
new SDM ( Teuchos::View, *PP_, p, recycledBlocks_+1 ) );
1232 if(recycleMethod_ ==
"harmvecs"){
1233 keff = getHarmonicVecsKryl(p, HH, *PPtmp);
1234 printer_->stream(
Debug) <<
"keff = " << keff << std::endl;
1237 PPtmp = Teuchos::rcp (
new SDM ( Teuchos::View, *PP_, p, keff ) );
1240 for (
int ii=0; ii<keff; ++ii) index[ii] = ii;
1241 Teuchos::RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1242 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1243 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1245 for (
int ii=0; ii < p; ++ii) index[ii] = ii;
1246 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1250 MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
1253 SDM HPtmp( Teuchos::View, *HP_, p+blockSize_, keff );
1254 HPtmp.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, *H_, *PPtmp, zero );
1258 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1259 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
1262 lwork =
static_cast<int> (Teuchos::ScalarTraits<ScalarType>::magnitude (work_[0]));
1263 work_.resize(lwork);
1264 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1265 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
1269 SDM Rtmp( Teuchos::View, *F_, keff, keff );
1270 for(
int ii=0;ii<keff;ii++) {
for(
int jj=ii;jj<keff;jj++) Rtmp(ii,jj) = HPtmp(ii,jj); }
1272 lapack.UNGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1273 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _UNGQR failed to construct the Q factor.");
1277 index.resize( p+blockSize_ );
1278 for (
int ii=0; ii < (p+blockSize_); ++ii) { index[ii] = ii; }
1279 Vtmp = MVT::CloneView( *V_, index );
1280 MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
1287 ipiv_.resize(Rtmp.numRows());
1288 lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
1289 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1291 lwork = Rtmp.numRows();
1292 work_.resize(lwork);
1293 lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1296 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
1302 template<
class ScalarType,
class MV,
class OP>
1304 const MagnitudeType one = Teuchos::ScalarTraits<ScalarType>::one();
1305 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1307 std::vector<MagnitudeType> d(keff);
1308 std::vector<ScalarType> dscalar(keff);
1309 std::vector<int> index(numBlocks_+1);
1318 if(recycledBlocks_ >= keff + p){
1321 for (
int ii=0; ii < p; ++ii) { index[ii] = keff+ii; }
1322 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1323 for (
int ii=0; ii < p; ++ii) { index[ii] =ii; }
1324 MVT::SetBlock(*V_, index, *Utmp);
1331 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1332 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1334 dscalar.resize(keff);
1335 MVT::MvNorm( *Utmp, d );
1336 for (
int i=0; i<keff; ++i) {
1338 dscalar[i] = (ScalarType)d[i];
1340 MVT::MvScale( *Utmp, dscalar );
1345 Teuchos::RCP<SDM> Gtmp = Teuchos::rcp(
new SDM( Teuchos::View, *G_, p+keff, p+keff-blockSize_ ) );
1348 for (
int i=0; i<keff; ++i)
1349 (*Gtmp)(i,i) = d[i];
1356 SDM PPtmp( Teuchos::View, *PP_, p+keff-blockSize_, recycledBlocks_+1 );
1357 keff_new = getHarmonicVecsAugKryl( keff, p-blockSize_, *Gtmp, oldState.
V, PPtmp );
1364 Teuchos::RCP<MV> U1tmp;
1366 index.resize( keff );
1367 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1368 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1369 index.resize( keff_new );
1370 for (
int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
1371 U1tmp = MVT::CloneViewNonConst( *U1_, index );
1372 SDM PPtmp( Teuchos::View, *PP_, keff, keff_new );
1373 MVT::MvTimesMatAddMv( one, *Utmp, PPtmp, zero, *U1tmp );
1378 index.resize(p-blockSize_);
1379 for (
int ii=0; ii < p-blockSize_; ii++) { index[ii] = ii; }
1380 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1381 SDM PPtmp( Teuchos::View, *PP_, p-blockSize_, keff_new, keff );
1382 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
1386 SDM HPtmp( Teuchos::View, *HP_, p+keff, keff_new );
1388 SDM PPtmp( Teuchos::View, *PP_, p-blockSize_+keff, keff_new );
1389 HPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*Gtmp,PPtmp,zero);
1393 int info = 0, lwork = -1;
1394 tau_.resize(keff_new);
1395 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1396 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
1398 lwork =
static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work_[0]));
1399 work_.resize(lwork);
1400 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1401 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
1405 SDM Ftmp( Teuchos::View, *F_, keff_new, keff_new );
1406 for(
int i=0;i<keff_new;i++) {
for(
int j=i;j<keff_new;j++) Ftmp(i,j) = HPtmp(i,j); }
1408 lapack.UNGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1409 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _UNGQR failed to construct the Q factor.");
1416 Teuchos::RCP<MV> C1tmp;
1419 for (
int i=0; i < keff; i++) { index[i] = i; }
1420 Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
1421 index.resize(keff_new);
1422 for (
int i=0; i < keff_new; i++) { index[i] = i; }
1423 C1tmp = MVT::CloneViewNonConst( *C1_, index );
1424 SDM PPtmp( Teuchos::View, *HP_, keff, keff_new );
1425 MVT::MvTimesMatAddMv( one, *Ctmp, PPtmp, zero, *C1tmp );
1430 for (
int i=0; i < p; ++i) { index[i] = i; }
1431 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1432 SDM PPtmp( Teuchos::View, *HP_, p, keff_new, keff, 0 );
1433 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *C1tmp );
1442 ipiv_.resize(Ftmp.numRows());
1443 lapack.GETRF(Ftmp.numRows(),Ftmp.numCols(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&info);
1444 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1447 lwork = Ftmp.numRows();
1448 work_.resize(lwork);
1449 lapack.GETRI(Ftmp.numRows(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1450 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRI failed to compute an LU factorization.");
1453 index.resize(keff_new);
1454 for (
int i=0; i < keff_new; i++) { index[i] = i; }
1455 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1456 MVT::MvTimesMatAddMv( one, *U1tmp, Ftmp, zero, *Utmp );
1461 if (keff != keff_new) {
1463 block_gcrodr_iter->setSize( keff, numBlocks_ );
1465 SDM b1( Teuchos::View, *G_, recycledBlocks_+2, 1, 0, recycledBlocks_ );
1471 template<
class ScalarType,
class MV,
class OP>
1474 int m2 = GG.numCols();
1475 bool xtraVec =
false;
1476 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1477 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1478 std::vector<int> index;
1481 std::vector<MagnitudeType> wr(m2), wi(m2);
1484 std::vector<MagnitudeType> w(m2);
1487 SDM vr(m2,m2,
false);
1490 std::vector<int> iperm(m2);
1493 builtRecycleSpace_ =
true;
1503 SDM A_tmp( keff+m+blockSize_, keff+m );
1508 for (
int i=0; i<keff; ++i) { index[i] = i; }
1509 Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
1510 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1511 SDM A11( Teuchos::View, A_tmp, keff, keff );
1512 MVT::MvTransMv( one, *Ctmp, *Utmp, A11 );
1515 SDM A21( Teuchos::View, A_tmp, m+blockSize_, keff, keff );
1516 index.resize(m+blockSize_);
1517 for (i=0; i < m+blockSize_; i++) { index[i] = i; }
1518 Teuchos::RCP<const MV> Vp = MVT::CloneView( *VV, index );
1519 MVT::MvTransMv( one, *Vp, *Utmp, A21 );
1522 for( i=keff; i<keff+m; i++)
1526 SDM A( m2, A_tmp.numCols() );
1527 A.multiply(
Teuchos::TRANS, Teuchos::NO_TRANS, one, GG, A_tmp, zero );
1535 char balanc=
'P', jobvl=
'N', jobvr=
'V', sense=
'N';
1536 int ld = A.numRows();
1538 int ldvl = ld, ldvr = ld;
1539 int info = 0,ilo = 0,ihi = 0;
1542 std::vector<ScalarType> beta(ld);
1543 std::vector<ScalarType> work(lwork);
1544 std::vector<MagnitudeType> rwork(lwork);
1545 std::vector<MagnitudeType> lscale(ld), rscale(ld);
1546 std::vector<MagnitudeType> rconde(ld), rcondv(ld);
1547 std::vector<int> iwork(ld+6);
1552 lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0],
1553 &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
1554 &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &rwork[0],
1555 &iwork[0], bwork, &info);
1556 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK GGEVX failed to compute eigensolutions.");
1560 for( i=0; i<ld; i++ )
1561 w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] ) / std::abs( beta[i] );
1563 this->sort(w,ld,iperm);
1565 bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
1568 for( i=0; i<recycledBlocks_; i++ )
1569 for( j=0; j<ld; j++ )
1570 PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
1572 if(scalarTypeIsComplex==
false) {
1575 if (wi[iperm[ld-recycledBlocks_]] != 0.0) {
1577 for ( i=ld-recycledBlocks_; i<ld; i++ )
1578 if (wi[iperm[i]] != 0.0) countImag++;
1580 if (countImag % 2) xtraVec =
true;
1584 if (wi[iperm[ld-recycledBlocks_]] > 0.0) {
1585 for( j=0; j<ld; j++ )
1586 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
1589 for( j=0; j<ld; j++ )
1590 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
1598 return recycledBlocks_+1;
1600 return recycledBlocks_;
1604 template<
class ScalarType,
class MV,
class OP>
1606 bool xtraVec =
false;
1607 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1608 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1611 std::vector<MagnitudeType> wr(m), wi(m);
1617 std::vector<MagnitudeType> w(m);
1620 std::vector<int> iperm(m);
1626 builtRecycleSpace_ =
true;
1630 Teuchos::RCP<SDM> harmRitzMatrix = Teuchos::rcp(
new SDM( m, blockSize_));
1633 for(
int i=0; i<=blockSize_-1; i++) (*harmRitzMatrix)[blockSize_-1-i][harmRitzMatrix->numRows()-1-i] = 1;
1636 lapack.GESV(m, blockSize_, HHt.values(), HHt.stride(), &iperm[0], harmRitzMatrix->values(), harmRitzMatrix->stride(), &info);
1638 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK GESV failed to compute a solution.");
1642 Teuchos::SerialDenseMatrix<int, ScalarType> H_lbl(Teuchos::View, HH, blockSize_, blockSize_, (HH).numRows()-blockSize_, (HH).numCols()-blockSize_ );
1643 Teuchos::SerialDenseMatrix<int, ScalarType> H_lbl_t( H_lbl,
Teuchos::TRANS );
1648 Teuchos::RCP<SDM> Htemp = Teuchos::null;
1649 Htemp = Teuchos::rcp(
new SDM(H_lbl_t.numRows(), H_lbl_t.numCols()));
1650 Htemp -> multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, H_lbl_t, H_lbl, zero);
1651 H_lbl_t.assign(*Htemp);
1653 Htemp = Teuchos::rcp(
new SDM(harmRitzMatrix -> numRows(), harmRitzMatrix -> numCols()));
1654 Htemp -> multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, *harmRitzMatrix, H_lbl_t, zero);
1655 harmRitzMatrix -> assign(*Htemp);
1658 int harmColIndex, HHColIndex;
1659 Htemp = Teuchos::rcp(
new SDM(Teuchos::Copy,HH,HH.numRows()-blockSize_,HH.numCols()));
1660 for(
int i = 0; i<blockSize_; i++){
1661 harmColIndex = harmRitzMatrix -> numCols() - i -1;
1663 for(
int j=0; j<m; j++) (*Htemp)[HHColIndex][j] += (*harmRitzMatrix)[harmColIndex][j];
1665 harmRitzMatrix = Htemp;
1673 std::vector<ScalarType> work(1);
1674 std::vector<MagnitudeType> rwork(2*m);
1680 lapack.GEEV(
'N',
'V', m, harmRitzMatrix->values(), harmRitzMatrix->stride(), &wr[0], &wi[0],
1681 vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
1683 lwork = std::abs (static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work[0])));
1684 work.resize( lwork );
1686 lapack.GEEV(
'N',
'V', m, harmRitzMatrix->values(), harmRitzMatrix->stride(), &wr[0], &wi[0],
1687 vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
1689 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK GEEV failed to compute eigensolutions.");
1692 for(
int i=0; i<m; ++i ) w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] );
1694 this->sort(w, m, iperm);
1696 bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
1699 for(
int i=0; i<recycledBlocks_; ++i )
1700 for(
int j=0; j<m; j++ )
1701 PP(j,i) = vr(j,iperm[i]);
1703 if(scalarTypeIsComplex==
false) {
1706 if (wi[iperm[recycledBlocks_-1]] != 0.0) {
1708 for (
int i=0; i<recycledBlocks_; ++i )
1709 if (wi[iperm[i]] != 0.0) countImag++;
1711 if (countImag % 2) xtraVec =
true;
1715 if (wi[iperm[recycledBlocks_-1]] > 0.0) {
1716 for(
int j=0; j<m; ++j )
1717 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
1720 for(
int j=0; j<m; ++j )
1721 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
1729 printer_->stream(
Debug) <<
"Recycled " << recycledBlocks_+1 <<
" vectors" << std::endl;
1730 return recycledBlocks_+1;
1733 printer_->stream(
Debug) <<
"Recycled " << recycledBlocks_ <<
" vectors" << std::endl;
1734 return recycledBlocks_;
1739 template<
class ScalarType,
class MV,
class OP>
1741 int l, r, j, i, flag;
1768 if (dlist[j] > dlist[j - 1]) j = j + 1;
1769 if (dlist[j - 1] > dK) {
1770 dlist[i - 1] = dlist[j - 1];
1771 iperm[i - 1] = iperm[j - 1];
1784 dlist[r] = dlist[0];
1785 iperm[r] = iperm[0];
1799 template<
class ScalarType,
class MV,
class OP>
1803 using Teuchos::rcp_const_cast;
1807 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1808 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1809 std::vector<int> index(numBlocks_+1);
1825 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1826 int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1830 std::vector<int> currIdx;
1832 if ( adaptiveBlockSize_ ) {
1833 blockSize_ = numCurrRHS;
1834 currIdx.resize( numCurrRHS );
1835 for (
int i=0; i<numCurrRHS; ++i)
1836 currIdx[i] = startPtr+i;
1839 currIdx.resize( blockSize_ );
1840 for (
int i=0; i<numCurrRHS; ++i)
1841 currIdx[i] = startPtr+i;
1842 for (
int i=numCurrRHS; i<blockSize_; ++i)
1847 problem_->setLSIndex( currIdx );
1853 loaDetected_ =
false;
1856 bool isConverged =
true;
1859 initializeStateStorage();
1862 Teuchos::ParameterList plist;
1864 while (numRHS2Solve > 0){
1874 printer_->stream(
Debug) <<
" Now solving RHS index " << currIdx[0] <<
" using recycled subspace of dimension " << keff << std::endl << std::endl;
1878 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1879 RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1880 RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1881 problem_->apply( *Utmp, *Ctmp );
1883 RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1887 SDM Ftmp( Teuchos::View, *F_, keff, keff );
1888 int rank = ortho_->normalize(*Ctmp, rcp(&Ftmp,
false));
1890 TEUCHOS_TEST_FOR_EXCEPTION(rank != keff,
BlockGCRODRSolMgrOrthoFailure,
"Belos::BlockGCRODRSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
1895 ipiv_.resize(Ftmp.numRows());
1896 lapack.GETRF(Ftmp.numRows(),Ftmp.numCols(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&info);
1897 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1900 int lwork = Ftmp.numRows();
1901 work_.resize(lwork);
1902 lapack.GETRI(Ftmp.numRows(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1903 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
1906 MVT::MvTimesMatAddMv( one, *Utmp, Ftmp, zero, *U1tmp );
1911 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1912 Ctmp = MVT::CloneViewNonConst( *C_, index );
1913 Utmp = MVT::CloneView( *U_, index );
1916 SDM Ctr(keff,blockSize_);
1917 problem_->computeCurrPrecResVec( &*R_ );
1918 MVT::MvTransMv( one, *Ctmp, *R_, Ctr );
1921 RCP<MV> update = MVT::Clone( *problem_->getCurrLHSVec(), blockSize_ );
1922 MVT::MvInit( *update, 0.0 );
1923 MVT::MvTimesMatAddMv( one, *Utmp, Ctr, one, *update );
1924 problem_->updateSolution( update,
true );
1927 MVT::MvTimesMatAddMv( -one, *Ctmp, Ctr, one, *R_ );
1935 if (V_ == Teuchos::null) {
1937 Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
1938 V_ = MVT::Clone( *rhsMV, (numBlocks_+1)*blockSize_ );
1942 if (MVT::GetNumberVecs(*V_) < (numBlocks_+1)*blockSize_ ) {
1943 Teuchos::RCP<const MV> tmp = V_;
1944 V_ = MVT::Clone( *tmp, (numBlocks_+1)*blockSize_ );
1949 printer_->stream(
Debug) <<
" No recycled subspace available for RHS index " << std::endl << std::endl;
1951 Teuchos::ParameterList primeList;
1954 primeList.set(
"Num Blocks",numBlocks_-1);
1955 primeList.set(
"Block Size",blockSize_);
1956 primeList.set(
"Recycled Blocks",0);
1957 primeList.set(
"Keep Hessenberg",
true);
1958 primeList.set(
"Initialize Hessenberg",
true);
1960 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1961 if (blockSize_*static_cast<ptrdiff_t>(numBlocks_) > dim) {
1962 ptrdiff_t tmpNumBlocks = 0;
1963 if (blockSize_ == 1)
1964 tmpNumBlocks = dim / blockSize_;
1966 tmpNumBlocks = ( dim - blockSize_) / blockSize_;
1967 printer_->stream(
Warnings) <<
"Belos::BlockGmresSolMgr::solve(): Warning! Requested Krylov subspace dimension is larger than operator dimension!" 1968 << std::endl <<
"The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << tmpNumBlocks << std::endl;
1969 primeList.set(
"Num Blocks",Teuchos::as<int>(tmpNumBlocks));
1973 primeList.set(
"Num Blocks",numBlocks_-1);
1976 Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > block_gmres_iter;
1980 block_gmres_iter->setSize( blockSize_, numBlocks_-1 );
1983 Teuchos::RCP<MV> V_0;
1984 if (currIdx[blockSize_-1] == -1) {
1985 V_0 = MVT::Clone( *(problem_->getInitPrecResVec()), blockSize_ );
1986 problem_->computeCurrPrecResVec( &*V_0 );
1989 V_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
1993 Teuchos::RCP<SDM > z_0 = Teuchos::rcp(
new SDM( blockSize_, blockSize_ ) );
1996 int rank = ortho_->normalize( *V_0, z_0 );
2005 block_gmres_iter->initializeGmres(newstate);
2007 bool primeConverged =
false;
2010 printer_->stream(
Debug) <<
" Preparing to Iterate!!!!" << std::endl << std::endl;
2011 block_gmres_iter->iterate();
2016 if ( convTest_->getStatus() ==
Passed ) {
2017 printer_->stream(
Debug) <<
"We converged during the prime the pump stage" << std::endl << std::endl;
2018 primeConverged = !(expConvTest_->getLOADetected());
2019 if ( expConvTest_->getLOADetected() ) {
2021 loaDetected_ =
true;
2022 printer_->stream(
Warnings) <<
"Belos::BlockGmresSolMgr::solve(): Warning! Solver has experienced a loss of accuracy!" << std::endl;
2028 else if( maxIterTest_->getStatus() ==
Passed ) {
2030 primeConverged =
false;
2036 printer_->stream(
Debug) <<
" We did not converge on priming cycle of Block GMRES. Therefore we recycle and restart. " << std::endl << std::endl;
2041 if (blockSize_ != 1) {
2042 printer_->stream(
Errors) <<
"Error! Caught std::exception in BlockGmresIter::iterate() at iteration " 2043 << block_gmres_iter->getNumIters() << std::endl
2044 << e.what() << std::endl;
2045 if (convTest_->getStatus() !=
Passed)
2046 primeConverged =
false;
2050 block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
2052 sTest_->checkStatus( &*block_gmres_iter );
2053 if (convTest_->getStatus() !=
Passed)
2054 isConverged =
false;
2057 catch (
const std::exception &e) {
2058 printer_->stream(
Errors) <<
"Error! Caught std::exception in BlockGmresIter::iterate() at iteration " 2059 << block_gmres_iter->getNumIters() << std::endl
2060 << e.what() << std::endl;
2068 RCP<MV> update = block_gmres_iter->getCurrentUpdate();
2069 problem_->updateSolution( update,
true );
2073 problem_->computeCurrPrecResVec( &*R_ );
2076 newstate = block_gmres_iter->getState();
2079 H_->assign(*(newstate.
H));
2088 V_ = rcp_const_cast<MV>(newstate.
V);
2089 newstate.
V.release();
2091 buildRecycleSpaceKryl(keff, block_gmres_iter);
2092 printer_->stream(
Debug) <<
"Generated recycled subspace using RHS index " << currIdx[0] <<
" of dimension " << keff << std::endl << std::endl;
2096 if (primeConverged) {
2117 problem_->setCurrLS();
2120 startPtr += numCurrRHS;
2121 numRHS2Solve -= numCurrRHS;
2122 if ( numRHS2Solve > 0 ) {
2123 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
2124 if ( adaptiveBlockSize_ ) {
2125 blockSize_ = numCurrRHS;
2126 currIdx.resize( numCurrRHS );
2127 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2130 currIdx.resize( blockSize_ );
2131 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2132 for (
int i=numCurrRHS; i<blockSize_; ++i) currIdx[i] = -1;
2135 problem_->setLSIndex( currIdx );
2138 currIdx.resize( numRHS2Solve );
2148 Teuchos::ParameterList blockgcrodrList;
2149 blockgcrodrList.set(
"Num Blocks",numBlocks_);
2150 blockgcrodrList.set(
"Block Size",blockSize_);
2151 blockgcrodrList.set(
"Recycled Blocks",keff);
2153 Teuchos::RCP<BlockGCRODRIter<ScalarType,MV,OP> > block_gcrodr_iter;
2157 index.resize( blockSize_ );
2158 for(
int ii = 0; ii < blockSize_; ii++) index[ii] = ii;
2159 Teuchos::RCP<MV> V0 = MVT::CloneViewNonConst( *V_, index );
2162 MVT::Assign(*R_,*V0);
2165 for(
int i=0; i < keff; i++){ index[i] = i;};
2166 B_ = rcp(
new SDM(Teuchos::View, *G_, keff, numBlocks_*blockSize_, 0, keff));
2167 H_ = rcp(
new SDM(Teuchos::View, *G_, (numBlocks_-1)*blockSize_ + blockSize_, (numBlocks_-1)*blockSize_, keff ,keff ));
2171 newstate.
U = MVT::CloneViewNonConst(*U_, index);
2172 newstate.
C = MVT::CloneViewNonConst(*C_, index);
2174 newstate.
curDim = blockSize_;
2175 block_gcrodr_iter -> initialize(newstate);
2177 int numRestarts = 0;
2181 block_gcrodr_iter -> iterate();
2186 if( convTest_->getStatus() ==
Passed ) {
2194 else if(maxIterTest_->getStatus() ==
Passed ){
2196 isConverged =
false;
2203 else if (block_gcrodr_iter->getCurSubspaceDim() == block_gcrodr_iter->getMaxSubspaceDim()){
2208 Teuchos::RCP<MV> update = block_gcrodr_iter->getCurrentUpdate();
2209 problem_->updateSolution(update,
true);
2210 buildRecycleSpaceAugKryl(block_gcrodr_iter);
2212 printer_->stream(
Debug) <<
" Generated new recycled subspace using RHS index " << currIdx[0] <<
" of dimension " << keff << std::endl << std::endl;
2214 if(numRestarts >= maxRestarts_) {
2215 isConverged =
false;
2221 printer_ -> stream(
Debug) <<
" Performing restart number " << numRestarts <<
" of " << maxRestarts_ << std::endl << std::endl;
2224 problem_->computeCurrPrecResVec( &*R_ );
2225 index.resize( blockSize_ );
2226 for (
int ii=0; ii<blockSize_; ++ii) index[ii] = ii;
2227 Teuchos::RCP<MV> V0 = MVT::CloneViewNonConst( *V_, index );
2228 MVT::SetBlock(*R_,index,*V0);
2232 index.resize( numBlocks_*blockSize_ );
2233 for (
int ii=0; ii<(numBlocks_*blockSize_); ++ii) index[ii] = ii;
2234 restartState.
V = MVT::CloneViewNonConst( *V_, index );
2235 index.resize( keff );
2236 for (
int ii=0; ii<keff; ++ii) index[ii] = ii;
2237 restartState.
U = MVT::CloneViewNonConst( *U_, index );
2238 restartState.
C = MVT::CloneViewNonConst( *C_, index );
2239 B_ = rcp(
new SDM(Teuchos::View, *G_, keff, (numBlocks_-1)*blockSize_, 0, keff));
2240 H_ = rcp(
new SDM(Teuchos::View, *G_, numBlocks_*blockSize_, (numBlocks_-1)*blockSize_, keff ,keff ));
2241 restartState.
B = B_;
2242 restartState.
H = H_;
2243 restartState.
curDim = blockSize_;
2244 block_gcrodr_iter->initialize(restartState);
2253 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Belos::BlockGCRODRSolMgr::solve(): Invalid return from BlockGCRODRIter::iterate().");
2259 block_gcrodr_iter->updateLSQR( block_gcrodr_iter->getCurSubspaceDim() );
2261 sTest_->checkStatus( &*block_gcrodr_iter );
2262 if (convTest_->getStatus() !=
Passed) isConverged =
false;
2265 catch(
const std::exception &e){
2266 printer_->stream(
Errors) <<
"Error! Caught exception in BlockGCRODRIter::iterate() at iteration " 2267 << block_gcrodr_iter->getNumIters() << std::endl
2268 << e.what() << std::endl;
2275 Teuchos::RCP<MV> update = block_gcrodr_iter->getCurrentUpdate();
2276 problem_->updateSolution( update,
true );
2295 problem_->setCurrLS();
2298 startPtr += numCurrRHS;
2299 numRHS2Solve -= numCurrRHS;
2300 if ( numRHS2Solve > 0 ) {
2301 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
2302 if ( adaptiveBlockSize_ ) {
2303 blockSize_ = numCurrRHS;
2304 currIdx.resize( numCurrRHS );
2305 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2308 currIdx.resize( blockSize_ );
2309 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2310 for (
int i=numCurrRHS; i<blockSize_; ++i) currIdx[i] = -1;
2313 problem_->setLSIndex( currIdx );
2316 currIdx.resize( numRHS2Solve );
2320 if (!builtRecycleSpace_) {
2321 buildRecycleSpaceAugKryl(block_gcrodr_iter);
2322 printer_->stream(
Debug) <<
" Generated new recycled subspace using RHS index " << currIdx[0] <<
" of dimension " << keff << std::endl << std::endl;
2330 #ifdef BELOS_TEUCHOS_TIME_MONITOR 2336 numIters_ = maxIterTest_->getNumIters();
2339 const std::vector<MagnitudeType>* pTestValues = NULL;
2340 pTestValues = impConvTest_->getTestValue();
2341 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
2342 "Belos::BlockGCRODRSolMgr::solve(): The implicit convergence test's " 2343 "getTestValue() method returned NULL. Please report this bug to the " 2344 "Belos developers.");
2345 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
2346 "Belos::BlockGCRODRSolMgr::solve(): The implicit convergence test's " 2347 "getTestValue() method returned a vector of length zero. Please report " 2348 "this bug to the Belos developers.");
2352 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Collection of types and exceptions used within the Belos solvers.
OperatorTraits< ScalarType, MV, OP > OPT
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
int getHarmonicVecsAugKryl(int keff, int m, const SDM &HH, const Teuchos::RCP< const MV > &VV, SDM &PP)
Teuchos::LAPACK< int, ScalarType > lapack
Class which manages the output and verbosity of the Belos solvers.
Structure to contain pointers to BlockGCRODRIter state variables.
static const bool adaptiveBlockSize_default_
virtual ~BlockGCRODRSolMgr()
Destructor.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > ortho_
Orthogonalization manager.
bool loaDetected_
Whether a loss of accuracy was detected during the solve.
OrthoManagerFactory< ScalarType, MV, OP > ortho_factory_type
ScaleType
The type of scaling to use on the residual norm value.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B
The projection of the Krylov subspace against the recycled subspace *
Exception thrown to signal error in a status test during Belos::StatusTest::checkStatus().
void buildRecycleSpaceAugKryl(Teuchos::RCP< BlockGCRODRIter< ScalarType, MV, OP > > gcrodr_iter)
This class implements the block GMRES iteration, where a block Krylov subspace is constructed...
Teuchos::RCP< const MV > V
The current Krylov basis.
BlockGCRODRIterOrthoFailure is thrown when the BlockGCRODRIter object is unable to compute independen...
Teuchos::ScalarTraits< MagnitudeType > SMT
Thrown when the linear problem was not set up correctly.
A factory class for generating StatusTestOutput objects.
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Set the parameters the solver should use to solve the linear problem.
bool isSet_
Whether setParameters() successfully finished setting parameters.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > z
The current right-hand side of the least squares system RY = Z.
An implementation of StatusTestResNorm using a family of residual norms.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
Teuchos::SerialDenseVector< int, ScalarType > SDV
Structure to contain pointers to GmresIteration state variables.
int getNumIters() const
Get the iteration count for the most recent call to solve().
Belos::StatusTest class for specifying a maximum number of iterations.
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > expConvTest_
Teuchos::RCP< Teuchos::ParameterList > params_
This solver's current parameter list.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Get a parameter list containing the valid parameters for this object.
BlockGCRODRSolMgrRecyclingFailure(const std::string &what_arg)
void reset(const ResetType type)
Performs a reset of the solver manager specified by the ResetType.
Thrown if any problem occurs in using or creating the recycle subspace.
A factory class for generating StatusTestOutput objects.
Teuchos::RCP< MV > V
The current Krylov basis.
BlockGCRODRSolMgrOrthoFailure(const std::string &what_arg)
int curDim
The current dimension of the reduction.
Traits class which defines basic operations on multivectors.
std::string description() const
A description of the Block GCRODR solver manager.
Belos::StatusTest for logically combining several status tests.
MultiVecTraits< ScalarType, MV > MVT
Teuchos::RCP< Teuchos::Time > timerSolve_
Timer for solve().
Teuchos::ScalarTraits< MagnitudeType > MT
void sort(std::vector< MagnitudeType > &dlist, int n, std::vector< int > &iperm)
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
BlockGCRODRSolMgrLAPACKFailure(const std::string &what_arg)
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
The current linear problem to solve.
A solver manager for the Block GCRO-DR (Block Recycling GMRES) linear solver.
Pure virtual base class which describes the basic interface for a solver manager. ...
int curDim
The current dimension of the reduction.
static const std::string recycleMethod_default_
void initializeStateStorage()
std::string recycleMethod_
MagnitudeType achievedTol_
Belos concrete class for performing the block GCRO-DR (block GMRES with recycling) iteration...
std::vector< ScalarType > work_
bool isLOADetected() const
Whether a loss of accuracy was detected during the most recent solve.
BlockGCRODRSolMgr()
Default constructor.
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem)
Set the linear problem to solve on the next call to solve().
MagnitudeType achievedTol() const
Get the residual for the most recent call to solve().
Teuchos::RCP< OutputManager< ScalarType > > printer_
ReturnType
Whether the Belos solve converged for all linear systems.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get current linear problem being solved for in this object.
The Belos::SolverManager is a templated virtual base class that defines the basic interface that any ...
Thrown when the solution manager was unable to orthogonalize the basis vectors.
Teuchos::RCP< MV > U
The recycled subspace and its projection.
Teuchos::ScalarTraits< ScalarType > SCT
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > impConvTest_
Belos concrete class for performing the block GMRES iteration.
Implementation of the Block GCRO-DR (Block Recycling GMRES) iteration.
Teuchos::RCP< std::ostream > outputStream_
void buildRecycleSpaceKryl(int &keff, Teuchos::RCP< BlockGmresIter< ScalarType, MV, OP > > block_gmres_iter)
int setTolerance(MagnitudeType tolerance)
Set the value of the tolerance.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > create(const Teuchos::RCP< OutputManager< ScalarType > > &printer, Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test, int mod, int printStates)
Create the StatusTestOutput object specified by the outputStyle.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > convTest_
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
MagnitudeType orthoKappa_
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
A class for extending the status testing capabilities of Belos via logical combinations.
int getHarmonicVecsKryl(int m, const SDM &HH, SDM &PP)
Class which defines basic traits for the operator type.
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
Parent class to all Belos exceptions.
Pure virtual base class which augments the basic interface for a Gmres linear solver iteration...
GmresIterationOrthoFailure is thrown when the GmresIteration object is unable to compute independent ...
Belos header file which uses auto-configuration information to include necessary C++ headers...
BlockGCRODRSolMgrLinearProblemFailure(const std::string &what_arg)
std::vector< ScalarType > tau_
Teuchos::RCP< const Teuchos::ParameterList > defaultParams_
Default parameter list.
Teuchos::SerialDenseMatrix< int, ScalarType > SDM
ReturnType solve()
Solve the current linear problem.
bool builtRecycleSpace_
Whether we have generated or regenerated a recycle space yet this solve.
Thrown when an LAPACK call fails.
ortho_factory_type orthoFactory_
Factory for creating MatOrthoManager subclass instances.
Belos concrete class for performing the block, flexible GMRES iteration.
OutputType
Style of output used to display status test information.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.