45 #ifndef BELOS_BLOCK_GCRODR_SOLMGR_HPP
46 #define BELOS_BLOCK_GCRODR_SOLMGR_HPP
64 #ifdef BELOS_TEUCHOS_TIME_MONITOR
66 #endif // BELOS_TEUCHOS_TIME_MONITOR
125 template<
class ScalarType,
class MV,
class OP>
224 "Belos::BlockGCRODRSolMgr::setProblem: The input LinearProblem cannot be null.");
227 if (! problem->isProblemSet()) {
228 const bool success = problem->setProblem();
230 "Belos::BlockGCRODRSolMgr::setProblem: Calling the input LinearProblem's setProblem() method failed. This likely means that the "
231 "LinearProblem has a missing (null) matrix A, solution vector X, or right-hand side vector B. Please set these items in the LinearProblem and try again.");
321 void sort (std::vector<MagnitudeType>& dlist,
int n, std::vector<int>& iperm);
426 template<
class ScalarType,
class MV,
class OP>
429 template<
class ScalarType,
class MV,
class OP>
436 template<
class ScalarType,
class MV,
class OP>
442 template<
class ScalarType,
class MV,
class OP>
451 "Belos::BlockGCRODR constructor: The solver manager's constructor needs "
452 "the linear problem argument 'problem' to be nonnull.");
462 template<
class ScalarType,
class MV,
class OP>
464 adaptiveBlockSize_ = adaptiveBlockSize_default_;
465 recycleMethod_ = recycleMethod_default_;
467 loaDetected_ =
false;
468 builtRecycleSpace_ =
false;
539 template<
class ScalarType,
class MV,
class OP>
541 std::ostringstream oss;
542 oss <<
"Belos::BlockGCRODRSolMgr<" << SCT::name() <<
", ...>";
544 oss <<
"Ortho Type='"<<orthoType_ ;
545 oss <<
", Num Blocks=" <<numBlocks_;
546 oss <<
", Block Size=" <<blockSize_;
547 oss <<
", Num Recycle Blocks=" << recycledBlocks_;
548 oss <<
", Max Restarts=" << maxRestarts_;
553 template<
class ScalarType,
class MV,
class OP>
557 using Teuchos::parameterList;
559 using Teuchos::rcpFromRef;
561 if (defaultParams_.is_null()) {
562 RCP<ParameterList> pl = parameterList ();
564 const MagnitudeType convTol = SMT::squareroot (SCT::magnitude (SCT::eps()));
565 const int maxRestarts = 1000;
566 const int maxIters = 5000;
567 const int blockSize = 2;
568 const int numBlocks = 100;
569 const int numRecycledBlocks = 25;
572 const int outputFreq = 1;
573 RCP<std::ostream> outputStream = rcpFromRef (std::cout);
574 const std::string impResScale (
"Norm of Preconditioned Initial Residual");
575 const std::string expResScale (
"Norm of Initial Residual");
576 const std::string timerLabel (
"Belos");
577 const std::string orthoType (
"ICGS");
578 RCP<const ParameterList> orthoParams = orthoFactory_.getDefaultParameters (orthoType);
588 pl->set (
"Convergence Tolerance", convTol,
589 "The tolerance that the solver needs to achieve in order for "
590 "the linear system(s) to be declared converged. The meaning "
591 "of this tolerance depends on the convergence test details.");
592 pl->set(
"Maximum Restarts", maxRestarts,
593 "The maximum number of restart cycles (not counting the first) "
594 "allowed for each set of right-hand sides solved.");
595 pl->set(
"Maximum Iterations", maxIters,
596 "The maximum number of iterations allowed for each set of "
597 "right-hand sides solved.");
598 pl->set(
"Block Size", blockSize,
599 "How many linear systems to solve at once.");
600 pl->set(
"Num Blocks", numBlocks,
601 "The maximum number of blocks allowed in the Krylov subspace "
602 "for each set of right-hand sides solved. This includes the "
603 "initial block. Total Krylov basis storage, not counting the "
604 "recycled subspace, is \"Num Blocks\" times \"Block Size\".");
605 pl->set(
"Num Recycled Blocks", numRecycledBlocks,
606 "The maximum number of vectors in the recycled subspace." );
607 pl->set(
"Verbosity", verbosity,
608 "What type(s) of solver information should be written "
609 "to the output stream.");
610 pl->set(
"Output Style", outputStyle,
611 "What style is used for the solver information to write "
612 "to the output stream.");
613 pl->set(
"Output Frequency", outputFreq,
614 "How often convergence information should be written "
615 "to the output stream.");
616 pl->set(
"Output Stream", outputStream,
617 "A reference-counted pointer to the output stream where all "
618 "solver output is sent.");
619 pl->set(
"Implicit Residual Scaling", impResScale,
620 "The type of scaling used in the implicit residual convergence test.");
621 pl->set(
"Explicit Residual Scaling", expResScale,
622 "The type of scaling used in the explicit residual convergence test.");
623 pl->set(
"Timer Label", timerLabel,
624 "The string to use as a prefix for the timer labels.");
625 pl->set(
"Orthogonalization", orthoType,
626 "The orthogonalization method to use. Valid options: " +
627 orthoFactory_.validNamesString());
628 pl->set (
"Orthogonalization Parameters", *orthoParams,
629 "Sublist of parameters specific to the orthogonalization method to use.");
630 pl->set(
"Orthogonalization Constant", orthoKappa,
631 "When using DGKS orthogonalization: the \"depTol\" constant, used "
632 "to determine whether another step of classical Gram-Schmidt is "
633 "necessary. Otherwise ignored. Nonpositive values are ignored.");
636 return defaultParams_;
639 template<
class ScalarType,
class MV,
class OP>
643 using Teuchos::isParameterType;
644 using Teuchos::getParameter;
647 using Teuchos::parameterList;
650 using Teuchos::rcp_dynamic_cast;
651 using Teuchos::rcpFromRef;
657 RCP<const ParameterList> defaultParams = getValidParameters();
664 params_ = parameterList (*defaultParams);
692 const int maxRestarts = params_->
get<
int> (
"Maximum Restarts");
694 "Belos::BlockGCRODRSolMgr: The \"Maximum Restarts\" parameter "
695 "must be nonnegative, but you specified a negative value of "
696 << maxRestarts <<
".");
698 const int maxIters = params_->get<
int> (
"Maximum Iterations");
699 TEUCHOS_TEST_FOR_EXCEPTION(maxIters <= 0, std::invalid_argument,
700 "Belos::BlockGCRODRSolMgr: The \"Maximum Iterations\" parameter "
701 "must be positive, but you specified a nonpositive value of "
704 const int numBlocks = params_->get<
int> (
"Num Blocks");
705 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument,
706 "Belos::BlockGCRODRSolMgr: The \"Num Blocks\" parameter must be "
707 "positive, but you specified a nonpositive value of " << numBlocks
710 const int blockSize = params_->get<
int> (
"Block Size");
711 TEUCHOS_TEST_FOR_EXCEPTION(blockSize <= 0, std::invalid_argument,
712 "Belos::BlockGCRODRSolMgr: The \"Block Size\" parameter must be "
713 "positive, but you specified a nonpositive value of " << blockSize
716 const int recycledBlocks = params_->get<
int> (
"Num Recycled Blocks");
717 TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks <= 0, std::invalid_argument,
718 "Belos::BlockGCRODRSolMgr: The \"Num Recycled Blocks\" parameter must "
719 "be positive, but you specified a nonpositive value of "
720 << recycledBlocks <<
".");
721 TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks >= numBlocks,
722 std::invalid_argument,
"Belos::BlockGCRODRSolMgr: The \"Num Recycled "
723 "Blocks\" parameter must be less than the \"Num Blocks\" parameter, "
724 "but you specified \"Num Recycled Blocks\" = " << recycledBlocks
725 <<
" and \"Num Blocks\" = " << numBlocks <<
".");
729 maxRestarts_ = maxRestarts;
730 maxIters_ = maxIters;
731 numBlocks_ = numBlocks;
732 blockSize_ = blockSize;
733 recycledBlocks_ = recycledBlocks;
740 std::string tempLabel = params_->get<std::string> (
"Timer Label");
741 const bool labelChanged = (tempLabel != label_);
743 #ifdef BELOS_TEUCHOS_TIME_MONITOR
744 std::string solveLabel = label_ +
": BlockGCRODRSolMgr total solve time";
745 if (timerSolve_.is_null()) {
748 }
else if (labelChanged) {
757 #endif // BELOS_TEUCHOS_TIME_MONITOR
763 if (params_->isParameter (
"Verbosity")) {
764 if (isParameterType<int> (*params_,
"Verbosity")) {
765 verbosity_ = params_->get<
int> (
"Verbosity");
768 verbosity_ = (int) getParameter<MsgType> (*params_,
"Verbosity");
775 if (params_->isParameter (
"Output Style")) {
776 if (isParameterType<int> (*params_,
"Output Style")) {
777 outputStyle_ = params_->get<
int> (
"Output Style");
780 outputStyle_ = (int) getParameter<OutputType> (*params_,
"Output Style");
808 if (params_->isParameter (
"Output Stream")) {
810 outputStream_ = getParameter<RCP<std::ostream> > (*params_,
"Output Stream");
812 catch (InvalidParameter&) {
813 outputStream_ = rcpFromRef (std::cout);
820 if (outputStream_.is_null()) {
825 outputFreq_ = params_->get<
int> (
"Output Frequency");
828 if (! outputTest_.is_null()) {
829 outputTest_->setOutputFrequency (outputFreq_);
837 if (printer_.is_null()) {
840 printer_->setVerbosity (verbosity_);
841 printer_->setOStream (outputStream_);
852 if (params_->isParameter (
"Orthogonalization")) {
853 const std::string& tempOrthoType =
854 params_->get<std::string> (
"Orthogonalization");
856 if (! orthoFactory_.isValidName (tempOrthoType)) {
857 std::ostringstream os;
858 os <<
"Belos::BlockGCRODRSolMgr: Invalid orthogonalization name \""
859 << tempOrthoType <<
"\". The following are valid options "
860 <<
"for the \"Orthogonalization\" name parameter: ";
861 orthoFactory_.printValidNames (os);
864 if (tempOrthoType != orthoType_) {
866 orthoType_ = tempOrthoType;
883 RCP<ParameterList> orthoParams = sublist (params,
"Orthogonalization Parameters",
true);
885 "Failed to get orthogonalization parameters. "
886 "Please report this bug to the Belos developers.");
912 ortho_ = orthoFactory_.makeMatOrthoManager (orthoType_, null, printer_,
913 label_, orthoParams);
925 if (params_->isParameter (
"Orthogonalization Constant")) {
928 if (orthoKappa > 0) {
929 orthoKappa_ = orthoKappa;
932 if (orthoType_ ==
"DGKS" && ! ortho_.is_null()) {
938 rcp_dynamic_cast<ortho_man_type>(ortho_)->setDepTol (orthoKappa_);
948 convTol_ = params_->get<
MagnitudeType> (
"Convergence Tolerance");
949 if (! impConvTest_.is_null()) {
952 if (! expConvTest_.is_null()) {
953 expConvTest_->setTolerance (convTol_);
957 if (params_->isParameter (
"Implicit Residual Scaling")) {
958 std::string tempImpResScale =
959 getParameter<std::string> (*params_,
"Implicit Residual Scaling");
962 if (impResScale_ != tempImpResScale) {
964 impResScale_ = tempImpResScale;
966 if (! impConvTest_.is_null()) {
979 if (params_->isParameter(
"Explicit Residual Scaling")) {
980 std::string tempExpResScale =
981 getParameter<std::string> (*params_,
"Explicit Residual Scaling");
984 if (expResScale_ != tempExpResScale) {
986 expResScale_ = tempExpResScale;
988 if (! expConvTest_.is_null()) {
1006 if (maxIterTest_.is_null()) {
1009 maxIterTest_->setMaxIters (maxIters_);
1014 if (impConvTest_.is_null()) {
1015 impConvTest_ =
rcp (
new StatusTestResNorm_t (convTol_));
1021 if (expConvTest_.is_null()) {
1022 expConvTest_ =
rcp (
new StatusTestResNorm_t (convTol_));
1023 expConvTest_->defineResForm (StatusTestResNorm_t::Explicit,
Belos::TwoNorm);
1029 if (convTest_.is_null()) {
1030 convTest_ =
rcp (
new StatusTestCombo_t (StatusTestCombo_t::SEQ,
1038 sTest_ =
rcp (
new StatusTestCombo_t (StatusTestCombo_t::OR,
1039 maxIterTest_, convTest_));
1043 outputTest_ = stoFactory.
create (printer_, sTest_, outputFreq_,
1047 std::string solverDesc =
"Block GCRODR ";
1048 outputTest_->setSolverDesc (solverDesc);
1055 template<
class ScalarType,
class MV,
class OP>
1066 int KrylSpaDim = (numBlocks_ - 1) * blockSize_;
1069 int augSpaDim = KrylSpaDim + recycledBlocks_ + 1;
1072 int augSpaImgDim = KrylSpaDim + blockSize_ + recycledBlocks_+1;
1097 U_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1101 if (MVT::GetNumberVecs(*U_) < recycledBlocks_+1) {
1103 U_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1109 C_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1113 if (MVT::GetNumberVecs(*C_) < recycledBlocks_+1) {
1115 C_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1121 U1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1125 if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1127 U1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1133 C1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1137 if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1139 C1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1145 R_ = MVT::Clone( *rhsMV, blockSize_ );
1155 if ( (G_->numRows() != augSpaImgDim) || (G_->numCols() != augSpaDim) )
1157 G_->reshape( augSpaImgDim, augSpaDim );
1159 G_->putScalar(zero);
1172 if ( (F_->numRows() != recycledBlocks_+1) || (F_->numCols() != recycledBlocks_+1) ){
1173 F_->reshape( recycledBlocks_+1, recycledBlocks_+1 );
1176 F_->putScalar(zero);
1183 if ( (PP_->numRows() != augSpaImgDim) || (PP_->numCols() != recycledBlocks_+1) ){
1184 PP_->reshape( augSpaImgDim, recycledBlocks_+1 );
1192 if ( (HP_->numRows() != augSpaImgDim) || (HP_->numCols() != augSpaDim) ){
1193 HP_->reshape( augSpaImgDim, augSpaDim );
1198 tau_.resize(recycledBlocks_+1);
1201 work_.resize(recycledBlocks_+1);
1204 ipiv_.resize(recycledBlocks_+1);
1208 template<
class ScalarType,
class MV,
class OP>
1214 int p = block_gmres_iter->getState().curDim;
1215 std::vector<int> index(keff);
1220 if(recycledBlocks_ >= p + blockSize_){
1224 for (
int ii=0; ii < p; ++ii) index[ii] = ii;
1226 MVT::SetBlock(*V_, index, *Utmp);
1232 if(recycleMethod_ ==
"harmvecs"){
1233 keff = getHarmonicVecsKryl(p, HH, *PPtmp);
1234 printer_->stream(
Debug) <<
"keff = " << keff << std::endl;
1240 for (
int ii=0; ii<keff; ++ii) index[ii] = ii;
1245 for (
int ii=0; ii < p; ++ii) index[ii] = ii;
1250 MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
1263 work_.resize(lwork);
1270 for(
int ii=0;ii<keff;ii++) {
for(
int jj=ii;jj<keff;jj++) Rtmp(ii,jj) = HPtmp(ii,jj); }
1277 index.resize( p+blockSize_ );
1278 for (
int ii=0; ii < (p+blockSize_); ++ii) { index[ii] = ii; }
1279 Vtmp = MVT::CloneView( *V_, index );
1280 MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
1292 work_.resize(lwork);
1296 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
1302 template<
class ScalarType,
class MV,
class OP>
1307 std::vector<MagnitudeType> d(keff);
1308 std::vector<ScalarType> dscalar(keff);
1309 std::vector<int> index(numBlocks_+1);
1318 if(recycledBlocks_ >= keff + p){
1321 for (
int ii=0; ii < p; ++ii) { index[ii] = keff+ii; }
1323 for (
int ii=0; ii < p; ++ii) { index[ii] =ii; }
1324 MVT::SetBlock(*V_, index, *Utmp);
1331 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1334 dscalar.resize(keff);
1335 MVT::MvNorm( *Utmp, d );
1336 for (
int i=0; i<keff; ++i) {
1338 dscalar[i] = (ScalarType)d[i];
1340 MVT::MvScale( *Utmp, dscalar );
1348 for (
int i=0; i<keff; ++i)
1349 (*Gtmp)(i,i) = d[i];
1357 keff_new = getHarmonicVecsAugKryl( keff, p-blockSize_, *Gtmp, oldState.
V, PPtmp );
1366 index.resize( keff );
1367 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1369 index.resize( keff_new );
1370 for (
int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
1371 U1tmp = MVT::CloneViewNonConst( *U1_, index );
1373 MVT::MvTimesMatAddMv( one, *Utmp, PPtmp, zero, *U1tmp );
1378 index.resize(p-blockSize_);
1379 for (
int ii=0; ii < p-blockSize_; ii++) { index[ii] = ii; }
1382 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
1393 int info = 0, lwork = -1;
1394 tau_.resize(keff_new);
1399 work_.resize(lwork);
1406 for(
int i=0;i<keff_new;i++) {
for(
int j=i;j<keff_new;j++) Ftmp(i,j) = HPtmp(i,j); }
1419 for (
int i=0; i < keff; i++) { index[i] = i; }
1421 index.resize(keff_new);
1422 for (
int i=0; i < keff_new; i++) { index[i] = i; }
1423 C1tmp = MVT::CloneViewNonConst( *C1_, index );
1425 MVT::MvTimesMatAddMv( one, *Ctmp, PPtmp, zero, *C1tmp );
1430 for (
int i=0; i < p; ++i) { index[i] = i; }
1433 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *C1tmp );
1448 work_.resize(lwork);
1453 index.resize(keff_new);
1454 for (
int i=0; i < keff_new; i++) { index[i] = i; }
1456 MVT::MvTimesMatAddMv( one, *U1tmp, Ftmp, zero, *Utmp );
1461 if (keff != keff_new) {
1463 block_gcrodr_iter->setSize( keff, numBlocks_ );
1471 template<
class ScalarType,
class MV,
class OP>
1475 bool xtraVec =
false;
1478 std::vector<int> index;
1481 std::vector<MagnitudeType> wr(m2), wi(m2);
1484 std::vector<MagnitudeType> w(m2);
1487 SDM vr(m2,m2,
false);
1490 std::vector<int> iperm(m2);
1493 builtRecycleSpace_ =
true;
1503 SDM A_tmp( keff+m+blockSize_, keff+m );
1508 for (
int i=0; i<keff; ++i) { index[i] = i; }
1512 MVT::MvTransMv( one, *Ctmp, *Utmp, A11 );
1516 index.resize(m+blockSize_);
1517 for (i=0; i < m+blockSize_; i++) { index[i] = i; }
1519 MVT::MvTransMv( one, *Vp, *Utmp, A21 );
1522 for( i=keff; i<keff+m; i++)
1535 char balanc=
'P', jobvl=
'N', jobvr=
'V', sense=
'N';
1536 int ld =
A.numRows();
1538 int ldvl = ld, ldvr = ld;
1539 int info = 0,ilo = 0,ihi = 0;
1542 std::vector<ScalarType> beta(ld);
1543 std::vector<ScalarType> work(lwork);
1544 std::vector<MagnitudeType> rwork(lwork);
1545 std::vector<MagnitudeType> lscale(ld), rscale(ld);
1546 std::vector<MagnitudeType> rconde(ld), rcondv(ld);
1547 std::vector<int> iwork(ld+6);
1552 lapack.GGEVX(balanc, jobvl, jobvr, sense, ld,
A.values(), ld, B.
values(), ld, &wr[0], &wi[0],
1553 &beta[0], vl, ldvl, vr.
values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
1554 &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &rwork[0],
1555 &iwork[0], bwork, &info);
1560 for( i=0; i<ld; i++ )
1563 this->sort(w,ld,iperm);
1568 for( i=0; i<recycledBlocks_; i++ )
1569 for( j=0; j<ld; j++ )
1570 PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
1572 if(scalarTypeIsComplex==
false) {
1575 if (wi[iperm[ld-recycledBlocks_]] != 0.0) {
1577 for ( i=ld-recycledBlocks_; i<ld; i++ )
1578 if (wi[iperm[i]] != 0.0) countImag++;
1580 if (countImag % 2) xtraVec =
true;
1584 if (wi[iperm[ld-recycledBlocks_]] > 0.0) {
1585 for( j=0; j<ld; j++ )
1586 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
1589 for( j=0; j<ld; j++ )
1590 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
1598 return recycledBlocks_+1;
1600 return recycledBlocks_;
1604 template<
class ScalarType,
class MV,
class OP>
1606 bool xtraVec =
false;
1611 std::vector<MagnitudeType> wr(m), wi(m);
1617 std::vector<MagnitudeType> w(m);
1620 std::vector<int> iperm(m);
1626 builtRecycleSpace_ =
true;
1633 for(
int i=0; i<=blockSize_-1; i++) (*harmRitzMatrix)[blockSize_-1-i][harmRitzMatrix->numRows()-1-i] = 1;
1636 lapack.GESV(m, blockSize_, HHt.
values(), HHt.
stride(), &iperm[0], harmRitzMatrix->values(), harmRitzMatrix->stride(), &info);
1653 Htemp =
Teuchos::rcp(
new SDM(harmRitzMatrix -> numRows(), harmRitzMatrix -> numCols()));
1655 harmRitzMatrix -> assign(*Htemp);
1658 int harmColIndex, HHColIndex;
1660 for(
int i = 0; i<blockSize_; i++){
1661 harmColIndex = harmRitzMatrix -> numCols() - i -1;
1663 for(
int j=0; j<m; j++) (*Htemp)[HHColIndex][j] += (*harmRitzMatrix)[harmColIndex][j];
1665 harmRitzMatrix = Htemp;
1673 std::vector<ScalarType> work(1);
1674 std::vector<MagnitudeType> rwork(2*m);
1680 lapack.GEEV(
'N',
'V', m, harmRitzMatrix->values(), harmRitzMatrix->stride(), &wr[0], &wi[0],
1681 vl, ldvl, vr.
values(), vr.
stride(), &work[0], lwork, &rwork[0], &info);
1684 work.resize( lwork );
1686 lapack.GEEV(
'N',
'V', m, harmRitzMatrix->values(), harmRitzMatrix->stride(), &wr[0], &wi[0],
1687 vl, ldvl, vr.
values(), vr.
stride(), &work[0], lwork, &rwork[0], &info);
1692 for(
int i=0; i<m; ++i ) w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] );
1694 this->sort(w, m, iperm);
1699 for(
int i=0; i<recycledBlocks_; ++i )
1700 for(
int j=0; j<m; j++ )
1701 PP(j,i) = vr(j,iperm[i]);
1703 if(scalarTypeIsComplex==
false) {
1706 if (wi[iperm[recycledBlocks_-1]] != 0.0) {
1708 for (
int i=0; i<recycledBlocks_; ++i )
1709 if (wi[iperm[i]] != 0.0) countImag++;
1711 if (countImag % 2) xtraVec =
true;
1715 if (wi[iperm[recycledBlocks_-1]] > 0.0) {
1716 for(
int j=0; j<m; ++j )
1717 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
1720 for(
int j=0; j<m; ++j )
1721 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
1729 printer_->stream(
Debug) <<
"Recycled " << recycledBlocks_+1 <<
" vectors" << std::endl;
1730 return recycledBlocks_+1;
1733 printer_->stream(
Debug) <<
"Recycled " << recycledBlocks_ <<
" vectors" << std::endl;
1734 return recycledBlocks_;
1739 template<
class ScalarType,
class MV,
class OP>
1741 int l, r, j, i, flag;
1768 if (dlist[j] > dlist[j - 1]) j = j + 1;
1769 if (dlist[j - 1] > dK) {
1770 dlist[i - 1] = dlist[j - 1];
1771 iperm[i - 1] = iperm[j - 1];
1784 dlist[r] = dlist[0];
1785 iperm[r] = iperm[0];
1799 template<
class ScalarType,
class MV,
class OP>
1803 using Teuchos::rcp_const_cast;
1809 std::vector<int> index(numBlocks_+1);
1825 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1826 int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1830 std::vector<int> currIdx;
1832 if ( adaptiveBlockSize_ ) {
1833 blockSize_ = numCurrRHS;
1834 currIdx.resize( numCurrRHS );
1835 for (
int i=0; i<numCurrRHS; ++i)
1836 currIdx[i] = startPtr+i;
1839 currIdx.resize( blockSize_ );
1840 for (
int i=0; i<numCurrRHS; ++i)
1841 currIdx[i] = startPtr+i;
1842 for (
int i=numCurrRHS; i<blockSize_; ++i)
1847 problem_->setLSIndex( currIdx );
1853 loaDetected_ =
false;
1856 bool isConverged =
true;
1859 initializeStateStorage();
1864 while (numRHS2Solve > 0){
1874 printer_->stream(
Debug) <<
" Now solving RHS index " << currIdx[0] <<
" using recycled subspace of dimension " << keff << std::endl << std::endl;
1878 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1879 RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1880 RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1881 problem_->apply( *Utmp, *Ctmp );
1883 RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1888 int rank = ortho_->normalize(*Ctmp,
rcp(&Ftmp,
false));
1901 work_.resize(lwork);
1906 MVT::MvTimesMatAddMv( one, *Utmp, Ftmp, zero, *U1tmp );
1911 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1912 Ctmp = MVT::CloneViewNonConst( *C_, index );
1913 Utmp = MVT::CloneView( *U_, index );
1916 SDM Ctr(keff,blockSize_);
1917 problem_->computeCurrPrecResVec( &*R_ );
1918 MVT::MvTransMv( one, *Ctmp, *R_, Ctr );
1921 RCP<MV> update = MVT::Clone( *problem_->getCurrLHSVec(), blockSize_ );
1922 MVT::MvInit( *update, 0.0 );
1923 MVT::MvTimesMatAddMv( one, *Utmp, Ctr, one, *update );
1924 problem_->updateSolution( update,
true );
1927 MVT::MvTimesMatAddMv( -one, *Ctmp, Ctr, one, *R_ );
1938 V_ = MVT::Clone( *rhsMV, (numBlocks_+1)*blockSize_ );
1942 if (MVT::GetNumberVecs(*V_) < (numBlocks_+1)*blockSize_ ) {
1944 V_ = MVT::Clone( *tmp, (numBlocks_+1)*blockSize_ );
1949 printer_->stream(
Debug) <<
" No recycled subspace available for RHS index " << std::endl << std::endl;
1954 primeList.
set(
"Num Blocks",numBlocks_-1);
1955 primeList.
set(
"Block Size",blockSize_);
1956 primeList.
set(
"Recycled Blocks",0);
1957 primeList.
set(
"Keep Hessenberg",
true);
1958 primeList.
set(
"Initialize Hessenberg",
true);
1960 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1961 if (blockSize_*static_cast<ptrdiff_t>(numBlocks_) > dim) {
1962 ptrdiff_t tmpNumBlocks = 0;
1963 if (blockSize_ == 1)
1964 tmpNumBlocks = dim / blockSize_;
1966 tmpNumBlocks = ( dim - blockSize_) / blockSize_;
1967 printer_->stream(
Warnings) <<
"Belos::BlockGCRODRSolMgr::solve(): Warning! Requested Krylov subspace dimension is larger than operator dimension!"
1968 << std::endl <<
"The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << tmpNumBlocks << std::endl;
1969 primeList.
set(
"Num Blocks",Teuchos::as<int>(tmpNumBlocks));
1973 primeList.
set(
"Num Blocks",numBlocks_-1);
1980 block_gmres_iter->setSize( blockSize_, numBlocks_-1 );
1984 if (currIdx[blockSize_-1] == -1) {
1985 V_0 = MVT::Clone( *(problem_->getInitPrecResVec()), blockSize_ );
1986 problem_->computeCurrPrecResVec( &*V_0 );
1989 V_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
1996 int rank = ortho_->normalize( *V_0, z_0 );
2005 block_gmres_iter->initializeGmres(newstate);
2007 bool primeConverged =
false;
2010 printer_->stream(
Debug) <<
" Preparing to Iterate!!!!" << std::endl << std::endl;
2011 block_gmres_iter->iterate();
2016 if ( convTest_->getStatus() ==
Passed ) {
2017 printer_->stream(
Debug) <<
"We converged during the prime the pump stage" << std::endl << std::endl;
2018 primeConverged = !(expConvTest_->getLOADetected());
2019 if ( expConvTest_->getLOADetected() ) {
2021 loaDetected_ =
true;
2022 printer_->stream(
Warnings) <<
"Belos::BlockGCRODRSolMgr::solve(): Warning! Solver has experienced a loss of accuracy!" << std::endl;
2028 else if( maxIterTest_->getStatus() ==
Passed ) {
2030 primeConverged =
false;
2036 printer_->stream(
Debug) <<
" We did not converge on priming cycle of Block GMRES. Therefore we recycle and restart. " << std::endl << std::endl;
2041 if (blockSize_ != 1) {
2042 printer_->stream(
Errors) <<
"Error! Caught std::exception in BlockGmresIter::iterate() at iteration "
2043 << block_gmres_iter->getNumIters() << std::endl
2044 << e.what() << std::endl;
2045 if (convTest_->getStatus() !=
Passed)
2046 primeConverged =
false;
2050 block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
2052 sTest_->checkStatus( &*block_gmres_iter );
2053 if (convTest_->getStatus() !=
Passed)
2054 isConverged =
false;
2059 achievedTol_ = MT::one();
2061 MVT::MvInit( *X, SCT::zero() );
2062 printer_->stream(
Warnings) <<
"Belos::BlockGCRODRSolMgr::solve(): Warning! NaN has been detected!"
2066 catch (
const std::exception &e) {
2067 printer_->stream(
Errors) <<
"Error! Caught std::exception in BlockGmresIter::iterate() at iteration "
2068 << block_gmres_iter->getNumIters() << std::endl
2069 << e.what() << std::endl;
2077 RCP<MV> update = block_gmres_iter->getCurrentUpdate();
2078 problem_->updateSolution( update,
true );
2082 problem_->computeCurrPrecResVec( &*R_ );
2085 newstate = block_gmres_iter->getState();
2088 H_->assign(*(newstate.
H));
2097 V_ = rcp_const_cast<MV>(newstate.
V);
2100 buildRecycleSpaceKryl(keff, block_gmres_iter);
2101 printer_->stream(
Debug) <<
"Generated recycled subspace using RHS index " << currIdx[0] <<
" of dimension " << keff << std::endl << std::endl;
2105 if (primeConverged) {
2126 problem_->setCurrLS();
2129 startPtr += numCurrRHS;
2130 numRHS2Solve -= numCurrRHS;
2131 if ( numRHS2Solve > 0 ) {
2132 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
2133 if ( adaptiveBlockSize_ ) {
2134 blockSize_ = numCurrRHS;
2135 currIdx.resize( numCurrRHS );
2136 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2139 currIdx.resize( blockSize_ );
2140 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2141 for (
int i=numCurrRHS; i<blockSize_; ++i) currIdx[i] = -1;
2144 problem_->setLSIndex( currIdx );
2147 currIdx.resize( numRHS2Solve );
2158 blockgcrodrList.
set(
"Num Blocks",numBlocks_);
2159 blockgcrodrList.
set(
"Block Size",blockSize_);
2160 blockgcrodrList.
set(
"Recycled Blocks",keff);
2166 index.resize( blockSize_ );
2167 for(
int ii = 0; ii < blockSize_; ii++) index[ii] = ii;
2171 MVT::Assign(*R_,*V0);
2174 for(
int i=0; i < keff; i++){ index[i] = i;};
2176 H_ =
rcp(
new SDM(
Teuchos::View, *G_, (numBlocks_-1)*blockSize_ + blockSize_, (numBlocks_-1)*blockSize_, keff ,keff ));
2180 newstate.
U = MVT::CloneViewNonConst(*U_, index);
2181 newstate.
C = MVT::CloneViewNonConst(*C_, index);
2183 newstate.
curDim = blockSize_;
2184 block_gcrodr_iter -> initialize(newstate);
2186 int numRestarts = 0;
2190 block_gcrodr_iter -> iterate();
2195 if( convTest_->getStatus() ==
Passed ) {
2203 else if(maxIterTest_->getStatus() ==
Passed ){
2205 isConverged =
false;
2212 else if (block_gcrodr_iter->getCurSubspaceDim() == block_gcrodr_iter->getMaxSubspaceDim()){
2218 problem_->updateSolution(update,
true);
2219 buildRecycleSpaceAugKryl(block_gcrodr_iter);
2221 printer_->stream(
Debug) <<
" Generated new recycled subspace using RHS index " << currIdx[0] <<
" of dimension " << keff << std::endl << std::endl;
2223 if(numRestarts >= maxRestarts_) {
2224 isConverged =
false;
2230 printer_ -> stream(
Debug) <<
" Performing restart number " << numRestarts <<
" of " << maxRestarts_ << std::endl << std::endl;
2233 problem_->computeCurrPrecResVec( &*R_ );
2234 index.resize( blockSize_ );
2235 for (
int ii=0; ii<blockSize_; ++ii) index[ii] = ii;
2237 MVT::SetBlock(*R_,index,*V0);
2241 index.resize( numBlocks_*blockSize_ );
2242 for (
int ii=0; ii<(numBlocks_*blockSize_); ++ii) index[ii] = ii;
2243 restartState.
V = MVT::CloneViewNonConst( *V_, index );
2244 index.resize( keff );
2245 for (
int ii=0; ii<keff; ++ii) index[ii] = ii;
2246 restartState.
U = MVT::CloneViewNonConst( *U_, index );
2247 restartState.
C = MVT::CloneViewNonConst( *C_, index );
2249 H_ =
rcp(
new SDM(
Teuchos::View, *G_, numBlocks_*blockSize_, (numBlocks_-1)*blockSize_, keff ,keff ));
2250 restartState.
B = B_;
2251 restartState.
H = H_;
2252 restartState.
curDim = blockSize_;
2253 block_gcrodr_iter->initialize(restartState);
2262 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Belos::BlockGCRODRSolMgr::solve(): Invalid return from BlockGCRODRIter::iterate().");
2268 block_gcrodr_iter->updateLSQR( block_gcrodr_iter->getCurSubspaceDim() );
2270 sTest_->checkStatus( &*block_gcrodr_iter );
2271 if (convTest_->getStatus() !=
Passed) isConverged =
false;
2274 catch(
const std::exception &e){
2275 printer_->stream(
Errors) <<
"Error! Caught exception in BlockGCRODRIter::iterate() at iteration "
2276 << block_gcrodr_iter->getNumIters() << std::endl
2277 << e.what() << std::endl;
2285 problem_->updateSolution( update,
true );
2304 problem_->setCurrLS();
2307 startPtr += numCurrRHS;
2308 numRHS2Solve -= numCurrRHS;
2309 if ( numRHS2Solve > 0 ) {
2310 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
2311 if ( adaptiveBlockSize_ ) {
2312 blockSize_ = numCurrRHS;
2313 currIdx.resize( numCurrRHS );
2314 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2317 currIdx.resize( blockSize_ );
2318 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2319 for (
int i=numCurrRHS; i<blockSize_; ++i) currIdx[i] = -1;
2322 problem_->setLSIndex( currIdx );
2325 currIdx.resize( numRHS2Solve );
2329 if (!builtRecycleSpace_) {
2330 buildRecycleSpaceAugKryl(block_gcrodr_iter);
2331 printer_->stream(
Debug) <<
" Generated new recycled subspace using RHS index " << currIdx[0] <<
" of dimension " << keff << std::endl << std::endl;
2339 #ifdef BELOS_TEUCHOS_TIME_MONITOR
2345 numIters_ = maxIterTest_->getNumIters();
2348 const std::vector<MagnitudeType>* pTestValues = NULL;
2349 pTestValues = impConvTest_->getTestValue();
2351 "Belos::BlockGCRODRSolMgr::solve(): The implicit convergence test's "
2352 "getTestValue() method returned NULL. Please report this bug to the "
2353 "Belos developers.");
2355 "Belos::BlockGCRODRSolMgr::solve(): The implicit convergence test's "
2356 "getTestValue() method returned a vector of length zero. Please report "
2357 "this bug to the Belos developers.");
2361 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
ScalarType * values() const
Collection of types and exceptions used within the Belos solvers.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get current linear problem being solved for in this object.
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_
int getNumIters() const
Get the iteration count for the most recent call to solve().
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)
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
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(ParameterList &l, const std::string &name)
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)
Teuchos::ScalarTraits< MagnitudeType > SMT
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
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.
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.
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.
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.
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_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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()
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)
bool is_null(const RCP< T > &p)
std::string recycleMethod_
MagnitudeType achievedTol_
MagnitudeType achievedTol() const
Get the residual for the most recent call to solve().
Belos concrete class for performing the block GCRO-DR (block GMRES with recycling) iteration...
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
std::vector< ScalarType > work_
BlockGCRODRSolMgr()
Default constructor.
A linear system to solve, and its associated information.
std::string description() const
A description of the Block GCRODR solver manager.
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)
Teuchos::RCP< OutputManager< ScalarType > > printer_
ReturnType
Whether the Belos solve converged for all linear systems.
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< 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)
static magnitudeType magnitude(T a)
OrdinalType numCols() const
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_
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Get a parameter list containing the valid parameters for this object.
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_
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
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.
OrdinalType stride() const
ortho_factory_type orthoFactory_
Factory for creating MatOrthoManager subclass instances.
OrdinalType numRows() const
Belos concrete class for performing the block, flexible GMRES iteration.
bool isLOADetected() const
Whether a loss of accuracy was detected during the most recent solve.
OutputType
Style of output used to display status test information.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.