45 #ifndef BELOS_BLOCK_GCRODR_SOLMGR_HPP 46 #define BELOS_BLOCK_GCRODR_SOLMGR_HPP 62 #include "Teuchos_BLAS.hpp" 63 #include "Teuchos_LAPACK.hpp" 64 #include "Teuchos_as.hpp" 65 #ifdef BELOS_TEUCHOS_TIME_MONITOR 66 #include "Teuchos_TimeMonitor.hpp" 67 #endif // BELOS_TEUCHOS_TIME_MONITOR 126 template<
class ScalarType,
class MV,
class OP>
132 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
133 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
MagnitudeType;
134 typedef Teuchos::ScalarTraits<MagnitudeType>
MT;
135 typedef Teuchos::ScalarTraits<MagnitudeType>
SMT;
137 typedef Teuchos::SerialDenseMatrix<int,ScalarType>
SDM;
138 typedef Teuchos::SerialDenseVector<int,ScalarType>
SDV;
177 const Teuchos::RCP<Teuchos::ParameterList> &pl);
224 TEUCHOS_TEST_FOR_EXCEPTION(problem.is_null(), std::invalid_argument,
225 "Belos::BlockGCRODRSolMgr::setProblem: The input LinearProblem cannot be null.");
228 if (! problem->isProblemSet()) {
229 const bool success = problem->setProblem();
230 TEUCHOS_TEST_FOR_EXCEPTION(success, std::runtime_error,
231 "Belos::BlockGCRODRSolMgr::setProblem: Calling the input LinearProblem's setProblem() method failed. This likely means that the " 232 "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.");
239 void setParameters(
const Teuchos::RCP<Teuchos::ParameterList> ¶ms );
318 const Teuchos::RCP<const MV>& VV,
322 void sort (std::vector<MagnitudeType>& dlist,
int n, std::vector<int>& iperm);
328 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >
problem_;
335 Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
sTest_;
349 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> >
ortho_;
401 Teuchos::RCP<SDM >
G_;
402 Teuchos::RCP<SDM >
H_;
403 Teuchos::RCP<SDM >
B_;
408 Teuchos::RCP<SDM >
F_;
427 template<
class ScalarType,
class MV,
class OP>
430 template<
class ScalarType,
class MV,
class OP>
437 template<
class ScalarType,
class MV,
class OP>
443 template<
class ScalarType,
class MV,
class OP>
446 const Teuchos::RCP<Teuchos::ParameterList> &pl ) {
451 TEUCHOS_TEST_FOR_EXCEPTION(problem.is_null(), std::invalid_argument,
452 "Belos::BlockGCRODR constructor: The solver manager's constructor needs " 453 "the linear problem argument 'problem' to be nonnull.");
463 template<
class ScalarType,
class MV,
class OP>
465 adaptiveBlockSize_ = adaptiveBlockSize_default_;
466 recycleMethod_ = recycleMethod_default_;
468 loaDetected_ =
false;
469 builtRecycleSpace_ =
false;
540 template<
class ScalarType,
class MV,
class OP>
542 std::ostringstream oss;
543 oss <<
"Belos::BlockGCRODRSolMgr<" << SCT::name() <<
", ...>";
545 oss <<
"Ortho Type='"<<orthoType_ ;
546 oss <<
", Num Blocks=" <<numBlocks_;
547 oss <<
", Block Size=" <<blockSize_;
548 oss <<
", Num Recycle Blocks=" << recycledBlocks_;
549 oss <<
", Max Restarts=" << maxRestarts_;
554 template<
class ScalarType,
class MV,
class OP>
555 Teuchos::RCP<const Teuchos::ParameterList>
557 using Teuchos::ParameterList;
558 using Teuchos::parameterList;
560 using Teuchos::rcpFromRef;
562 if (defaultParams_.is_null()) {
563 RCP<ParameterList> pl = parameterList ();
565 const MagnitudeType convTol = SMT::squareroot (SCT::magnitude (SCT::eps()));
566 const int maxRestarts = 1000;
567 const int maxIters = 5000;
568 const int blockSize = 2;
569 const int numBlocks = 100;
570 const int numRecycledBlocks = 25;
573 const int outputFreq = 1;
574 RCP<std::ostream> outputStream = rcpFromRef (std::cout);
575 const std::string impResScale (
"Norm of Preconditioned Initial Residual");
576 const std::string expResScale (
"Norm of Initial Residual");
577 const std::string timerLabel (
"Belos");
578 const std::string orthoType (
"ICGS");
579 RCP<const ParameterList> orthoParams = orthoFactory_.getDefaultParameters (orthoType);
589 pl->set (
"Convergence Tolerance", convTol,
590 "The tolerance that the solver needs to achieve in order for " 591 "the linear system(s) to be declared converged. The meaning " 592 "of this tolerance depends on the convergence test details.");
593 pl->set(
"Maximum Restarts", maxRestarts,
594 "The maximum number of restart cycles (not counting the first) " 595 "allowed for each set of right-hand sides solved.");
596 pl->set(
"Maximum Iterations", maxIters,
597 "The maximum number of iterations allowed for each set of " 598 "right-hand sides solved.");
599 pl->set(
"Block Size", blockSize,
600 "How many linear systems to solve at once.");
601 pl->set(
"Num Blocks", numBlocks,
602 "The maximum number of blocks allowed in the Krylov subspace " 603 "for each set of right-hand sides solved. This includes the " 604 "initial block. Total Krylov basis storage, not counting the " 605 "recycled subspace, is \"Num Blocks\" times \"Block Size\".");
606 pl->set(
"Num Recycled Blocks", numRecycledBlocks,
607 "The maximum number of vectors in the recycled subspace." );
608 pl->set(
"Verbosity", verbosity,
609 "What type(s) of solver information should be written " 610 "to the output stream.");
611 pl->set(
"Output Style", outputStyle,
612 "What style is used for the solver information to write " 613 "to the output stream.");
614 pl->set(
"Output Frequency", outputFreq,
615 "How often convergence information should be written " 616 "to the output stream.");
617 pl->set(
"Output Stream", outputStream,
618 "A reference-counted pointer to the output stream where all " 619 "solver output is sent.");
620 pl->set(
"Implicit Residual Scaling", impResScale,
621 "The type of scaling used in the implicit residual convergence test.");
622 pl->set(
"Explicit Residual Scaling", expResScale,
623 "The type of scaling used in the explicit residual convergence test.");
624 pl->set(
"Timer Label", timerLabel,
625 "The string to use as a prefix for the timer labels.");
626 pl->set(
"Orthogonalization", orthoType,
627 "The orthogonalization method to use. Valid options: " +
628 orthoFactory_.validNamesString());
629 pl->set (
"Orthogonalization Parameters", *orthoParams,
630 "Sublist of parameters specific to the orthogonalization method to use.");
631 pl->set(
"Orthogonalization Constant", orthoKappa,
632 "When using DGKS orthogonalization: the \"depTol\" constant, used " 633 "to determine whether another step of classical Gram-Schmidt is " 634 "necessary. Otherwise ignored. Nonpositive values are ignored.");
637 return defaultParams_;
640 template<
class ScalarType,
class MV,
class OP>
644 using Teuchos::isParameterType;
645 using Teuchos::getParameter;
647 using Teuchos::ParameterList;
648 using Teuchos::parameterList;
651 using Teuchos::rcp_dynamic_cast;
652 using Teuchos::rcpFromRef;
653 using Teuchos::Exceptions::InvalidParameter;
654 using Teuchos::Exceptions::InvalidParameterName;
655 using Teuchos::Exceptions::InvalidParameterType;
658 RCP<const ParameterList> defaultParams = getValidParameters();
660 if (params.is_null()) {
665 params_ = parameterList (*defaultParams);
677 params->validateParametersAndSetDefaults (*defaultParams);
693 const int maxRestarts = params_->get<
int> (
"Maximum Restarts");
694 TEUCHOS_TEST_FOR_EXCEPTION(maxRestarts <= 0, std::invalid_argument,
695 "Belos::BlockGCRODRSolMgr: The \"Maximum Restarts\" parameter " 696 "must be nonnegative, but you specified a negative value of " 697 << maxRestarts <<
".");
699 const int maxIters = params_->get<
int> (
"Maximum Iterations");
700 TEUCHOS_TEST_FOR_EXCEPTION(maxIters <= 0, std::invalid_argument,
701 "Belos::BlockGCRODRSolMgr: The \"Maximum Iterations\" parameter " 702 "must be positive, but you specified a nonpositive value of " 705 const int numBlocks = params_->get<
int> (
"Num Blocks");
706 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument,
707 "Belos::BlockGCRODRSolMgr: The \"Num Blocks\" parameter must be " 708 "positive, but you specified a nonpositive value of " << numBlocks
711 const int blockSize = params_->get<
int> (
"Block Size");
712 TEUCHOS_TEST_FOR_EXCEPTION(blockSize <= 0, std::invalid_argument,
713 "Belos::BlockGCRODRSolMgr: The \"Block Size\" parameter must be " 714 "positive, but you specified a nonpositive value of " << blockSize
717 const int recycledBlocks = params_->get<
int> (
"Num Recycled Blocks");
718 TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks <= 0, std::invalid_argument,
719 "Belos::BlockGCRODRSolMgr: The \"Num Recycled Blocks\" parameter must " 720 "be positive, but you specified a nonpositive value of " 721 << recycledBlocks <<
".");
722 TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks >= numBlocks,
723 std::invalid_argument,
"Belos::BlockGCRODRSolMgr: The \"Num Recycled " 724 "Blocks\" parameter must be less than the \"Num Blocks\" parameter, " 725 "but you specified \"Num Recycled Blocks\" = " << recycledBlocks
726 <<
" and \"Num Blocks\" = " << numBlocks <<
".");
730 maxRestarts_ = maxRestarts;
731 maxIters_ = maxIters;
732 numBlocks_ = numBlocks;
733 blockSize_ = blockSize;
734 recycledBlocks_ = recycledBlocks;
741 std::string tempLabel = params_->get<std::string> (
"Timer Label");
742 const bool labelChanged = (tempLabel != label_);
744 #ifdef BELOS_TEUCHOS_TIME_MONITOR 745 std::string solveLabel = label_ +
": BlockGCRODRSolMgr total solve time";
746 if (timerSolve_.is_null()) {
748 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
749 }
else if (labelChanged) {
755 Teuchos::TimeMonitor::clearCounter (solveLabel);
756 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
758 #endif // BELOS_TEUCHOS_TIME_MONITOR 764 if (params_->isParameter (
"Verbosity")) {
765 if (isParameterType<int> (*params_,
"Verbosity")) {
766 verbosity_ = params_->get<
int> (
"Verbosity");
769 verbosity_ = (int) getParameter<MsgType> (*params_,
"Verbosity");
776 if (params_->isParameter (
"Output Style")) {
777 if (isParameterType<int> (*params_,
"Output Style")) {
778 outputStyle_ = params_->get<
int> (
"Output Style");
781 outputStyle_ = (int) getParameter<OutputType> (*params_,
"Output Style");
809 if (params_->isParameter (
"Output Stream")) {
811 outputStream_ = getParameter<RCP<std::ostream> > (*params_,
"Output Stream");
813 catch (InvalidParameter&) {
814 outputStream_ = rcpFromRef (std::cout);
821 if (outputStream_.is_null()) {
822 outputStream_ = rcp (
new Teuchos::oblackholestream);
826 outputFreq_ = params_->get<
int> (
"Output Frequency");
829 if (! outputTest_.is_null()) {
830 outputTest_->setOutputFrequency (outputFreq_);
838 if (printer_.is_null()) {
841 printer_->setVerbosity (verbosity_);
842 printer_->setOStream (outputStream_);
853 if (params_->isParameter (
"Orthogonalization")) {
854 const std::string& tempOrthoType =
855 params_->get<std::string> (
"Orthogonalization");
857 if (! orthoFactory_.isValidName (tempOrthoType)) {
858 std::ostringstream os;
859 os <<
"Belos::BlockGCRODRSolMgr: Invalid orthogonalization name \"" 860 << tempOrthoType <<
"\". The following are valid options " 861 <<
"for the \"Orthogonalization\" name parameter: ";
862 orthoFactory_.printValidNames (os);
863 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument, os.str());
865 if (tempOrthoType != orthoType_) {
867 orthoType_ = tempOrthoType;
884 RCP<ParameterList> orthoParams = sublist (params,
"Orthogonalization Parameters",
true);
885 TEUCHOS_TEST_FOR_EXCEPTION(orthoParams.is_null(), std::logic_error,
886 "Failed to get orthogonalization parameters. " 887 "Please report this bug to the Belos developers.");
913 ortho_ = orthoFactory_.makeMatOrthoManager (orthoType_, null, printer_,
914 label_, orthoParams);
926 if (params_->isParameter (
"Orthogonalization Constant")) {
929 if (orthoKappa > 0) {
930 orthoKappa_ = orthoKappa;
933 if (orthoType_ ==
"DGKS" && ! ortho_.is_null()) {
939 rcp_dynamic_cast<ortho_man_type>(ortho_)->setDepTol (orthoKappa_);
949 convTol_ = params_->get<
MagnitudeType> (
"Convergence Tolerance");
950 if (! impConvTest_.is_null()) {
953 if (! expConvTest_.is_null()) {
954 expConvTest_->setTolerance (convTol_);
958 if (params_->isParameter (
"Implicit Residual Scaling")) {
959 std::string tempImpResScale =
960 getParameter<std::string> (*params_,
"Implicit Residual Scaling");
963 if (impResScale_ != tempImpResScale) {
965 impResScale_ = tempImpResScale;
967 if (! impConvTest_.is_null()) {
980 if (params_->isParameter(
"Explicit Residual Scaling")) {
981 std::string tempExpResScale =
982 getParameter<std::string> (*params_,
"Explicit Residual Scaling");
985 if (expResScale_ != tempExpResScale) {
987 expResScale_ = tempExpResScale;
989 if (! expConvTest_.is_null()) {
1007 if (maxIterTest_.is_null()) {
1010 maxIterTest_->setMaxIters (maxIters_);
1015 if (impConvTest_.is_null()) {
1016 impConvTest_ = rcp (
new StatusTestResNorm_t (convTol_));
1022 if (expConvTest_.is_null()) {
1023 expConvTest_ = rcp (
new StatusTestResNorm_t (convTol_));
1024 expConvTest_->defineResForm (StatusTestResNorm_t::Explicit,
Belos::TwoNorm);
1030 if (convTest_.is_null()) {
1031 convTest_ = rcp (
new StatusTestCombo_t (StatusTestCombo_t::SEQ,
1039 sTest_ = rcp (
new StatusTestCombo_t (StatusTestCombo_t::OR,
1040 maxIterTest_, convTest_));
1044 outputTest_ = stoFactory.
create (printer_, sTest_, outputFreq_,
1048 std::string solverDesc =
"Block GCRODR ";
1049 outputTest_->setSolverDesc (solverDesc);
1056 template<
class ScalarType,
class MV,
class OP>
1061 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1064 Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
1067 int KrylSpaDim = (numBlocks_ - 1) * blockSize_;
1070 int augSpaDim = KrylSpaDim + recycledBlocks_ + 1;
1073 int augSpaImgDim = KrylSpaDim + blockSize_ + recycledBlocks_+1;
1097 if (U_ == Teuchos::null) {
1098 U_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1102 if (MVT::GetNumberVecs(*U_) < recycledBlocks_+1) {
1103 Teuchos::RCP<const MV> tmp = U_;
1104 U_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1109 if (C_ == Teuchos::null) {
1110 C_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1114 if (MVT::GetNumberVecs(*C_) < recycledBlocks_+1) {
1115 Teuchos::RCP<const MV> tmp = C_;
1116 C_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1121 if (U1_ == Teuchos::null) {
1122 U1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1126 if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1127 Teuchos::RCP<const MV> tmp = U1_;
1128 U1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1133 if (C1_ == Teuchos::null) {
1134 C1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1138 if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1139 Teuchos::RCP<const MV> tmp = C1_;
1140 C1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1145 if (R_ == Teuchos::null){
1146 R_ = MVT::Clone( *rhsMV, blockSize_ );
1152 if (G_ == Teuchos::null){
1153 G_ = Teuchos::rcp(
new SDM( augSpaImgDim, augSpaDim ) );
1156 if ( (G_->numRows() != augSpaImgDim) || (G_->numCols() != augSpaDim) )
1158 G_->reshape( augSpaImgDim, augSpaDim );
1160 G_->putScalar(zero);
1164 if (H_ == Teuchos::null){
1165 H_ = Teuchos::rcp (
new SDM ( Teuchos::View, *G_, KrylSpaDim + blockSize_, KrylSpaDim, recycledBlocks_+1 ,recycledBlocks_+1 ) );
1169 if (F_ == Teuchos::null){
1170 F_ = Teuchos::rcp(
new SDM( recycledBlocks_+1, recycledBlocks_+1 ) );
1173 if ( (F_->numRows() != recycledBlocks_+1) || (F_->numCols() != recycledBlocks_+1) ){
1174 F_->reshape( recycledBlocks_+1, recycledBlocks_+1 );
1177 F_->putScalar(zero);
1180 if (PP_ == Teuchos::null){
1181 PP_ = Teuchos::rcp(
new SDM( augSpaImgDim, recycledBlocks_+1 ) );
1184 if ( (PP_->numRows() != augSpaImgDim) || (PP_->numCols() != recycledBlocks_+1) ){
1185 PP_->reshape( augSpaImgDim, recycledBlocks_+1 );
1190 if (HP_ == Teuchos::null)
1191 HP_ = Teuchos::rcp(
new SDM( augSpaImgDim, augSpaDim ) );
1193 if ( (HP_->numRows() != augSpaImgDim) || (HP_->numCols() != augSpaDim) ){
1194 HP_->reshape( augSpaImgDim, augSpaDim );
1199 tau_.resize(recycledBlocks_+1);
1202 work_.resize(recycledBlocks_+1);
1205 ipiv_.resize(recycledBlocks_+1);
1209 template<
class ScalarType,
class MV,
class OP>
1212 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1213 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1215 int p = block_gmres_iter->getState().curDim;
1216 std::vector<int> index(keff);
1220 SDM HH(Teuchos::Copy, *H_, p+blockSize_, p);
1221 if(recycledBlocks_ >= p + blockSize_){
1225 for (
int ii=0; ii < p; ++ii) index[ii] = ii;
1226 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1227 MVT::SetBlock(*V_, index, *Utmp);
1232 Teuchos::RCP<SDM > PPtmp = Teuchos::rcp (
new SDM ( Teuchos::View, *PP_, p, recycledBlocks_+1 ) );
1233 if(recycleMethod_ ==
"harmvecs"){
1234 keff = getHarmonicVecsKryl(p, HH, *PPtmp);
1235 printer_->stream(
Debug) <<
"keff = " << keff << std::endl;
1238 PPtmp = Teuchos::rcp (
new SDM ( Teuchos::View, *PP_, p, keff ) );
1241 for (
int ii=0; ii<keff; ++ii) index[ii] = ii;
1242 Teuchos::RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1243 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1244 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1246 for (
int ii=0; ii < p; ++ii) index[ii] = ii;
1247 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1251 MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
1254 SDM HPtmp( Teuchos::View, *HP_, p+blockSize_, keff );
1255 HPtmp.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, *H_, *PPtmp, zero );
1259 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1260 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
1263 lwork =
static_cast<int> (Teuchos::ScalarTraits<ScalarType>::magnitude (work_[0]));
1264 work_.resize(lwork);
1265 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1266 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
1270 SDM Rtmp( Teuchos::View, *F_, keff, keff );
1271 for(
int ii=0;ii<keff;ii++) {
for(
int jj=ii;jj<keff;jj++) Rtmp(ii,jj) = HPtmp(ii,jj); }
1273 lapack.UNGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1274 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _UNGQR failed to construct the Q factor.");
1278 index.resize( p+blockSize_ );
1279 for (
int ii=0; ii < (p+blockSize_); ++ii) { index[ii] = ii; }
1280 Vtmp = MVT::CloneView( *V_, index );
1281 MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
1288 ipiv_.resize(Rtmp.numRows());
1289 lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
1290 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1292 lwork = Rtmp.numRows();
1293 work_.resize(lwork);
1294 lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1297 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
1303 template<
class ScalarType,
class MV,
class OP>
1305 const MagnitudeType one = Teuchos::ScalarTraits<ScalarType>::one();
1306 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1308 std::vector<MagnitudeType> d(keff);
1309 std::vector<ScalarType> dscalar(keff);
1310 std::vector<int> index(numBlocks_+1);
1319 if(recycledBlocks_ >= keff + p){
1322 for (
int ii=0; ii < p; ++ii) { index[ii] = keff+ii; }
1323 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1324 for (
int ii=0; ii < p; ++ii) { index[ii] =ii; }
1325 MVT::SetBlock(*V_, index, *Utmp);
1332 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1333 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1335 dscalar.resize(keff);
1336 MVT::MvNorm( *Utmp, d );
1337 for (
int i=0; i<keff; ++i) {
1339 dscalar[i] = (ScalarType)d[i];
1341 MVT::MvScale( *Utmp, dscalar );
1346 Teuchos::RCP<SDM> Gtmp = Teuchos::rcp(
new SDM( Teuchos::View, *G_, p+keff, p+keff-blockSize_ ) );
1349 for (
int i=0; i<keff; ++i)
1350 (*Gtmp)(i,i) = d[i];
1357 SDM PPtmp( Teuchos::View, *PP_, p+keff-blockSize_, recycledBlocks_+1 );
1358 keff_new = getHarmonicVecsAugKryl( keff, p-blockSize_, *Gtmp, oldState.
V, PPtmp );
1365 Teuchos::RCP<MV> U1tmp;
1367 index.resize( keff );
1368 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1369 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1370 index.resize( keff_new );
1371 for (
int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
1372 U1tmp = MVT::CloneViewNonConst( *U1_, index );
1373 SDM PPtmp( Teuchos::View, *PP_, keff, keff_new );
1374 MVT::MvTimesMatAddMv( one, *Utmp, PPtmp, zero, *U1tmp );
1379 index.resize(p-blockSize_);
1380 for (
int ii=0; ii < p-blockSize_; ii++) { index[ii] = ii; }
1381 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1382 SDM PPtmp( Teuchos::View, *PP_, p-blockSize_, keff_new, keff );
1383 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
1387 SDM HPtmp( Teuchos::View, *HP_, p+keff, keff_new );
1389 SDM PPtmp( Teuchos::View, *PP_, p-blockSize_+keff, keff_new );
1390 HPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*Gtmp,PPtmp,zero);
1394 int info = 0, lwork = -1;
1395 tau_.resize(keff_new);
1396 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1397 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
1399 lwork =
static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work_[0]));
1400 work_.resize(lwork);
1401 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1402 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
1406 SDM Ftmp( Teuchos::View, *F_, keff_new, keff_new );
1407 for(
int i=0;i<keff_new;i++) {
for(
int j=i;j<keff_new;j++) Ftmp(i,j) = HPtmp(i,j); }
1409 lapack.UNGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1410 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _UNGQR failed to construct the Q factor.");
1417 Teuchos::RCP<MV> C1tmp;
1420 for (
int i=0; i < keff; i++) { index[i] = i; }
1421 Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
1422 index.resize(keff_new);
1423 for (
int i=0; i < keff_new; i++) { index[i] = i; }
1424 C1tmp = MVT::CloneViewNonConst( *C1_, index );
1425 SDM PPtmp( Teuchos::View, *HP_, keff, keff_new );
1426 MVT::MvTimesMatAddMv( one, *Ctmp, PPtmp, zero, *C1tmp );
1431 for (
int i=0; i < p; ++i) { index[i] = i; }
1432 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1433 SDM PPtmp( Teuchos::View, *HP_, p, keff_new, keff, 0 );
1434 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *C1tmp );
1443 ipiv_.resize(Ftmp.numRows());
1444 lapack.GETRF(Ftmp.numRows(),Ftmp.numCols(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&info);
1445 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1448 lwork = Ftmp.numRows();
1449 work_.resize(lwork);
1450 lapack.GETRI(Ftmp.numRows(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1451 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRI failed to compute an LU factorization.");
1454 index.resize(keff_new);
1455 for (
int i=0; i < keff_new; i++) { index[i] = i; }
1456 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1457 MVT::MvTimesMatAddMv( one, *U1tmp, Ftmp, zero, *Utmp );
1462 if (keff != keff_new) {
1464 block_gcrodr_iter->setSize( keff, numBlocks_ );
1466 SDM b1( Teuchos::View, *G_, recycledBlocks_+2, 1, 0, recycledBlocks_ );
1472 template<
class ScalarType,
class MV,
class OP>
1475 int m2 = GG.numCols();
1476 bool xtraVec =
false;
1477 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1478 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1479 std::vector<int> index;
1482 std::vector<MagnitudeType> wr(m2), wi(m2);
1485 std::vector<MagnitudeType> w(m2);
1488 SDM vr(m2,m2,
false);
1491 std::vector<int> iperm(m2);
1494 builtRecycleSpace_ =
true;
1504 SDM A_tmp( keff+m+blockSize_, keff+m );
1509 for (
int i=0; i<keff; ++i) { index[i] = i; }
1510 Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
1511 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1512 SDM A11( Teuchos::View, A_tmp, keff, keff );
1513 MVT::MvTransMv( one, *Ctmp, *Utmp, A11 );
1516 SDM A21( Teuchos::View, A_tmp, m+blockSize_, keff, keff );
1517 index.resize(m+blockSize_);
1518 for (i=0; i < m+blockSize_; i++) { index[i] = i; }
1519 Teuchos::RCP<const MV> Vp = MVT::CloneView( *VV, index );
1520 MVT::MvTransMv( one, *Vp, *Utmp, A21 );
1523 for( i=keff; i<keff+m; i++)
1527 SDM A( m2, A_tmp.numCols() );
1528 A.multiply(
Teuchos::TRANS, Teuchos::NO_TRANS, one, GG, A_tmp, zero );
1536 char balanc=
'P', jobvl=
'N', jobvr=
'V', sense=
'N';
1537 int ld = A.numRows();
1539 int ldvl = ld, ldvr = ld;
1540 int info = 0,ilo = 0,ihi = 0;
1543 std::vector<ScalarType> beta(ld);
1544 std::vector<ScalarType> work(lwork);
1545 std::vector<MagnitudeType> rwork(lwork);
1546 std::vector<MagnitudeType> lscale(ld), rscale(ld);
1547 std::vector<MagnitudeType> rconde(ld), rcondv(ld);
1548 std::vector<int> iwork(ld+6);
1553 lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0],
1554 &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
1555 &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &rwork[0],
1556 &iwork[0], bwork, &info);
1557 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK GGEVX failed to compute eigensolutions.");
1561 for( i=0; i<ld; i++ )
1562 w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] ) / std::abs( beta[i] );
1564 this->sort(w,ld,iperm);
1566 bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
1569 for( i=0; i<recycledBlocks_; i++ )
1570 for( j=0; j<ld; j++ )
1571 PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
1573 if(scalarTypeIsComplex==
false) {
1576 if (wi[iperm[ld-recycledBlocks_]] != 0.0) {
1578 for ( i=ld-recycledBlocks_; i<ld; i++ )
1579 if (wi[iperm[i]] != 0.0) countImag++;
1581 if (countImag % 2) xtraVec =
true;
1585 if (wi[iperm[ld-recycledBlocks_]] > 0.0) {
1586 for( j=0; j<ld; j++ )
1587 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
1590 for( j=0; j<ld; j++ )
1591 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
1599 return recycledBlocks_+1;
1601 return recycledBlocks_;
1605 template<
class ScalarType,
class MV,
class OP>
1607 bool xtraVec =
false;
1608 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1609 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1612 std::vector<MagnitudeType> wr(m), wi(m);
1618 std::vector<MagnitudeType> w(m);
1621 std::vector<int> iperm(m);
1627 builtRecycleSpace_ =
true;
1631 Teuchos::RCP<SDM> harmRitzMatrix = Teuchos::rcp(
new SDM( m, blockSize_));
1634 for(
int i=0; i<=blockSize_-1; i++) (*harmRitzMatrix)[blockSize_-1-i][harmRitzMatrix->numRows()-1-i] = 1;
1637 lapack.GESV(m, blockSize_, HHt.values(), HHt.stride(), &iperm[0], harmRitzMatrix->values(), harmRitzMatrix->stride(), &info);
1639 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK GESV failed to compute a solution.");
1643 Teuchos::SerialDenseMatrix<int, ScalarType> H_lbl(Teuchos::View, HH, blockSize_, blockSize_, (HH).numRows()-blockSize_, (HH).numCols()-blockSize_ );
1644 Teuchos::SerialDenseMatrix<int, ScalarType> H_lbl_t( H_lbl,
Teuchos::TRANS );
1649 Teuchos::RCP<SDM> Htemp = Teuchos::null;
1650 Htemp = Teuchos::rcp(
new SDM(H_lbl_t.numRows(), H_lbl_t.numCols()));
1651 Htemp -> multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, H_lbl_t, H_lbl, zero);
1652 H_lbl_t.assign(*Htemp);
1654 Htemp = Teuchos::rcp(
new SDM(harmRitzMatrix -> numRows(), harmRitzMatrix -> numCols()));
1655 Htemp -> multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, *harmRitzMatrix, H_lbl_t, zero);
1656 harmRitzMatrix -> assign(*Htemp);
1659 int harmColIndex, HHColIndex;
1660 Htemp = Teuchos::rcp(
new SDM(Teuchos::Copy,HH,HH.numRows()-blockSize_,HH.numCols()));
1661 for(
int i = 0; i<blockSize_; i++){
1662 harmColIndex = harmRitzMatrix -> numCols() - i -1;
1664 for(
int j=0; j<m; j++) (*Htemp)[HHColIndex][j] += (*harmRitzMatrix)[harmColIndex][j];
1666 harmRitzMatrix = Htemp;
1674 std::vector<ScalarType> work(1);
1675 std::vector<MagnitudeType> rwork(2*m);
1681 lapack.GEEV(
'N',
'V', m, harmRitzMatrix->values(), harmRitzMatrix->stride(), &wr[0], &wi[0],
1682 vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
1684 lwork = std::abs (static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work[0])));
1685 work.resize( lwork );
1687 lapack.GEEV(
'N',
'V', m, harmRitzMatrix->values(), harmRitzMatrix->stride(), &wr[0], &wi[0],
1688 vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
1690 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK GEEV failed to compute eigensolutions.");
1693 for(
int i=0; i<m; ++i ) w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] );
1695 this->sort(w, m, iperm);
1697 bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
1700 for(
int i=0; i<recycledBlocks_; ++i )
1701 for(
int j=0; j<m; j++ )
1702 PP(j,i) = vr(j,iperm[i]);
1704 if(scalarTypeIsComplex==
false) {
1707 if (wi[iperm[recycledBlocks_-1]] != 0.0) {
1709 for (
int i=0; i<recycledBlocks_; ++i )
1710 if (wi[iperm[i]] != 0.0) countImag++;
1712 if (countImag % 2) xtraVec =
true;
1716 if (wi[iperm[recycledBlocks_-1]] > 0.0) {
1717 for(
int j=0; j<m; ++j )
1718 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
1721 for(
int j=0; j<m; ++j )
1722 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
1730 printer_->stream(
Debug) <<
"Recycled " << recycledBlocks_+1 <<
" vectors" << std::endl;
1731 return recycledBlocks_+1;
1734 printer_->stream(
Debug) <<
"Recycled " << recycledBlocks_ <<
" vectors" << std::endl;
1735 return recycledBlocks_;
1740 template<
class ScalarType,
class MV,
class OP>
1742 int l, r, j, i, flag;
1769 if (dlist[j] > dlist[j - 1]) j = j + 1;
1770 if (dlist[j - 1] > dK) {
1771 dlist[i - 1] = dlist[j - 1];
1772 iperm[i - 1] = iperm[j - 1];
1785 dlist[r] = dlist[0];
1786 iperm[r] = iperm[0];
1800 template<
class ScalarType,
class MV,
class OP>
1804 using Teuchos::rcp_const_cast;
1808 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1809 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1810 std::vector<int> index(numBlocks_+1);
1826 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1827 int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1831 std::vector<int> currIdx;
1833 if ( adaptiveBlockSize_ ) {
1834 blockSize_ = numCurrRHS;
1835 currIdx.resize( numCurrRHS );
1836 for (
int i=0; i<numCurrRHS; ++i)
1837 currIdx[i] = startPtr+i;
1840 currIdx.resize( blockSize_ );
1841 for (
int i=0; i<numCurrRHS; ++i)
1842 currIdx[i] = startPtr+i;
1843 for (
int i=numCurrRHS; i<blockSize_; ++i)
1848 problem_->setLSIndex( currIdx );
1854 loaDetected_ =
false;
1857 bool isConverged =
true;
1860 initializeStateStorage();
1863 Teuchos::ParameterList plist;
1865 while (numRHS2Solve > 0){
1875 printer_->stream(
Debug) <<
" Now solving RHS index " << currIdx[0] <<
" using recycled subspace of dimension " << keff << std::endl << std::endl;
1879 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1880 RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1881 RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1882 problem_->apply( *Utmp, *Ctmp );
1884 RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1888 SDM Ftmp( Teuchos::View, *F_, keff, keff );
1889 int rank = ortho_->normalize(*Ctmp, rcp(&Ftmp,
false));
1891 TEUCHOS_TEST_FOR_EXCEPTION(rank != keff,
BlockGCRODRSolMgrOrthoFailure,
"Belos::BlockGCRODRSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
1896 ipiv_.resize(Ftmp.numRows());
1897 lapack.GETRF(Ftmp.numRows(),Ftmp.numCols(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&info);
1898 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1901 int lwork = Ftmp.numRows();
1902 work_.resize(lwork);
1903 lapack.GETRI(Ftmp.numRows(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1904 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
1907 MVT::MvTimesMatAddMv( one, *Utmp, Ftmp, zero, *U1tmp );
1912 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1913 Ctmp = MVT::CloneViewNonConst( *C_, index );
1914 Utmp = MVT::CloneView( *U_, index );
1917 SDM Ctr(keff,blockSize_);
1918 problem_->computeCurrPrecResVec( &*R_ );
1919 MVT::MvTransMv( one, *Ctmp, *R_, Ctr );
1922 RCP<MV> update = MVT::Clone( *problem_->getCurrLHSVec(), blockSize_ );
1923 MVT::MvInit( *update, 0.0 );
1924 MVT::MvTimesMatAddMv( one, *Utmp, Ctr, one, *update );
1925 problem_->updateSolution( update,
true );
1928 MVT::MvTimesMatAddMv( -one, *Ctmp, Ctr, one, *R_ );
1936 if (V_ == Teuchos::null) {
1938 Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
1939 V_ = MVT::Clone( *rhsMV, (numBlocks_+1)*blockSize_ );
1943 if (MVT::GetNumberVecs(*V_) < (numBlocks_+1)*blockSize_ ) {
1944 Teuchos::RCP<const MV> tmp = V_;
1945 V_ = MVT::Clone( *tmp, (numBlocks_+1)*blockSize_ );
1950 printer_->stream(
Debug) <<
" No recycled subspace available for RHS index " << std::endl << std::endl;
1952 Teuchos::ParameterList primeList;
1955 primeList.set(
"Num Blocks",numBlocks_-1);
1956 primeList.set(
"Block Size",blockSize_);
1957 primeList.set(
"Recycled Blocks",0);
1958 primeList.set(
"Keep Hessenberg",
true);
1959 primeList.set(
"Initialize Hessenberg",
true);
1961 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1962 if (blockSize_*static_cast<ptrdiff_t>(numBlocks_) > dim) {
1963 ptrdiff_t tmpNumBlocks = 0;
1964 if (blockSize_ == 1)
1965 tmpNumBlocks = dim / blockSize_;
1967 tmpNumBlocks = ( dim - blockSize_) / blockSize_;
1968 printer_->stream(
Warnings) <<
"Belos::BlockGmresSolMgr::solve(): Warning! Requested Krylov subspace dimension is larger than operator dimension!" 1969 << std::endl <<
"The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << tmpNumBlocks << std::endl;
1970 primeList.set(
"Num Blocks",Teuchos::as<int>(tmpNumBlocks));
1974 primeList.set(
"Num Blocks",numBlocks_-1);
1977 Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > block_gmres_iter;
1981 block_gmres_iter->setSize( blockSize_, numBlocks_-1 );
1984 Teuchos::RCP<MV> V_0;
1985 if (currIdx[blockSize_-1] == -1) {
1986 V_0 = MVT::Clone( *(problem_->getInitPrecResVec()), blockSize_ );
1987 problem_->computeCurrPrecResVec( &*V_0 );
1990 V_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
1994 Teuchos::RCP<SDM > z_0 = Teuchos::rcp(
new SDM( blockSize_, blockSize_ ) );
1997 int rank = ortho_->normalize( *V_0, z_0 );
2006 block_gmres_iter->initializeGmres(newstate);
2008 bool primeConverged =
false;
2011 printer_->stream(
Debug) <<
" Preparing to Iterate!!!!" << std::endl << std::endl;
2012 block_gmres_iter->iterate();
2017 if ( convTest_->getStatus() ==
Passed ) {
2018 printer_->stream(
Debug) <<
"We converged during the prime the pump stage" << std::endl << std::endl;
2019 primeConverged = !(expConvTest_->getLOADetected());
2020 if ( expConvTest_->getLOADetected() ) {
2022 loaDetected_ =
true;
2023 printer_->stream(
Warnings) <<
"Belos::BlockGmresSolMgr::solve(): Warning! Solver has experienced a loss of accuracy!" << std::endl;
2029 else if( maxIterTest_->getStatus() ==
Passed ) {
2031 primeConverged =
false;
2037 printer_->stream(
Debug) <<
" We did not converge on priming cycle of Block GMRES. Therefore we recycle and restart. " << std::endl << std::endl;
2042 if (blockSize_ != 1) {
2043 printer_->stream(
Errors) <<
"Error! Caught std::exception in BlockGmresIter::iterate() at iteration " 2044 << block_gmres_iter->getNumIters() << std::endl
2045 << e.what() << std::endl;
2046 if (convTest_->getStatus() !=
Passed)
2047 primeConverged =
false;
2051 block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
2053 sTest_->checkStatus( &*block_gmres_iter );
2054 if (convTest_->getStatus() !=
Passed)
2055 isConverged =
false;
2058 catch (
const std::exception &e) {
2059 printer_->stream(
Errors) <<
"Error! Caught std::exception in BlockGmresIter::iterate() at iteration " 2060 << block_gmres_iter->getNumIters() << std::endl
2061 << e.what() << std::endl;
2069 RCP<MV> update = block_gmres_iter->getCurrentUpdate();
2070 problem_->updateSolution( update,
true );
2074 problem_->computeCurrPrecResVec( &*R_ );
2077 newstate = block_gmres_iter->getState();
2080 H_->assign(*(newstate.
H));
2089 V_ = rcp_const_cast<MV>(newstate.
V);
2090 newstate.
V.release();
2092 buildRecycleSpaceKryl(keff, block_gmres_iter);
2093 printer_->stream(
Debug) <<
"Generated recycled subspace using RHS index " << currIdx[0] <<
" of dimension " << keff << std::endl << std::endl;
2097 if (primeConverged) {
2118 problem_->setCurrLS();
2121 startPtr += numCurrRHS;
2122 numRHS2Solve -= numCurrRHS;
2123 if ( numRHS2Solve > 0 ) {
2124 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
2125 if ( adaptiveBlockSize_ ) {
2126 blockSize_ = numCurrRHS;
2127 currIdx.resize( numCurrRHS );
2128 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2131 currIdx.resize( blockSize_ );
2132 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2133 for (
int i=numCurrRHS; i<blockSize_; ++i) currIdx[i] = -1;
2136 problem_->setLSIndex( currIdx );
2139 currIdx.resize( numRHS2Solve );
2149 Teuchos::ParameterList blockgcrodrList;
2150 blockgcrodrList.set(
"Num Blocks",numBlocks_);
2151 blockgcrodrList.set(
"Block Size",blockSize_);
2152 blockgcrodrList.set(
"Recycled Blocks",keff);
2154 Teuchos::RCP<BlockGCRODRIter<ScalarType,MV,OP> > block_gcrodr_iter;
2158 index.resize( blockSize_ );
2159 for(
int ii = 0; ii < blockSize_; ii++) index[ii] = ii;
2160 Teuchos::RCP<MV> V0 = MVT::CloneViewNonConst( *V_, index );
2163 MVT::Assign(*R_,*V0);
2166 for(
int i=0; i < keff; i++){ index[i] = i;};
2167 B_ = rcp(
new SDM(Teuchos::View, *G_, keff, numBlocks_*blockSize_, 0, keff));
2168 H_ = rcp(
new SDM(Teuchos::View, *G_, (numBlocks_-1)*blockSize_ + blockSize_, (numBlocks_-1)*blockSize_, keff ,keff ));
2172 newstate.
U = MVT::CloneViewNonConst(*U_, index);
2173 newstate.
C = MVT::CloneViewNonConst(*C_, index);
2175 newstate.
curDim = blockSize_;
2176 block_gcrodr_iter -> initialize(newstate);
2178 int numRestarts = 0;
2182 block_gcrodr_iter -> iterate();
2187 if( convTest_->getStatus() ==
Passed ) {
2195 else if(maxIterTest_->getStatus() ==
Passed ){
2197 isConverged =
false;
2204 else if (block_gcrodr_iter->getCurSubspaceDim() == block_gcrodr_iter->getMaxSubspaceDim()){
2209 Teuchos::RCP<MV> update = block_gcrodr_iter->getCurrentUpdate();
2210 problem_->updateSolution(update,
true);
2211 buildRecycleSpaceAugKryl(block_gcrodr_iter);
2213 printer_->stream(
Debug) <<
" Generated new recycled subspace using RHS index " << currIdx[0] <<
" of dimension " << keff << std::endl << std::endl;
2215 if(numRestarts >= maxRestarts_) {
2216 isConverged =
false;
2222 printer_ -> stream(
Debug) <<
" Performing restart number " << numRestarts <<
" of " << maxRestarts_ << std::endl << std::endl;
2225 problem_->computeCurrPrecResVec( &*R_ );
2226 index.resize( blockSize_ );
2227 for (
int ii=0; ii<blockSize_; ++ii) index[ii] = ii;
2228 Teuchos::RCP<MV> V0 = MVT::CloneViewNonConst( *V_, index );
2229 MVT::SetBlock(*R_,index,*V0);
2233 index.resize( numBlocks_*blockSize_ );
2234 for (
int ii=0; ii<(numBlocks_*blockSize_); ++ii) index[ii] = ii;
2235 restartState.
V = MVT::CloneViewNonConst( *V_, index );
2236 index.resize( keff );
2237 for (
int ii=0; ii<keff; ++ii) index[ii] = ii;
2238 restartState.
U = MVT::CloneViewNonConst( *U_, index );
2239 restartState.
C = MVT::CloneViewNonConst( *C_, index );
2240 B_ = rcp(
new SDM(Teuchos::View, *G_, keff, (numBlocks_-1)*blockSize_, 0, keff));
2241 H_ = rcp(
new SDM(Teuchos::View, *G_, numBlocks_*blockSize_, (numBlocks_-1)*blockSize_, keff ,keff ));
2242 restartState.
B = B_;
2243 restartState.
H = H_;
2244 restartState.
curDim = blockSize_;
2245 block_gcrodr_iter->initialize(restartState);
2254 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Belos::BlockGCRODRSolMgr::solve(): Invalid return from BlockGCRODRIter::iterate().");
2260 block_gcrodr_iter->updateLSQR( block_gcrodr_iter->getCurSubspaceDim() );
2262 sTest_->checkStatus( &*block_gcrodr_iter );
2263 if (convTest_->getStatus() !=
Passed) isConverged =
false;
2266 catch(
const std::exception &e){
2267 printer_->stream(
Errors) <<
"Error! Caught exception in BlockGCRODRIter::iterate() at iteration " 2268 << block_gcrodr_iter->getNumIters() << std::endl
2269 << e.what() << std::endl;
2276 Teuchos::RCP<MV> update = block_gcrodr_iter->getCurrentUpdate();
2277 problem_->updateSolution( update,
true );
2296 problem_->setCurrLS();
2299 startPtr += numCurrRHS;
2300 numRHS2Solve -= numCurrRHS;
2301 if ( numRHS2Solve > 0 ) {
2302 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
2303 if ( adaptiveBlockSize_ ) {
2304 blockSize_ = numCurrRHS;
2305 currIdx.resize( numCurrRHS );
2306 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2309 currIdx.resize( blockSize_ );
2310 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2311 for (
int i=numCurrRHS; i<blockSize_; ++i) currIdx[i] = -1;
2314 problem_->setLSIndex( currIdx );
2317 currIdx.resize( numRHS2Solve );
2321 if (!builtRecycleSpace_) {
2322 buildRecycleSpaceAugKryl(block_gcrodr_iter);
2323 printer_->stream(
Debug) <<
" Generated new recycled subspace using RHS index " << currIdx[0] <<
" of dimension " << keff << std::endl << std::endl;
2331 #ifdef BELOS_TEUCHOS_TIME_MONITOR 2337 numIters_ = maxIterTest_->getNumIters();
2340 const std::vector<MagnitudeType>* pTestValues = NULL;
2341 pTestValues = impConvTest_->getTestValue();
2342 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
2343 "Belos::BlockGCRODRSolMgr::solve(): The implicit convergence test's " 2344 "getTestValue() method returned NULL. Please report this bug to the " 2345 "Belos developers.");
2346 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
2347 "Belos::BlockGCRODRSolMgr::solve(): The implicit convergence test's " 2348 "getTestValue() method returned a vector of length zero. Please report " 2349 "this bug to the Belos developers.");
2353 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.