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.");
322 void sort (std::vector<MagnitudeType>& dlist,
int n, std::vector<int>& iperm);
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>
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);
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")) {
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){
1157 if ( (G_->numRows() != augSpaImgDim) || (G_->numCols() != augSpaDim) )
1159 G_->reshape( augSpaImgDim, augSpaDim );
1161 G_->putScalar(zero);
1165 if (H_ == Teuchos::null){
1170 if (F_ == Teuchos::null){
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){
1185 if ( (PP_->numRows() != augSpaImgDim) || (PP_->numCols() != recycledBlocks_+1) ){
1186 PP_->reshape( augSpaImgDim, recycledBlocks_+1 );
1191 if (HP_ == Teuchos::null)
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>
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 );
1265 work_.resize(lwork);
1272 for(
int ii=0;ii<keff;ii++) {
for(
int jj=ii;jj<keff;jj++) Rtmp(ii,jj) = HPtmp(ii,jj); }
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 );
1294 work_.resize(lwork);
1298 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
1304 template<
class ScalarType,
class MV,
class OP>
1309 std::vector<MagnitudeType> d(keff);
1310 std::vector<ScalarType> dscalar(keff);
1311 std::vector<int> index(numBlocks_+1);
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];
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; }
1384 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
1395 int info = 0, lwork = -1;
1396 tau_.resize(keff_new);
1401 work_.resize(lwork);
1408 for(
int i=0;i<keff_new;i++) {
for(
int j=i;j<keff_new;j++) Ftmp(i,j) = HPtmp(i,j); }
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 );
1450 work_.resize(lwork);
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_ );
1473 template<
class ScalarType,
class MV,
class OP>
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++)
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;
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);
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>
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);
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);
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>
1737 int l, r, j, i, flag;
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.
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_
ScalarType * values() const
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.
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.
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_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. ...
OrdinalType numRows() const
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_
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_
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().
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)
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.
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.
OrdinalType stride() const
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 ...
OrdinalType numCols() const
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.
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.