45 #ifndef BELOS_BLOCK_GCRODR_SOLMGR_HPP 46 #define BELOS_BLOCK_GCRODR_SOLMGR_HPP 65 #ifdef BELOS_TEUCHOS_TIME_MONITOR 67 #endif // BELOS_TEUCHOS_TIME_MONITOR 126 template<
class ScalarType,
class MV,
class OP>
225 "Belos::BlockGCRODRSolMgr::setProblem: The input LinearProblem cannot be null.");
228 if (! problem->isProblemSet()) {
229 const bool success = problem->setProblem();
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.");
254 problem_->setProblem();
292 void initializeStateStorage();
310 int getHarmonicVecsKryl (
int m,
const SDM& HH, SDM& PP);
315 int getHarmonicVecsAugKryl (
int keff,
322 void sort (std::vector<MagnitudeType>& dlist,
int n, std::vector<int>& iperm);
342 ortho_factory_type orthoFactory_;
367 static const bool adaptiveBlockSize_default_;
368 static const std::string recycleMethod_default_;
371 MagnitudeType convTol_, orthoKappa_, achievedTol_;
372 int blockSize_, maxRestarts_, maxIters_, numIters_;
373 int verbosity_, outputStyle_, outputFreq_;
374 bool adaptiveBlockSize_;
375 std::string orthoType_, recycleMethod_;
376 std::string impResScale_, expResScale_;
384 int numBlocks_, recycledBlocks_;
406 std::vector<ScalarType> tau_;
407 std::vector<ScalarType> work_;
409 std::vector<int> ipiv_;
421 bool builtRecycleSpace_;
427 template<
class ScalarType,
class MV,
class OP>
428 const bool BlockGCRODRSolMgr<ScalarType,MV,OP>::adaptiveBlockSize_default_ =
true;
430 template<
class ScalarType,
class MV,
class OP>
431 const std::string BlockGCRODRSolMgr<ScalarType,MV,OP>::recycleMethod_default_ =
"harmvecs";
437 template<
class ScalarType,
class MV,
class OP>
443 template<
class ScalarType,
class MV,
class OP>
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>
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 (
"DGKS");
579 RCP<const ParameterList> orthoParams = orthoFactory_.getDefaultParameters (orthoType);
586 const MagnitudeType orthoKappa = -SMT::one();
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.");
627 pl->set(
"Orthogonalization", orthoType,
628 "The orthogonalization method to use. Valid options: " +
629 orthoFactory_.validNamesString());
630 pl->set (
"Orthogonalization Parameters", *orthoParams,
631 "Sublist of parameters specific to the orthogonalization method to use.");
632 pl->set(
"Orthogonalization Constant", orthoKappa,
633 "When using DGKS orthogonalization: the \"depTol\" constant, used " 634 "to determine whether another step of classical Gram-Schmidt is " 635 "necessary. Otherwise ignored. Nonpositive values are ignored.");
638 return defaultParams_;
641 template<
class ScalarType,
class MV,
class OP>
645 using Teuchos::isParameterType;
646 using Teuchos::getParameter;
649 using Teuchos::parameterList;
652 using Teuchos::rcp_dynamic_cast;
653 using Teuchos::rcpFromRef;
659 RCP<const ParameterList> defaultParams = getValidParameters();
666 params_ = parameterList (*defaultParams);
694 const int maxRestarts = params_->
get<
int> (
"Maximum Restarts");
696 "Belos::BlockGCRODRSolMgr: The \"Maximum Restarts\" parameter " 697 "must be nonnegative, but you specified a negative value of " 698 << maxRestarts <<
".");
700 const int maxIters = params_->get<
int> (
"Maximum Iterations");
702 "Belos::BlockGCRODRSolMgr: The \"Maximum Iterations\" parameter " 703 "must be positive, but you specified a nonpositive value of " 706 const int numBlocks = params_->get<
int> (
"Num Blocks");
708 "Belos::BlockGCRODRSolMgr: The \"Num Blocks\" parameter must be " 709 "positive, but you specified a nonpositive value of " << numBlocks
712 const int blockSize = params_->get<
int> (
"Block Size");
714 "Belos::BlockGCRODRSolMgr: The \"Block Size\" parameter must be " 715 "positive, but you specified a nonpositive value of " << blockSize
718 const int recycledBlocks = params_->get<
int> (
"Num Recycled Blocks");
720 "Belos::BlockGCRODRSolMgr: The \"Num Recycled Blocks\" parameter must " 721 "be positive, but you specified a nonpositive value of " 722 << recycledBlocks <<
".");
724 std::invalid_argument,
"Belos::BlockGCRODRSolMgr: The \"Num Recycled " 725 "Blocks\" parameter must be less than the \"Num Blocks\" parameter, " 726 "but you specified \"Num Recycled Blocks\" = " << recycledBlocks
727 <<
" and \"Num Blocks\" = " << numBlocks <<
".");
731 maxRestarts_ = maxRestarts;
732 maxIters_ = maxIters;
733 numBlocks_ = numBlocks;
734 blockSize_ = blockSize;
735 recycledBlocks_ = recycledBlocks;
742 std::string tempLabel = params_->get<std::string> (
"Timer Label");
743 const bool labelChanged = (tempLabel != label_);
745 #ifdef BELOS_TEUCHOS_TIME_MONITOR 746 std::string solveLabel = label_ +
": BlockGCRODRSolMgr total solve time";
747 if (timerSolve_.is_null()) {
750 }
else if (labelChanged) {
759 #endif // BELOS_TEUCHOS_TIME_MONITOR 765 if (params_->isParameter (
"Verbosity")) {
766 if (isParameterType<int> (*params_,
"Verbosity")) {
767 verbosity_ = params_->get<
int> (
"Verbosity");
770 verbosity_ = (int) getParameter<MsgType> (*params_,
"Verbosity");
777 if (params_->isParameter (
"Output Style")) {
778 if (isParameterType<int> (*params_,
"Output Style")) {
779 outputStyle_ = params_->get<
int> (
"Output Style");
782 outputStyle_ = (int) getParameter<OutputType> (*params_,
"Output Style");
810 if (params_->isParameter (
"Output Stream")) {
812 outputStream_ = getParameter<RCP<std::ostream> > (*params_,
"Output Stream");
814 catch (InvalidParameter&) {
815 outputStream_ = rcpFromRef (std::cout);
822 if (outputStream_.is_null()) {
827 outputFreq_ = params_->get<
int> (
"Output Frequency");
830 if (! outputTest_.is_null()) {
831 outputTest_->setOutputFrequency (outputFreq_);
839 if (printer_.is_null()) {
842 printer_->setVerbosity (verbosity_);
843 printer_->setOStream (outputStream_);
854 if (params_->isParameter (
"Orthogonalization")) {
855 const std::string& tempOrthoType =
856 params_->get<std::string> (
"Orthogonalization");
858 if (! orthoFactory_.isValidName (tempOrthoType)) {
859 std::ostringstream os;
860 os <<
"Belos::BlockGCRODRSolMgr: Invalid orthogonalization name \"" 861 << tempOrthoType <<
"\". The following are valid options " 862 <<
"for the \"Orthogonalization\" name parameter: ";
863 orthoFactory_.printValidNames (os);
866 if (tempOrthoType != orthoType_) {
868 orthoType_ = tempOrthoType;
885 RCP<ParameterList> orthoParams = sublist (params,
"Orthogonalization Parameters",
true);
887 "Failed to get orthogonalization parameters. " 888 "Please report this bug to the Belos developers.");
914 ortho_ = orthoFactory_.makeMatOrthoManager (orthoType_, null, printer_,
915 label_, orthoParams);
927 if (params_->isParameter (
"Orthogonalization Constant")) {
928 const MagnitudeType orthoKappa =
929 params_->get<MagnitudeType> (
"Orthogonalization Constant");
930 if (orthoKappa > 0) {
931 orthoKappa_ = orthoKappa;
934 if (orthoType_ ==
"DGKS" && ! ortho_.is_null()) {
940 rcp_dynamic_cast<ortho_man_type>(ortho_)->setDepTol (orthoKappa_);
950 convTol_ = params_->get<MagnitudeType> (
"Convergence Tolerance");
951 if (! impConvTest_.is_null()) {
954 if (! expConvTest_.is_null()) {
955 expConvTest_->setTolerance (convTol_);
959 if (params_->isParameter (
"Implicit Residual Scaling")) {
960 std::string tempImpResScale =
961 getParameter<std::string> (*params_,
"Implicit Residual Scaling");
964 if (impResScale_ != tempImpResScale) {
966 impResScale_ = tempImpResScale;
968 if (! impConvTest_.is_null()) {
981 if (params_->isParameter(
"Explicit Residual Scaling")) {
982 std::string tempExpResScale =
983 getParameter<std::string> (*params_,
"Explicit Residual Scaling");
986 if (expResScale_ != tempExpResScale) {
988 expResScale_ = tempExpResScale;
990 if (! expConvTest_.is_null()) {
1008 if (maxIterTest_.is_null()) {
1011 maxIterTest_->setMaxIters (maxIters_);
1016 if (impConvTest_.is_null()) {
1017 impConvTest_ =
rcp (
new StatusTestResNorm_t (convTol_));
1023 if (expConvTest_.is_null()) {
1024 expConvTest_ =
rcp (
new StatusTestResNorm_t (convTol_));
1025 expConvTest_->defineResForm (StatusTestResNorm_t::Explicit,
Belos::TwoNorm);
1031 if (convTest_.is_null()) {
1032 convTest_ =
rcp (
new StatusTestCombo_t (StatusTestCombo_t::SEQ,
1040 sTest_ =
rcp (
new StatusTestCombo_t (StatusTestCombo_t::OR,
1041 maxIterTest_, convTest_));
1045 outputTest_ = stoFactory.
create (printer_, sTest_, outputFreq_,
1049 std::string solverDesc =
"Block GCRODR ";
1050 outputTest_->setSolverDesc (solverDesc);
1057 template<
class ScalarType,
class MV,
class OP>
1068 int KrylSpaDim = (numBlocks_ - 1) * blockSize_;
1071 int augSpaDim = KrylSpaDim + recycledBlocks_ + 1;
1074 int augSpaImgDim = KrylSpaDim + blockSize_ + recycledBlocks_+1;
1098 if (U_ == Teuchos::null) {
1099 U_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1103 if (MVT::GetNumberVecs(*U_) < recycledBlocks_+1) {
1105 U_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1110 if (C_ == Teuchos::null) {
1111 C_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1115 if (MVT::GetNumberVecs(*C_) < recycledBlocks_+1) {
1117 C_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1122 if (U1_ == Teuchos::null) {
1123 U1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1127 if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1129 U1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1134 if (C1_ == Teuchos::null) {
1135 C1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1139 if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1141 C1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1146 if (R_ == Teuchos::null){
1147 R_ = MVT::Clone( *rhsMV, blockSize_ );
1153 if (G_ == Teuchos::null){
1154 G_ =
Teuchos::rcp(
new SDM( augSpaImgDim, augSpaDim ) );
1157 if ( (G_->numRows() != augSpaImgDim) || (G_->numCols() != augSpaDim) )
1159 G_->reshape( augSpaImgDim, augSpaDim );
1161 G_->putScalar(zero);
1165 if (H_ == Teuchos::null){
1166 H_ =
Teuchos::rcp (
new SDM (
Teuchos::View, *G_, KrylSpaDim + blockSize_, KrylSpaDim, recycledBlocks_+1 ,recycledBlocks_+1 ) );
1170 if (F_ == Teuchos::null){
1171 F_ =
Teuchos::rcp(
new SDM( recycledBlocks_+1, recycledBlocks_+1 ) );
1174 if ( (F_->numRows() != recycledBlocks_+1) || (F_->numCols() != recycledBlocks_+1) ){
1175 F_->reshape( recycledBlocks_+1, recycledBlocks_+1 );
1178 F_->putScalar(zero);
1181 if (PP_ == Teuchos::null){
1182 PP_ =
Teuchos::rcp(
new SDM( augSpaImgDim, recycledBlocks_+1 ) );
1185 if ( (PP_->numRows() != augSpaImgDim) || (PP_->numCols() != recycledBlocks_+1) ){
1186 PP_->reshape( augSpaImgDim, recycledBlocks_+1 );
1191 if (HP_ == Teuchos::null)
1192 HP_ =
Teuchos::rcp(
new SDM( augSpaImgDim, augSpaDim ) );
1194 if ( (HP_->numRows() != augSpaImgDim) || (HP_->numCols() != augSpaDim) ){
1195 HP_->reshape( augSpaImgDim, augSpaDim );
1200 tau_.resize(recycledBlocks_+1);
1203 work_.resize(recycledBlocks_+1);
1206 ipiv_.resize(recycledBlocks_+1);
1210 template<
class ScalarType,
class MV,
class OP>
1211 void BlockGCRODRSolMgr<ScalarType,MV,OP>::buildRecycleSpaceKryl(
int& keff,
Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > block_gmres_iter){
1216 int p = block_gmres_iter->getState().curDim;
1217 std::vector<int> index(keff);
1222 if(recycledBlocks_ >= p + blockSize_){
1226 for (
int ii=0; ii < p; ++ii) index[ii] = ii;
1228 MVT::SetBlock(*V_, index, *Utmp);
1234 if(recycleMethod_ ==
"harmvecs"){
1235 keff = getHarmonicVecsKryl(p, HH, *PPtmp);
1236 printer_->stream(
Debug) <<
"keff = " << keff << std::endl;
1242 for (
int ii=0; ii<keff; ++ii) index[ii] = ii;
1247 for (
int ii=0; ii < p; ++ii) index[ii] = ii;
1252 MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
1260 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1261 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
1265 work_.resize(lwork);
1266 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1267 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
1272 for(
int ii=0;ii<keff;ii++) {
for(
int jj=ii;jj<keff;jj++) Rtmp(ii,jj) = HPtmp(ii,jj); }
1274 lapack.UNGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1275 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _UNGQR failed to construct the Q factor.");
1279 index.resize( p+blockSize_ );
1280 for (
int ii=0; ii < (p+blockSize_); ++ii) { index[ii] = ii; }
1281 Vtmp = MVT::CloneView( *V_, index );
1282 MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
1289 ipiv_.resize(Rtmp.numRows());
1290 lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
1291 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1293 lwork = Rtmp.numRows();
1294 work_.resize(lwork);
1295 lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1298 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
1304 template<
class ScalarType,
class MV,
class OP>
1305 void BlockGCRODRSolMgr<ScalarType,MV,OP>::buildRecycleSpaceAugKryl(
Teuchos::RCP<BlockGCRODRIter<ScalarType,MV,OP> > block_gcrodr_iter){
1309 std::vector<MagnitudeType> d(keff);
1310 std::vector<ScalarType> dscalar(keff);
1311 std::vector<int> index(numBlocks_+1);
1314 BlockGCRODRIterState<ScalarType,MV> oldState = block_gcrodr_iter->getState();
1315 int p = oldState.curDim;
1320 if(recycledBlocks_ >= keff + p){
1323 for (
int ii=0; ii < p; ++ii) { index[ii] = keff+ii; }
1325 for (
int ii=0; ii < p; ++ii) { index[ii] =ii; }
1326 MVT::SetBlock(*V_, index, *Utmp);
1333 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1336 dscalar.resize(keff);
1337 MVT::MvNorm( *Utmp, d );
1338 for (
int i=0; i<keff; ++i) {
1340 dscalar[i] = (ScalarType)d[i];
1342 MVT::MvScale( *Utmp, dscalar );
1350 for (
int i=0; i<keff; ++i)
1351 (*Gtmp)(i,i) = d[i];
1358 SDM PPtmp(
Teuchos::View, *PP_, p+keff-blockSize_, recycledBlocks_+1 );
1359 keff_new = getHarmonicVecsAugKryl( keff, p-blockSize_, *Gtmp, oldState.V, PPtmp );
1368 index.resize( keff );
1369 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1371 index.resize( keff_new );
1372 for (
int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
1373 U1tmp = MVT::CloneViewNonConst( *U1_, index );
1375 MVT::MvTimesMatAddMv( one, *Utmp, PPtmp, zero, *U1tmp );
1380 index.resize(p-blockSize_);
1381 for (
int ii=0; ii < p-blockSize_; ii++) { index[ii] = ii; }
1383 SDM PPtmp(
Teuchos::View, *PP_, p-blockSize_, keff_new, keff );
1384 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
1390 SDM PPtmp(
Teuchos::View, *PP_, p-blockSize_+keff, keff_new );
1395 int info = 0, lwork = -1;
1396 tau_.resize(keff_new);
1397 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1398 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
1401 work_.resize(lwork);
1402 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1403 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
1408 for(
int i=0;i<keff_new;i++) {
for(
int j=i;j<keff_new;j++) Ftmp(i,j) = HPtmp(i,j); }
1410 lapack.UNGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1411 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _UNGQR failed to construct the Q factor.");
1421 for (
int i=0; i < keff; i++) { index[i] = i; }
1423 index.resize(keff_new);
1424 for (
int i=0; i < keff_new; i++) { index[i] = i; }
1425 C1tmp = MVT::CloneViewNonConst( *C1_, index );
1427 MVT::MvTimesMatAddMv( one, *Ctmp, PPtmp, zero, *C1tmp );
1432 for (
int i=0; i < p; ++i) { index[i] = i; }
1435 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *C1tmp );
1444 ipiv_.resize(Ftmp.numRows());
1445 lapack.GETRF(Ftmp.numRows(),Ftmp.numCols(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&info);
1446 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1449 lwork = Ftmp.numRows();
1450 work_.resize(lwork);
1451 lapack.GETRI(Ftmp.numRows(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1452 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRI failed to compute an LU factorization.");
1455 index.resize(keff_new);
1456 for (
int i=0; i < keff_new; i++) { index[i] = i; }
1458 MVT::MvTimesMatAddMv( one, *U1tmp, Ftmp, zero, *Utmp );
1463 if (keff != keff_new) {
1465 block_gcrodr_iter->setSize( keff, numBlocks_ );
1467 SDM b1(
Teuchos::View, *G_, recycledBlocks_+2, 1, 0, recycledBlocks_ );
1473 template<
class ScalarType,
class MV,
class OP>
1474 int BlockGCRODRSolMgr<ScalarType,MV,OP>::getHarmonicVecsAugKryl(
int keff,
int m,
const SDM& GG,
const Teuchos::RCP<const MV>& VV, SDM& PP){
1476 int m2 = GG.numCols();
1477 bool xtraVec =
false;
1480 std::vector<int> index;
1483 std::vector<MagnitudeType> wr(m2), wi(m2);
1486 std::vector<MagnitudeType> w(m2);
1489 SDM vr(m2,m2,
false);
1492 std::vector<int> iperm(m2);
1495 builtRecycleSpace_ =
true;
1505 SDM A_tmp( keff+m+blockSize_, keff+m );
1510 for (
int i=0; i<keff; ++i) { index[i] = i; }
1514 MVT::MvTransMv( one, *Ctmp, *Utmp, A11 );
1518 index.resize(m+blockSize_);
1519 for (i=0; i < m+blockSize_; i++) { index[i] = i; }
1521 MVT::MvTransMv( one, *Vp, *Utmp, A21 );
1524 for( i=keff; i<keff+m; i++)
1528 SDM A( m2, A_tmp.numCols() );
1537 char balanc=
'P', jobvl=
'N', jobvr=
'V', sense=
'N';
1538 int ld = A.numRows();
1540 int ldvl = ld, ldvr = ld;
1541 int info = 0,ilo = 0,ihi = 0;
1542 MagnitudeType abnrm = 0.0, bbnrm = 0.0;
1544 std::vector<ScalarType> beta(ld);
1545 std::vector<ScalarType> work(lwork);
1546 std::vector<MagnitudeType> rwork(lwork);
1547 std::vector<MagnitudeType> lscale(ld), rscale(ld);
1548 std::vector<MagnitudeType> rconde(ld), rcondv(ld);
1549 std::vector<int> iwork(ld+6);
1554 lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0],
1555 &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
1556 &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &rwork[0],
1557 &iwork[0], bwork, &info);
1558 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK GGEVX failed to compute eigensolutions.");
1562 for( i=0; i<ld; i++ )
1565 this->sort(w,ld,iperm);
1570 for( i=0; i<recycledBlocks_; i++ )
1571 for( j=0; j<ld; j++ )
1572 PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
1574 if(scalarTypeIsComplex==
false) {
1577 if (wi[iperm[ld-recycledBlocks_]] != 0.0) {
1579 for ( i=ld-recycledBlocks_; i<ld; i++ )
1580 if (wi[iperm[i]] != 0.0) countImag++;
1582 if (countImag % 2) xtraVec =
true;
1586 if (wi[iperm[ld-recycledBlocks_]] > 0.0) {
1587 for( j=0; j<ld; j++ )
1588 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
1591 for( j=0; j<ld; j++ )
1592 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
1600 return recycledBlocks_+1;
1602 return recycledBlocks_;
1606 template<
class ScalarType,
class MV,
class OP>
1607 int BlockGCRODRSolMgr<ScalarType,MV,OP>::getHarmonicVecsKryl(
int m,
const SDM& HH, SDM& PP){
1608 bool xtraVec =
false;
1613 std::vector<MagnitudeType> wr(m), wi(m);
1619 std::vector<MagnitudeType> w(m);
1622 std::vector<int> iperm(m);
1626 std::vector<ScalarType> work(lwork);
1627 std::vector<MagnitudeType> rwork(lwork);
1633 builtRecycleSpace_ =
true;
1640 for(
int i=0; i<=blockSize_-1; i++) (*harmRitzMatrix)[blockSize_-1-i][harmRitzMatrix->numRows()-1-i] = 1;
1643 lapack.GESV(m, blockSize_, HHt.values(), HHt.stride(), &iperm[0], harmRitzMatrix->values(), harmRitzMatrix->stride(), &info);
1645 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK GESV failed to compute a solution.");
1656 Htemp =
Teuchos::rcp(
new SDM(H_lbl_t.numRows(), H_lbl_t.numCols()));
1658 H_lbl_t.assign(*Htemp);
1660 Htemp =
Teuchos::rcp(
new SDM(harmRitzMatrix -> numRows(), harmRitzMatrix -> numCols()));
1662 harmRitzMatrix -> assign(*Htemp);
1665 int harmColIndex, HHColIndex;
1667 for(
int i = 0; i<blockSize_; i++){
1668 harmColIndex = harmRitzMatrix -> numCols() - i -1;
1670 for(
int j=0; j<m; j++) (*Htemp)[HHColIndex][j] += (*harmRitzMatrix)[harmColIndex][j];
1672 harmRitzMatrix = Htemp;
1682 lapack.GEEV(
'N',
'V', m, harmRitzMatrix -> values(), harmRitzMatrix -> stride(), &wr[0], &wi[0],
1683 vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
1685 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK GEEV failed to compute eigensolutions.");
1688 for(
int i=0; i<m; ++i ) w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] );
1690 this->sort(w, m, iperm);
1695 for(
int i=0; i<recycledBlocks_; ++i )
1696 for(
int j=0; j<m; j++ )
1697 PP(j,i) = vr(j,iperm[i]);
1699 if(scalarTypeIsComplex==
false) {
1702 if (wi[iperm[recycledBlocks_-1]] != 0.0) {
1704 for (
int i=0; i<recycledBlocks_; ++i )
1705 if (wi[iperm[i]] != 0.0) countImag++;
1707 if (countImag % 2) xtraVec =
true;
1711 if (wi[iperm[recycledBlocks_-1]] > 0.0) {
1712 for(
int j=0; j<m; ++j )
1713 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
1716 for(
int j=0; j<m; ++j )
1717 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
1725 printer_->stream(
Debug) <<
"Recycled " << recycledBlocks_+1 <<
" vectors" << std::endl;
1726 return recycledBlocks_+1;
1729 printer_->stream(
Debug) <<
"Recycled " << recycledBlocks_ <<
" vectors" << std::endl;
1730 return recycledBlocks_;
1735 template<
class ScalarType,
class MV,
class OP>
1736 void BlockGCRODRSolMgr<ScalarType,MV,OP>::sort(std::vector<MagnitudeType>& dlist,
int n, std::vector<int>& iperm) {
1737 int l, r, j, i, flag;
1739 MagnitudeType dRR, dK;
1764 if (dlist[j] > dlist[j - 1]) j = j + 1;
1765 if (dlist[j - 1] > dK) {
1766 dlist[i - 1] = dlist[j - 1];
1767 iperm[i - 1] = iperm[j - 1];
1780 dlist[r] = dlist[0];
1781 iperm[r] = iperm[0];
1795 template<
class ScalarType,
class MV,
class OP>
1799 using Teuchos::rcp_const_cast;
1805 std::vector<int> index(numBlocks_+1);
1821 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1822 int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1826 std::vector<int> currIdx;
1828 if ( adaptiveBlockSize_ ) {
1829 blockSize_ = numCurrRHS;
1830 currIdx.resize( numCurrRHS );
1831 for (
int i=0; i<numCurrRHS; ++i)
1832 currIdx[i] = startPtr+i;
1835 currIdx.resize( blockSize_ );
1836 for (
int i=0; i<numCurrRHS; ++i)
1837 currIdx[i] = startPtr+i;
1838 for (
int i=numCurrRHS; i<blockSize_; ++i)
1843 problem_->setLSIndex( currIdx );
1849 loaDetected_ =
false;
1852 bool isConverged =
true;
1855 initializeStateStorage();
1860 while (numRHS2Solve > 0){
1870 printer_->stream(
Debug) <<
" Now solving RHS index " << currIdx[0] <<
" using recycled subspace of dimension " << keff << std::endl << std::endl;
1874 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1875 RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1876 RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1877 problem_->apply( *Utmp, *Ctmp );
1879 RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1884 int rank = ortho_->normalize(*Ctmp,
rcp(&Ftmp,
false));
1897 work_.resize(lwork);
1902 MVT::MvTimesMatAddMv( one, *Utmp, Ftmp, zero, *U1tmp );
1907 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1908 Ctmp = MVT::CloneViewNonConst( *C_, index );
1909 Utmp = MVT::CloneView( *U_, index );
1912 SDM Ctr(keff,blockSize_);
1913 problem_->computeCurrPrecResVec( &*R_ );
1914 MVT::MvTransMv( one, *Ctmp, *R_, Ctr );
1917 RCP<MV> update = MVT::Clone( *problem_->getCurrLHSVec(), blockSize_ );
1918 MVT::MvInit( *update, 0.0 );
1919 MVT::MvTimesMatAddMv( one, *Utmp, Ctr, one, *update );
1920 problem_->updateSolution( update,
true );
1923 MVT::MvTimesMatAddMv( -one, *Ctmp, Ctr, one, *R_ );
1931 if (V_ == Teuchos::null) {
1934 V_ = MVT::Clone( *rhsMV, (numBlocks_+1)*blockSize_ );
1938 if (MVT::GetNumberVecs(*V_) < (numBlocks_+1)*blockSize_ ) {
1940 V_ = MVT::Clone( *tmp, (numBlocks_+1)*blockSize_ );
1945 printer_->stream(
Debug) <<
" No recycled subspace available for RHS index " << std::endl << std::endl;
1950 primeList.
set(
"Num Blocks",numBlocks_-1);
1951 primeList.
set(
"Block Size",blockSize_);
1952 primeList.
set(
"Recycled Blocks",0);
1953 primeList.
set(
"Keep Hessenberg",
true);
1954 primeList.
set(
"Initialize Hessenberg",
true);
1956 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1957 if (blockSize_*static_cast<ptrdiff_t>(numBlocks_) > dim) {
1958 ptrdiff_t tmpNumBlocks = 0;
1959 if (blockSize_ == 1)
1960 tmpNumBlocks = dim / blockSize_;
1962 tmpNumBlocks = ( dim - blockSize_) / blockSize_;
1963 printer_->stream(
Warnings) <<
"Belos::BlockGmresSolMgr::solve(): Warning! Requested Krylov subspace dimension is larger than operator dimension!" 1964 << std::endl <<
"The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << tmpNumBlocks << std::endl;
1965 primeList.
set(
"Num Blocks",Teuchos::as<int>(tmpNumBlocks));
1969 primeList.
set(
"Num Blocks",numBlocks_-1);
1976 block_gmres_iter->setSize( blockSize_, numBlocks_-1 );
1980 if (currIdx[blockSize_-1] == -1) {
1981 V_0 = MVT::Clone( *(problem_->getInitPrecResVec()), blockSize_ );
1982 problem_->computeCurrPrecResVec( &*V_0 );
1985 V_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
1992 int rank = ortho_->normalize( *V_0, z_0 );
2001 block_gmres_iter->initializeGmres(newstate);
2003 bool primeConverged =
false;
2006 printer_->stream(
Debug) <<
" Preparing to Iterate!!!!" << std::endl << std::endl;
2007 block_gmres_iter->iterate();
2012 if ( convTest_->getStatus() ==
Passed ) {
2013 printer_->stream(
Debug) <<
"We converged during the prime the pump stage" << std::endl << std::endl;
2014 primeConverged = !(expConvTest_->getLOADetected());
2015 if ( expConvTest_->getLOADetected() ) {
2017 loaDetected_ =
true;
2018 printer_->stream(
Warnings) <<
"Belos::BlockGmresSolMgr::solve(): Warning! Solver has experienced a loss of accuracy!" << std::endl;
2024 else if( maxIterTest_->getStatus() ==
Passed ) {
2026 primeConverged =
false;
2032 printer_->stream(
Debug) <<
" We did not converge on priming cycle of Block GMRES. Therefore we recycle and restart. " << std::endl << std::endl;
2037 if (blockSize_ != 1) {
2038 printer_->stream(
Errors) <<
"Error! Caught std::exception in BlockGmresIter::iterate() at iteration " 2039 << block_gmres_iter->getNumIters() << std::endl
2040 << e.what() << std::endl;
2041 if (convTest_->getStatus() !=
Passed)
2042 primeConverged =
false;
2046 block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
2048 sTest_->checkStatus( &*block_gmres_iter );
2049 if (convTest_->getStatus() !=
Passed)
2050 isConverged =
false;
2053 catch (
const std::exception &e) {
2054 printer_->stream(
Errors) <<
"Error! Caught std::exception in BlockGmresIter::iterate() at iteration " 2055 << block_gmres_iter->getNumIters() << std::endl
2056 << e.what() << std::endl;
2064 RCP<MV> update = block_gmres_iter->getCurrentUpdate();
2065 problem_->updateSolution( update,
true );
2069 problem_->computeCurrPrecResVec( &*R_ );
2072 newstate = block_gmres_iter->getState();
2075 H_->assign(*(newstate.
H));
2084 V_ = rcp_const_cast<MV>(newstate.
V);
2087 buildRecycleSpaceKryl(keff, block_gmres_iter);
2088 printer_->stream(
Debug) <<
"Generated recycled subspace using RHS index " << currIdx[0] <<
" of dimension " << keff << std::endl << std::endl;
2092 if (primeConverged) {
2113 problem_->setCurrLS();
2116 startPtr += numCurrRHS;
2117 numRHS2Solve -= numCurrRHS;
2118 if ( numRHS2Solve > 0 ) {
2119 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
2120 if ( adaptiveBlockSize_ ) {
2121 blockSize_ = numCurrRHS;
2122 currIdx.resize( numCurrRHS );
2123 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2126 currIdx.resize( blockSize_ );
2127 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2128 for (
int i=numCurrRHS; i<blockSize_; ++i) currIdx[i] = -1;
2131 problem_->setLSIndex( currIdx );
2134 currIdx.resize( numRHS2Solve );
2145 blockgcrodrList.
set(
"Num Blocks",numBlocks_);
2146 blockgcrodrList.
set(
"Block Size",blockSize_);
2147 blockgcrodrList.
set(
"Recycled Blocks",keff);
2153 index.resize( blockSize_ );
2154 for(
int ii = 0; ii < blockSize_; ii++) index[ii] = ii;
2158 MVT::Assign(*R_,*V0);
2161 for(
int i=0; i < keff; i++){ index[i] = i;};
2163 H_ =
rcp(
new SDM(
Teuchos::View, *G_, (numBlocks_-1)*blockSize_ + blockSize_, (numBlocks_-1)*blockSize_, keff ,keff ));
2167 newstate.
U = MVT::CloneViewNonConst(*U_, index);
2168 newstate.
C = MVT::CloneViewNonConst(*C_, index);
2170 newstate.
curDim = blockSize_;
2171 block_gcrodr_iter -> initialize(newstate);
2173 int numRestarts = 0;
2177 block_gcrodr_iter -> iterate();
2182 if( convTest_->getStatus() ==
Passed ) {
2190 else if(maxIterTest_->getStatus() ==
Passed ){
2192 isConverged =
false;
2199 else if (block_gcrodr_iter->getCurSubspaceDim() == block_gcrodr_iter->getMaxSubspaceDim()){
2205 problem_->updateSolution(update,
true);
2206 buildRecycleSpaceAugKryl(block_gcrodr_iter);
2208 printer_->stream(
Debug) <<
" Generated new recycled subspace using RHS index " << currIdx[0] <<
" of dimension " << keff << std::endl << std::endl;
2210 if(numRestarts >= maxRestarts_) {
2211 isConverged =
false;
2217 printer_ -> stream(
Debug) <<
" Performing restart number " << numRestarts <<
" of " << maxRestarts_ << std::endl << std::endl;
2220 problem_->computeCurrPrecResVec( &*R_ );
2221 index.resize( blockSize_ );
2222 for (
int ii=0; ii<blockSize_; ++ii) index[ii] = ii;
2224 MVT::SetBlock(*R_,index,*V0);
2228 index.resize( numBlocks_*blockSize_ );
2229 for (
int ii=0; ii<(numBlocks_*blockSize_); ++ii) index[ii] = ii;
2230 restartState.
V = MVT::CloneViewNonConst( *V_, index );
2231 index.resize( keff );
2232 for (
int ii=0; ii<keff; ++ii) index[ii] = ii;
2233 restartState.
U = MVT::CloneViewNonConst( *U_, index );
2234 restartState.
C = MVT::CloneViewNonConst( *C_, index );
2236 H_ =
rcp(
new SDM(
Teuchos::View, *G_, numBlocks_*blockSize_, (numBlocks_-1)*blockSize_, keff ,keff ));
2237 restartState.
B = B_;
2238 restartState.
H = H_;
2239 restartState.
curDim = blockSize_;
2240 block_gcrodr_iter->initialize(restartState);
2249 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Belos::BlockGCRODRSolMgr::solve(): Invalid return from BlockGCRODRIter::iterate().");
2255 block_gcrodr_iter->updateLSQR( block_gcrodr_iter->getCurSubspaceDim() );
2257 sTest_->checkStatus( &*block_gcrodr_iter );
2258 if (convTest_->getStatus() !=
Passed) isConverged =
false;
2261 catch(
const std::exception &e){
2262 printer_->stream(
Errors) <<
"Error! Caught exception in BlockGCRODRIter::iterate() at iteration " 2263 << block_gcrodr_iter->getNumIters() << std::endl
2264 << e.what() << std::endl;
2272 problem_->updateSolution( update,
true );
2291 problem_->setCurrLS();
2294 startPtr += numCurrRHS;
2295 numRHS2Solve -= numCurrRHS;
2296 if ( numRHS2Solve > 0 ) {
2297 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
2298 if ( adaptiveBlockSize_ ) {
2299 blockSize_ = numCurrRHS;
2300 currIdx.resize( numCurrRHS );
2301 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2304 currIdx.resize( blockSize_ );
2305 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2306 for (
int i=numCurrRHS; i<blockSize_; ++i) currIdx[i] = -1;
2309 problem_->setLSIndex( currIdx );
2312 currIdx.resize( numRHS2Solve );
2316 if (!builtRecycleSpace_) {
2317 buildRecycleSpaceAugKryl(block_gcrodr_iter);
2318 printer_->stream(
Debug) <<
" Generated new recycled subspace using RHS index " << currIdx[0] <<
" of dimension " << keff << std::endl << std::endl;
2326 #ifdef BELOS_TEUCHOS_TIME_MONITOR 2332 numIters_ = maxIterTest_->getNumIters();
2335 const std::vector<MagnitudeType>* pTestValues = NULL;
2336 pTestValues = impConvTest_->getTestValue();
2338 "Belos::BlockGCRODRSolMgr::solve(): The implicit convergence test's " 2339 "getTestValue() method returned NULL. Please report this bug to the " 2340 "Belos developers.");
2342 "Belos::BlockGCRODRSolMgr::solve(): The implicit convergence test's " 2343 "getTestValue() method returned a vector of length zero. Please report " 2344 "this bug to the Belos developers.");
2348 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.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
Class which manages the output and verbosity of the Belos solvers.
Structure to contain pointers to BlockGCRODRIter state variables.
ScalarType * values() const
virtual ~BlockGCRODRSolMgr()
Destructor.
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().
This class implements the block GMRES iteration, where a block Krylov subspace is constructed...
Teuchos::RCP< const MV > V
The current Krylov basis.
T & get(const std::string &name, T def_value)
BlockGCRODRIterOrthoFailure is thrown when the BlockGCRODRIter object is unable to compute independen...
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Thrown when the linear problem was not set up correctly.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
A factory class for generating StatusTestOutput objects.
static magnitudeType real(T a)
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Set the parameters the solver should use to solve the linear problem.
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.
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< 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.
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
BlockGCRODRSolMgrLAPACKFailure(const std::string &what_arg)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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. ...
OrdinalType numRows() const
int curDim
The current dimension of the reduction.
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
Belos concrete class for performing the block GCRO-DR (block GMRES with recycling) iteration...
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().
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
MagnitudeType achievedTol() const
Get the residual for the most recent call to solve().
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::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
Belos concrete class for performing the block GMRES iteration.
Implementation of the Block GCRO-DR (Block Recycling GMRES) iteration.
static magnitudeType magnitude(T a)
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.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
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.
Class which defines basic traits for the operator type.
OrdinalType stride() const
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 ...
OrdinalType numCols() const
Belos header file which uses auto-configuration information to include necessary C++ headers...
BlockGCRODRSolMgrLinearProblemFailure(const std::string &what_arg)
ReturnType solve()
Solve the current linear problem.
Thrown when an LAPACK call fails.
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.