42 #ifndef BELOS_GCRODR_SOLMGR_HPP 43 #define BELOS_GCRODR_SOLMGR_HPP 65 #ifdef BELOS_TEUCHOS_TIME_MONITOR 67 #endif // BELOS_TEUCHOS_TIME_MONITOR 68 #if defined(HAVE_TEUCHOSCORE_CXX11) 69 # include <type_traits> 70 #endif // defined(HAVE_TEUCHOSCORE_CXX11) 154 template<
class ScalarType,
class MV,
class OP,
155 const bool lapackSupportsScalarType =
179 template<
class ScalarType,
class MV,
class OP>
183 #if defined(HAVE_TEUCHOSCORE_CXX11) 184 # if defined(HAVE_TEUCHOS_COMPLEX) 185 static_assert (std::is_same<ScalarType, std::complex<float> >::value ||
186 std::is_same<ScalarType, std::complex<double> >::value ||
187 std::is_same<ScalarType, float>::value ||
188 std::is_same<ScalarType, double>::value,
189 "Belos::GCRODRSolMgr: ScalarType must be one of the four " 190 "types (S,D,C,Z) supported by LAPACK.");
192 static_assert (std::is_same<ScalarType, float>::value ||
193 std::is_same<ScalarType, double>::value,
194 "Belos::GCRODRSolMgr: ScalarType must be float or double. " 195 "Complex arithmetic support is currently disabled. To " 196 "enable it, set Teuchos_ENABLE_COMPLEX=ON.");
197 # endif // defined(HAVE_TEUCHOS_COMPLEX) 198 #endif // defined(HAVE_TEUCHOSCORE_CXX11) 308 return Teuchos::tuple(timerSolve_);
352 bool set = problem_->setProblem();
354 throw "Could not set problem.";
398 std::string description()
const override;
408 void initializeStateStorage();
417 int getHarmonicVecs1(
int m,
426 int getHarmonicVecs2(
int keff,
int m,
432 void sort(std::vector<MagnitudeType>& dlist,
int n, std::vector<int>& iperm);
460 static constexpr
double orthoKappa_default_ = 0.0;
461 static constexpr
int maxRestarts_default_ = 100;
462 static constexpr
int maxIters_default_ = 1000;
463 static constexpr
int numBlocks_default_ = 50;
464 static constexpr
int blockSize_default_ = 1;
465 static constexpr
int recycledBlocks_default_ = 5;
468 static constexpr
int outputFreq_default_ = -1;
469 static constexpr
const char * impResScale_default_ =
"Norm of Preconditioned Initial Residual";
470 static constexpr
const char * expResScale_default_ =
"Norm of Initial Residual";
471 static constexpr
const char * label_default_ =
"Belos";
472 static constexpr
const char * orthoType_default_ =
"DGKS";
473 static constexpr std::ostream * outputStream_default_ = &std::cout;
528 template<
class ScalarType,
class MV,
class OP>
538 template<
class ScalarType,
class MV,
class OP>
550 problem == Teuchos::null, std::invalid_argument,
551 "Belos::GCRODRSolMgr constructor: The solver manager's " 552 "constructor needs the linear problem argument 'problem' " 567 template<
class ScalarType,
class MV,
class OP>
569 outputStream_ =
Teuchos::rcp(outputStream_default_,
false);
571 orthoKappa_ = orthoKappa_default_;
572 maxRestarts_ = maxRestarts_default_;
573 maxIters_ = maxIters_default_;
574 numBlocks_ = numBlocks_default_;
575 recycledBlocks_ = recycledBlocks_default_;
576 verbosity_ = verbosity_default_;
577 outputStyle_ = outputStyle_default_;
578 outputFreq_ = outputFreq_default_;
579 orthoType_ = orthoType_default_;
580 impResScale_ = impResScale_default_;
581 expResScale_ = expResScale_default_;
582 label_ = label_default_;
584 builtRecycleSpace_ =
false;
600 template<
class ScalarType,
class MV,
class OP>
605 using Teuchos::isParameterType;
606 using Teuchos::getParameter;
609 using Teuchos::parameterList;
612 using Teuchos::rcp_dynamic_cast;
613 using Teuchos::rcpFromRef;
620 RCP<const ParameterList> defaultParams = getValidParameters();
638 if (params_.is_null()) {
639 params_ = parameterList (*defaultParams);
647 if (params_ != params) {
653 params_ = parameterList (*params);
688 params_->validateParametersAndSetDefaults (*defaultParams);
693 maxRestarts_ = params->
get(
"Maximum Restarts", maxRestarts_default_);
696 params_->set (
"Maximum Restarts", maxRestarts_);
701 maxIters_ = params->
get (
"Maximum Iterations", maxIters_default_);
704 params_->set (
"Maximum Iterations", maxIters_);
705 if (! maxIterTest_.is_null())
706 maxIterTest_->setMaxIters (maxIters_);
711 numBlocks_ = params->
get (
"Num Blocks", numBlocks_default_);
713 "Belos::GCRODRSolMgr: The \"Num Blocks\" parameter must " 714 "be strictly positive, but you specified a value of " 715 << numBlocks_ <<
".");
717 params_->set (
"Num Blocks", numBlocks_);
722 recycledBlocks_ = params->
get (
"Num Recycled Blocks",
723 recycledBlocks_default_);
725 "Belos::GCRODRSolMgr: The \"Num Recycled Blocks\" " 726 "parameter must be strictly positive, but you specified " 727 "a value of " << recycledBlocks_ <<
".");
729 "Belos::GCRODRSolMgr: The \"Num Recycled Blocks\" " 730 "parameter must be less than the \"Num Blocks\" " 731 "parameter, but you specified \"Num Recycled Blocks\" " 732 "= " << recycledBlocks_ <<
" and \"Num Blocks\" = " 733 << numBlocks_ <<
".");
735 params_->set(
"Num Recycled Blocks", recycledBlocks_);
742 std::string tempLabel = params->
get (
"Timer Label", label_default_);
745 if (tempLabel != label_) {
747 params_->set (
"Timer Label", label_);
748 std::string solveLabel = label_ +
": GCRODRSolMgr total solve time";
749 #ifdef BELOS_TEUCHOS_TIME_MONITOR 752 if (ortho_ != Teuchos::null) {
753 ortho_->setLabel( label_ );
760 if (isParameterType<int> (*params,
"Verbosity")) {
761 verbosity_ = params->
get (
"Verbosity", verbosity_default_);
763 verbosity_ = (int) getParameter<Belos::MsgType> (*params,
"Verbosity");
766 params_->set (
"Verbosity", verbosity_);
769 if (! printer_.is_null())
770 printer_->setVerbosity (verbosity_);
775 if (isParameterType<int> (*params,
"Output Style")) {
776 outputStyle_ = params->
get (
"Output Style", outputStyle_default_);
778 outputStyle_ = (int) getParameter<OutputType> (*params,
"Output Style");
782 params_->set (
"Output Style", outputStyle_);
802 outputStream_ = getParameter<RCP<std::ostream> > (*params,
"Output Stream");
803 }
catch (InvalidParameter&) {
804 outputStream_ = rcpFromRef (std::cout);
811 if (outputStream_.is_null()) {
815 params_->set (
"Output Stream", outputStream_);
818 if (! printer_.is_null()) {
819 printer_->setOStream (outputStream_);
826 outputFreq_ = params->
get (
"Output Frequency", outputFreq_default_);
830 params_->set(
"Output Frequency", outputFreq_);
831 if (! outputTest_.is_null())
832 outputTest_->setOutputFrequency (outputFreq_);
839 if (printer_.is_null()) {
850 bool changedOrthoType =
false;
852 const std::string& tempOrthoType =
853 params->
get (
"Orthogonalization", orthoType_default_);
856 std::ostringstream os;
857 os <<
"Belos::GCRODRSolMgr: Invalid orthogonalization name \"" 858 << tempOrthoType <<
"\". The following are valid options " 859 <<
"for the \"Orthogonalization\" name parameter: ";
861 throw std::invalid_argument (os.str());
863 if (tempOrthoType != orthoType_) {
864 changedOrthoType =
true;
865 orthoType_ = tempOrthoType;
867 params_->set (
"Orthogonalization", orthoType_);
883 RCP<ParameterList> orthoParams;
886 using Teuchos::sublist;
888 const std::string paramName (
"Orthogonalization Parameters");
891 orthoParams = sublist (params_, paramName,
true);
892 }
catch (InvalidParameter&) {
899 orthoParams = sublist (params_, paramName,
true);
903 "Failed to get orthogonalization parameters. " 904 "Please report this bug to the Belos developers.");
909 if (ortho_.is_null() || changedOrthoType) {
915 label_, orthoParams);
924 RCP<PLA> pla = rcp_dynamic_cast<PLA> (ortho_);
930 label_, orthoParams);
932 pla->setParameterList (orthoParams);
944 if (params->
isParameter (
"Orthogonalization Constant")) {
947 orthoKappa = params->
get (
"Orthogonalization Constant", orthoKappa);
950 orthoKappa = params->
get (
"Orthogonalization Constant", orthoKappa_default_);
953 if (orthoKappa > 0) {
954 orthoKappa_ = orthoKappa;
956 params_->set(
"Orthogonalization Constant", orthoKappa_);
958 if (orthoType_ ==
"DGKS" && ! ortho_.is_null()) {
965 rcp_dynamic_cast<ortho_man_type>(ortho_)->setDepTol (orthoKappa_);
975 if (params->
isParameter(
"Convergence Tolerance")) {
977 convTol_ = params->
get (
"Convergence Tolerance",
985 params_->set (
"Convergence Tolerance", convTol_);
986 if (! impConvTest_.is_null())
987 impConvTest_->setTolerance (convTol_);
988 if (! expConvTest_.is_null())
989 expConvTest_->setTolerance (convTol_);
993 if (params->
isParameter (
"Implicit Residual Scaling")) {
994 std::string tempImpResScale =
995 getParameter<std::string> (*params,
"Implicit Residual Scaling");
998 if (impResScale_ != tempImpResScale) {
1000 impResScale_ = tempImpResScale;
1003 params_->set(
"Implicit Residual Scaling", impResScale_);
1013 if (! impConvTest_.is_null()) {
1019 impConvTest_ = null;
1026 if (params->
isParameter(
"Explicit Residual Scaling")) {
1027 std::string tempExpResScale =
1028 getParameter<std::string> (*params,
"Explicit Residual Scaling");
1031 if (expResScale_ != tempExpResScale) {
1033 expResScale_ = tempExpResScale;
1036 params_->set(
"Explicit Residual Scaling", expResScale_);
1039 if (! expConvTest_.is_null()) {
1045 expConvTest_ = null;
1056 if (maxIterTest_.is_null())
1061 if (impConvTest_.is_null()) {
1062 impConvTest_ =
rcp (
new StatusTestResNorm_t (convTol_));
1068 if (expConvTest_.is_null()) {
1069 expConvTest_ =
rcp (
new StatusTestResNorm_t (convTol_));
1070 expConvTest_->defineResForm (StatusTestResNorm_t::Explicit,
Belos::TwoNorm);
1076 if (convTest_.is_null()) {
1077 convTest_ =
rcp (
new StatusTestCombo_t (StatusTestCombo_t::SEQ,
1085 sTest_ =
rcp (
new StatusTestCombo_t (StatusTestCombo_t::OR,
1091 outputTest_ = stoFactory.
create (printer_, sTest_, outputFreq_,
1095 std::string solverDesc =
" GCRODR ";
1096 outputTest_->setSolverDesc( solverDesc );
1099 if (timerSolve_.is_null()) {
1100 std::string solveLabel = label_ +
": GCRODRSolMgr total solve time";
1101 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1111 template<
class ScalarType,
class MV,
class OP>
1116 using Teuchos::parameterList;
1119 static RCP<const ParameterList> validPL;
1121 RCP<ParameterList> pl = parameterList ();
1125 "The relative residual tolerance that needs to be achieved by the\n" 1126 "iterative solver in order for the linear system to be declared converged.");
1127 pl->set(
"Maximum Restarts", static_cast<int>(maxRestarts_default_),
1128 "The maximum number of cycles allowed for each\n" 1129 "set of RHS solved.");
1130 pl->set(
"Maximum Iterations", static_cast<int>(maxIters_default_),
1131 "The maximum number of iterations allowed for each\n" 1132 "set of RHS solved.");
1136 pl->set(
"Block Size", static_cast<int>(blockSize_default_),
1137 "Block Size Parameter -- currently must be 1 for GCRODR");
1138 pl->set(
"Num Blocks", static_cast<int>(numBlocks_default_),
1139 "The maximum number of vectors allowed in the Krylov subspace\n" 1140 "for each set of RHS solved.");
1141 pl->set(
"Num Recycled Blocks", static_cast<int>(recycledBlocks_default_),
1142 "The maximum number of vectors in the recycled subspace." );
1143 pl->set(
"Verbosity", static_cast<int>(verbosity_default_),
1144 "What type(s) of solver information should be outputted\n" 1145 "to the output stream.");
1146 pl->set(
"Output Style", static_cast<int>(outputStyle_default_),
1147 "What style is used for the solver information outputted\n" 1148 "to the output stream.");
1149 pl->set(
"Output Frequency", static_cast<int>(outputFreq_default_),
1150 "How often convergence information should be outputted\n" 1151 "to the output stream.");
1152 pl->set(
"Output Stream",
Teuchos::rcp(outputStream_default_,
false),
1153 "A reference-counted pointer to the output stream where all\n" 1154 "solver output is sent.");
1155 pl->set(
"Implicit Residual Scaling", static_cast<const char *>(impResScale_default_),
1156 "The type of scaling used in the implicit residual convergence test.");
1157 pl->set(
"Explicit Residual Scaling", static_cast<const char *>(expResScale_default_),
1158 "The type of scaling used in the explicit residual convergence test.");
1159 pl->set(
"Timer Label", static_cast<const char *>(label_default_),
1160 "The string to use as a prefix for the timer labels.");
1163 pl->set(
"Orthogonalization", static_cast<const char *>(orthoType_default_),
1164 "The type of orthogonalization to use. Valid options: " +
1166 RCP<const ParameterList> orthoParams =
1168 pl->
set (
"Orthogonalization Parameters", *orthoParams,
1169 "Parameters specific to the type of orthogonalization used.");
1171 pl->
set(
"Orthogonalization Constant",static_cast<MagnitudeType>(orthoKappa_default_),
1172 "When using DGKS orthogonalization: the \"depTol\" constant, used " 1173 "to determine whether another step of classical Gram-Schmidt is " 1174 "necessary. Otherwise ignored.");
1181 template<
class ScalarType,
class MV,
class OP>
1188 if (rhsMV == Teuchos::null) {
1196 "Belos::GCRODRSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
1199 if (U_ == Teuchos::null) {
1200 U_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1204 if (MVT::GetNumberVecs(*U_) < recycledBlocks_+1) {
1206 U_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1211 if (C_ == Teuchos::null) {
1212 C_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1216 if (MVT::GetNumberVecs(*C_) < recycledBlocks_+1) {
1218 C_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1223 if (V_ == Teuchos::null) {
1224 V_ = MVT::Clone( *rhsMV, numBlocks_+1 );
1228 if (MVT::GetNumberVecs(*V_) < numBlocks_+1) {
1230 V_ = MVT::Clone( *tmp, numBlocks_+1 );
1235 if (U1_ == Teuchos::null) {
1236 U1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1240 if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1242 U1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1247 if (C1_ == Teuchos::null) {
1248 C1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1252 if (MVT::GetNumberVecs(*C1_) < recycledBlocks_+1) {
1254 C1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1259 if (r_ == Teuchos::null)
1260 r_ = MVT::Clone( *rhsMV, 1 );
1263 tau_.resize(recycledBlocks_+1);
1266 work_.resize(recycledBlocks_+1);
1269 ipiv_.resize(recycledBlocks_+1);
1272 if (H2_ == Teuchos::null)
1275 if ( (H2_->numRows() != numBlocks_+recycledBlocks_+2) || (H2_->numCols() != numBlocks_+recycledBlocks_+1) )
1276 H2_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
1278 H2_->putScalar(zero);
1281 if (R_ == Teuchos::null)
1284 if ( (R_->numRows() != recycledBlocks_+1) || (R_->numCols() != recycledBlocks_+1) )
1285 R_->reshape( recycledBlocks_+1, recycledBlocks_+1 );
1287 R_->putScalar(zero);
1290 if (PP_ == Teuchos::null)
1293 if ( (PP_->numRows() != numBlocks_+recycledBlocks_+2) || (PP_->numCols() != recycledBlocks_+1) )
1294 PP_->reshape( numBlocks_+recycledBlocks_+2, recycledBlocks_+1 );
1298 if (HP_ == Teuchos::null)
1301 if ( (HP_->numRows() != numBlocks_+recycledBlocks_+2) || (HP_->numCols() != numBlocks_+recycledBlocks_+1) )
1302 HP_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
1310 template<
class ScalarType,
class MV,
class OP>
1318 if (!isSet_) { setParameters( params_ ); }
1322 std::vector<int> index(numBlocks_+1);
1329 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1330 std::vector<int> currIdx(1);
1334 problem_->setLSIndex( currIdx );
1337 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1338 if (static_cast<ptrdiff_t>(numBlocks_) > dim) {
1339 numBlocks_ = Teuchos::as<int>(dim);
1341 "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
1342 " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
1343 params_->set(
"Num Blocks", numBlocks_);
1347 bool isConverged =
true;
1350 initializeStateStorage();
1356 plist.
set(
"Num Blocks",numBlocks_);
1357 plist.
set(
"Recycled Blocks",recycledBlocks_);
1362 RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_iter;
1365 int prime_iterations = 0;
1369 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1373 while ( numRHS2Solve > 0 ) {
1376 builtRecycleSpace_ =
false;
1379 outputTest_->reset();
1387 "Belos::GCRODRSolMgr::solve(): Requested size of recycled subspace is not consistent with the current recycle subspace.");
1389 printer_->stream(
Debug) <<
" Now solving RHS index " << currIdx[0] <<
" using recycled subspace of dimension " << keff << std::endl << std::endl;
1392 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1393 RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1394 RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1395 problem_->apply( *Utmp, *Ctmp );
1397 RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1402 int rank = ortho_->normalize(*Ctmp,
rcp(&Rtmp,
false));
1415 work_.resize(lwork);
1420 MVT::MvTimesMatAddMv( one, *Utmp, Rtmp, zero, *U1tmp );
1425 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1426 Ctmp = MVT::CloneViewNonConst( *C_, index );
1427 Utmp = MVT::CloneView( *U_, index );
1431 problem_->computeCurrPrecResVec( &*r_ );
1432 MVT::MvTransMv( one, *Ctmp, *r_, Ctr );
1435 RCP<MV> update = MVT::Clone( *problem_->getCurrLHSVec(), 1 );
1436 MVT::MvInit( *update, 0.0 );
1437 MVT::MvTimesMatAddMv( one, *Utmp, Ctr, one, *update );
1438 problem_->updateSolution( update,
true );
1441 MVT::MvTimesMatAddMv( -one, *Ctmp, Ctr, one, *r_ );
1444 prime_iterations = 0;
1450 printer_->stream(
Debug) <<
" No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
1455 primeList.
set(
"Num Blocks",numBlocks_);
1456 primeList.
set(
"Recycled Blocks",0);
1459 RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_prime_iter;
1463 problem_->computeCurrPrecResVec( &*r_ );
1464 index.resize( 1 ); index[0] = 0;
1465 RCP<MV> v0 = MVT::CloneViewNonConst( *V_, index );
1466 MVT::SetBlock(*r_,index,*v0);
1470 index.resize( numBlocks_+1 );
1471 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1472 newstate.
V = MVT::CloneViewNonConst( *V_, index );
1473 newstate.
U = Teuchos::null;
1474 newstate.
C = Teuchos::null;
1476 newstate.
B = Teuchos::null;
1478 gcrodr_prime_iter->initialize(newstate);
1481 bool primeConverged =
false;
1483 gcrodr_prime_iter->iterate();
1486 if ( convTest_->getStatus() ==
Passed ) {
1488 primeConverged =
true;
1493 gcrodr_prime_iter->updateLSQR( gcrodr_prime_iter->getCurSubspaceDim() );
1496 sTest_->checkStatus( &*gcrodr_prime_iter );
1497 if (convTest_->getStatus() ==
Passed)
1498 primeConverged =
true;
1500 catch (
const std::exception &e) {
1501 printer_->stream(
Errors) <<
"Error! Caught exception in GCRODRIter::iterate() at iteration " 1502 << gcrodr_prime_iter->getNumIters() << std::endl
1503 << e.what() << std::endl;
1507 prime_iterations = gcrodr_prime_iter->getNumIters();
1510 RCP<MV> update = gcrodr_prime_iter->getCurrentUpdate();
1511 problem_->updateSolution( update,
true );
1514 newstate = gcrodr_prime_iter->getState();
1522 if (recycledBlocks_ < p+1) {
1526 keff = getHarmonicVecs1( p, *newstate.
H, *PPtmp );
1531 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1532 RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1533 RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1534 RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1536 for (
int ii=0; ii < p; ++ii) { index[ii] = ii; }
1537 RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1541 MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
1558 HPtmp.
stride (), &tau_[0], &work_[0], lwork, &info);
1561 " LAPACK's _GEQRF failed to compute a workspace size.");
1570 work_.resize (lwork);
1572 HPtmp.
stride (), &tau_[0], &work_[0], lwork, &info);
1575 " LAPACK's _GEQRF failed to compute a QR factorization.");
1580 for (
int ii = 0; ii < keff; ++ii) {
1581 for (
int jj = ii; jj < keff; ++jj) {
1582 Rtmp(ii,jj) = HPtmp(ii,jj);
1594 "LAPACK's _UNGQR failed to construct the Q factor.");
1599 index.resize (p + 1);
1600 for (
int ii = 0; ii < (p+1); ++ii) {
1603 Vtmp = MVT::CloneView( *V_, index );
1604 MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
1615 "LAPACK's _GETRF failed to compute an LU factorization.");
1625 work_.resize(lwork);
1629 "LAPACK's _GETRI failed to invert triangular matrix.");
1632 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
1634 printer_->stream(
Debug)
1635 <<
" Generated recycled subspace using RHS index " << currIdx[0]
1636 <<
" of dimension " << keff << std::endl << std::endl;
1641 if (primeConverged) {
1643 problem_->setCurrLS();
1647 if (numRHS2Solve > 0) {
1649 problem_->setLSIndex (currIdx);
1652 currIdx.resize (numRHS2Solve);
1662 gcrodr_iter->setSize( keff, numBlocks_ );
1665 gcrodr_iter->resetNumIters(prime_iterations);
1668 outputTest_->resetNumCalls();
1671 problem_->computeCurrPrecResVec( &*r_ );
1672 index.resize( 1 ); index[0] = 0;
1673 RCP<MV> v0 = MVT::CloneViewNonConst( *V_, index );
1674 MVT::SetBlock(*r_,index,*v0);
1678 index.resize( numBlocks_+1 );
1679 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1680 newstate.
V = MVT::CloneViewNonConst( *V_, index );
1681 index.resize( keff );
1682 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1683 newstate.
C = MVT::CloneViewNonConst( *C_, index );
1684 newstate.
U = MVT::CloneViewNonConst( *U_, index );
1688 gcrodr_iter->initialize(newstate);
1691 int numRestarts = 0;
1696 gcrodr_iter->iterate();
1703 if ( convTest_->getStatus() ==
Passed ) {
1712 else if ( maxIterTest_->getStatus() ==
Passed ) {
1714 isConverged =
false;
1722 else if ( gcrodr_iter->getCurSubspaceDim() == gcrodr_iter->getMaxSubspaceDim() ) {
1727 RCP<MV> update = gcrodr_iter->getCurrentUpdate();
1728 problem_->updateSolution( update,
true );
1730 buildRecycleSpace2(gcrodr_iter);
1732 printer_->stream(
Debug)
1733 <<
" Generated new recycled subspace using RHS index " 1734 << currIdx[0] <<
" of dimension " << keff << std::endl
1738 if (numRestarts >= maxRestarts_) {
1739 isConverged =
false;
1744 printer_->stream(
Debug)
1745 <<
" Performing restart number " << numRestarts <<
" of " 1746 << maxRestarts_ << std::endl << std::endl;
1749 problem_->computeCurrPrecResVec( &*r_ );
1750 index.resize( 1 ); index[0] = 0;
1751 RCP<MV> v00 = MVT::CloneViewNonConst( *V_, index );
1752 MVT::SetBlock(*r_,index,*v00);
1756 index.resize( numBlocks_+1 );
1757 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1758 restartState.
V = MVT::CloneViewNonConst( *V_, index );
1759 index.resize( keff );
1760 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1761 restartState.
U = MVT::CloneViewNonConst( *U_, index );
1762 restartState.
C = MVT::CloneViewNonConst( *C_, index );
1766 gcrodr_iter->initialize(restartState);
1780 true, std::logic_error,
"Belos::GCRODRSolMgr::solve: " 1781 "Invalid return from GCRODRIter::iterate().");
1786 gcrodr_iter->updateLSQR( gcrodr_iter->getCurSubspaceDim() );
1789 sTest_->checkStatus( &*gcrodr_iter );
1790 if (convTest_->getStatus() !=
Passed)
1791 isConverged =
false;
1794 catch (
const std::exception& e) {
1796 <<
"Error! Caught exception in GCRODRIter::iterate() at iteration " 1797 << gcrodr_iter->getNumIters() << std::endl << e.what() << std::endl;
1804 RCP<MV> update = gcrodr_iter->getCurrentUpdate();
1805 problem_->updateSolution( update,
true );
1808 problem_->setCurrLS();
1813 if (!builtRecycleSpace_) {
1814 buildRecycleSpace2(gcrodr_iter);
1815 printer_->stream(
Debug)
1816 <<
" Generated new recycled subspace using RHS index " << currIdx[0]
1817 <<
" of dimension " << keff << std::endl << std::endl;
1822 if (numRHS2Solve > 0) {
1824 problem_->setLSIndex (currIdx);
1827 currIdx.resize (numRHS2Solve);
1835 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1841 #endif // BELOS_TEUCHOS_TIME_MONITOR 1844 numIters_ = maxIterTest_->getNumIters ();
1856 const std::vector<MagnitudeType>* pTestValues = expConvTest_->getTestValue();
1857 if (pTestValues == NULL || pTestValues->size() < 1) {
1858 pTestValues = impConvTest_->getTestValue();
1861 "Belos::GCRODRSolMgr::solve(): The implicit convergence test's getTestValue() " 1862 "method returned NULL. Please report this bug to the Belos developers.");
1864 "Belos::GCRODRSolMgr::solve(): The implicit convergence test's getTestValue() " 1865 "method returned a vector of length zero. Please report this bug to the " 1866 "Belos developers.");
1871 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1878 template<
class ScalarType,
class MV,
class OP>
1884 std::vector<MagnitudeType> d(keff);
1885 std::vector<ScalarType> dscalar(keff);
1886 std::vector<int> index(numBlocks_+1);
1898 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1901 dscalar.resize(keff);
1902 MVT::MvNorm( *Utmp, d );
1903 for (
int i=0; i<keff; ++i) {
1905 dscalar[i] = (ScalarType)d[i];
1907 MVT::MvScale( *Utmp, dscalar );
1914 for (
int i=0; i<keff; ++i) {
1915 (*H2tmp)(i,i) = d[i];
1923 keff_new = getHarmonicVecs2( keff, p, *H2tmp, oldState.
V, PPtmp );
1932 index.resize( keff );
1933 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1935 index.resize( keff_new );
1936 for (
int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
1937 U1tmp = MVT::CloneViewNonConst( *U1_, index );
1939 MVT::MvTimesMatAddMv( one, *Utmp, PPtmp, zero, *U1tmp );
1945 for (
int ii=0; ii < p; ii++) { index[ii] = ii; }
1948 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
1959 int info = 0, lwork = -1;
1960 tau_.resize (keff_new);
1962 HPtmp.
stride (), &tau_[0], &work_[0], lwork, &info);
1965 "LAPACK's _GEQRF failed to compute a workspace size.");
1972 work_.resize (lwork);
1974 HPtmp.
stride (), &tau_[0], &work_[0], lwork, &info);
1977 "LAPACK's _GEQRF failed to compute a QR factorization.");
1982 for(
int i=0;i<keff_new;i++) {
for(
int j=i;j<keff_new;j++) Rtmp(i,j) = HPtmp(i,j); }
1993 "LAPACK's _UNGQR failed to construct the Q factor.");
2003 for (
int i=0; i < keff; i++) { index[i] = i; }
2005 index.resize(keff_new);
2006 for (
int i=0; i < keff_new; i++) { index[i] = i; }
2007 C1tmp = MVT::CloneViewNonConst( *C1_, index );
2009 MVT::MvTimesMatAddMv( one, *Ctmp, PPtmp, zero, *C1tmp );
2013 index.resize( p+1 );
2014 for (
int i=0; i < p+1; ++i) { index[i] = i; }
2017 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *C1tmp );
2032 work_.resize(lwork);
2037 index.resize(keff_new);
2038 for (
int i=0; i < keff_new; i++) { index[i] = i; }
2040 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
2044 if (keff != keff_new) {
2046 gcrodr_iter->setSize( keff, numBlocks_ );
2056 template<
class ScalarType,
class MV,
class OP>
2061 bool xtraVec =
false;
2066 std::vector<MagnitudeType> wr(m), wi(m);
2072 std::vector<MagnitudeType> w(m);
2075 std::vector<int> iperm(m);
2079 std::vector<ScalarType> work(lwork);
2080 std::vector<MagnitudeType> rwork(lwork);
2086 builtRecycleSpace_ =
true;
2096 ScalarType d = HH(m, m-1) * HH(m, m-1);
2098 for( i=0; i<m; ++i )
2099 harmHH(i, m-1) += d * e_m[i];
2107 lapack.GEEV(
'N',
'V', m, harmHH.values(), harmHH.stride(), &wr[0], &wi[0],
2108 vl, ldvl, vr.
values(), vr.
stride(), &work[0], lwork, &rwork[0], &info);
2112 for( i=0; i<m; ++i )
2116 this->sort(w, m, iperm);
2121 for( i=0; i<recycledBlocks_; ++i ) {
2122 for( j=0; j<m; j++ ) {
2123 PP(j,i) = vr(j,iperm[i]);
2127 if(!scalarTypeIsComplex) {
2130 if (wi[iperm[recycledBlocks_-1]] != 0.0) {
2132 for ( i=0; i<recycledBlocks_; ++i ) {
2133 if (wi[iperm[i]] != 0.0)
2142 if (wi[iperm[recycledBlocks_-1]] > 0.0) {
2143 for( j=0; j<m; ++j ) {
2144 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
2148 for( j=0; j<m; ++j ) {
2149 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
2158 return recycledBlocks_+1;
2161 return recycledBlocks_;
2167 template<
class ScalarType,
class MV,
class OP>
2174 bool xtraVec =
false;
2177 std::vector<int> index;
2180 std::vector<MagnitudeType> wr(m2), wi(m2);
2183 std::vector<MagnitudeType> w(m2);
2189 std::vector<int> iperm(m2);
2192 builtRecycleSpace_ =
true;
2205 index.resize(keffloc);
2206 for (i=0; i<keffloc; ++i) { index[i] = i; }
2210 MVT::MvTransMv( one, *Ctmp, *Utmp, A11 );
2215 for (i=0; i < m+1; i++) { index[i] = i; }
2217 MVT::MvTransMv( one, *Vp, *Utmp, A21 );
2220 for( i=keffloc; i<keffloc+m; i++ ) {
2234 char balanc=
'P', jobvl=
'N', jobvr=
'V', sense=
'N';
2235 int ld =
A.numRows();
2237 int ldvl = ld, ldvr = ld;
2238 int info = 0,ilo = 0,ihi = 0;
2241 std::vector<ScalarType> beta(ld);
2242 std::vector<ScalarType> work(lwork);
2243 std::vector<MagnitudeType> rwork(lwork);
2244 std::vector<MagnitudeType> lscale(ld), rscale(ld);
2245 std::vector<MagnitudeType> rconde(ld), rcondv(ld);
2246 std::vector<int> iwork(ld+6);
2251 lapack.GGEVX(balanc, jobvl, jobvr, sense, ld,
A.values(), ld,
B.values(), ld, &wr[0], &wi[0],
2252 &beta[0], vl, ldvl, vr.
values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
2253 &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &rwork[0],
2254 &iwork[0], bwork, &info);
2259 for( i=0; i<ld; i++ ) {
2265 this->sort(w,ld,iperm);
2270 for( i=0; i<recycledBlocks_; i++ ) {
2271 for( j=0; j<ld; j++ ) {
2272 PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
2276 if(!scalarTypeIsComplex) {
2279 if (wi[iperm[ld-recycledBlocks_]] != 0.0) {
2281 for ( i=ld-recycledBlocks_; i<ld; i++ ) {
2282 if (wi[iperm[i]] != 0.0)
2291 if (wi[iperm[ld-recycledBlocks_]] > 0.0) {
2292 for( j=0; j<ld; j++ ) {
2293 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
2297 for( j=0; j<ld; j++ ) {
2298 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
2307 return recycledBlocks_+1;
2310 return recycledBlocks_;
2317 template<
class ScalarType,
class MV,
class OP>
2319 int l, r, j, i, flag;
2348 if (dlist[j] > dlist[j - 1]) j = j + 1;
2350 if (dlist[j - 1] > dK) {
2351 dlist[i - 1] = dlist[j - 1];
2352 iperm[i - 1] = iperm[j - 1];
2366 dlist[r] = dlist[0];
2367 iperm[r] = iperm[0];
2382 template<
class ScalarType,
class MV,
class OP>
2384 std::ostringstream out;
2387 out <<
"Ortho Type: \"" << orthoType_ <<
"\"";
2388 out <<
", Num Blocks: " <<numBlocks_;
2389 out <<
", Num Recycle Blocks: " << recycledBlocks_;
2390 out <<
", Max Restarts: " << maxRestarts_;
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
Collection of types and exceptions used within the Belos solvers.
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
Teuchos::ScalarTraits< MagnitudeType > MT
Partial specialization for ScalarType types for which Teuchos::LAPACK has a valid implementation...
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
Teuchos::RCP< OutputManager< ScalarType > > printer_
GCRODRSolMgrLAPACKFailure(const std::string &what_arg)
Class which manages the output and verbosity of the Belos solvers.
ScalarType * values() const
bool is_null(const boost::shared_ptr< T > &p)
static const bool requiresLapack
ScaleType
The type of scaling to use on the residual norm value.
Exception thrown to signal error in a status test during Belos::StatusTest::checkStatus().
T & get(ParameterList &l, const std::string &name)
GCRODRSolMgrLAPACKFailure is thrown when a nonzero value is retuned from an LAPACK call...
Teuchos::RCP< std::ostream > outputStream_
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
bool is_null(const std::shared_ptr< T > &p)
MultiVecTraits< ScalarType, MV > MVT
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
Teuchos::ScalarTraits< ScalarType > SCT
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
A factory class for generating StatusTestOutput objects.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H_
An implementation of StatusTestResNorm using a family of residual norms.
GCRODRSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i.e.
GCRODRSolMgrRecyclingFailure(const std::string &what_arg)
This class implements the GCRODR iteration, where a single-std::vector Krylov subspace is constructed...
std::vector< ScalarType > tau_
std::ostream & printValidNames(std::ostream &out) const
Print all recognized MatOrthoManager names to the given ostream.
Belos::StatusTest class for specifying a maximum number of iterations.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > PP_
int curDim
The current dimension of the reduction.
static std::string name()
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > HP_
A factory class for generating StatusTestOutput objects.
Teuchos::RCP< Belos::MatOrthoManager< Scalar, MV, OP > > makeMatOrthoManager(const std::string &ortho, const Teuchos::RCP< const OP > &M, const Teuchos::RCP< OutputManager< Scalar > > &outMan, const std::string &label, const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Return an instance of the specified MatOrthoManager subclass.
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
std::string validNamesString() const
List (as a string) of recognized MatOrthoManager names.
A Belos::StatusTest class for specifying a maximum number of iterations.
GCRODRSolMgrOrthoFailure(const std::string &what_arg)
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
ResetType
How to reset the solver.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Pure virtual base class which describes the basic interface for a solver manager. ...
static constexpr double convTol
Default convergence tolerance.
OrdinalType numRows() const
Teuchos::RCP< MV > V
The current Krylov basis.
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)
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > convTest_
Teuchos::RCP< MV > U
The recycled subspace and its projection.
std::vector< ScalarType > work_
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
Teuchos::RCP< const Teuchos::ParameterList > getDefaultParameters(const std::string &name) const
Default parameters for the given MatOrthoManager subclass.
A linear system to solve, and its associated information.
GCRODRSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate orthonorm...
Class which describes the linear problem to be solved by the iterative solver.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H2_
GCRODRIterOrthoFailure is thrown when the GCRODRIter object is unable to compute independent directio...
Implementation of the GCRODR (Recycling GMRES) iterative linear solver.
bool isType(const std::string &name) const
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
bool isValidName(const std::string &name) const
Whether this factory recognizes the MatOrthoManager with the given name.
ReturnType
Whether the Belos solve converged for all linear systems.
Teuchos::RCP< Teuchos::Time > timerSolve_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B_
Belos concrete class for performing the GCRO-DR iteration.
Structure to contain pointers to GCRODRIter state variables.
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
GCRODRSolMgrRecyclingFailure is thrown when any problem occurs in using/creating the recycling subspa...
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > impConvTest_
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
static magnitudeType magnitude(T a)
virtual ~GCRODRSolMgr()
Destructor.
GCRODRSolMgrLinearProblemFailure(const std::string &what_arg)
GCRODRSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
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< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B
The projection of the Krylov subspace against the recycled subspace.
OrthoManagerFactory< ScalarType, MV, OP > ortho_factory_type
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Get current linear problem being solved for in this object.
A class for extending the status testing capabilities of Belos via logical combinations.
bool isParameter(const std::string &name) const
int getNumIters() const override
Get the iteration count for the most recent call to solve().
Details::SolverManagerRequiresLapack< ScalarType, MV, OP, requiresLapack > base_type
Class which defines basic traits for the operator type.
OperatorTraits< ScalarType, MV, OP > OPT
OrdinalType stride() const
Teuchos::LAPACK< int, ScalarType > lapack
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
Parent class to all Belos exceptions.
MagnitudeType orthoKappa_
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
Base class for Belos::SolverManager subclasses which normally can only compile with ScalarType types ...
Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > ortho_
Orthogonalization manager.
OrdinalType numCols() const
Belos header file which uses auto-configuration information to include necessary C++ headers...
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve...
Teuchos::RCP< Teuchos::ParameterList > params_
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > R_
Belos concrete class for performing the block, flexible GMRES iteration.