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) 303 return Teuchos::tuple(timerSolve_);
347 bool set = problem_->setProblem();
349 throw "Could not set problem.";
393 std::string description()
const;
403 void initializeStateStorage();
412 int getHarmonicVecs1(
int m,
421 int getHarmonicVecs2(
int keff,
int m,
427 void sort(std::vector<MagnitudeType>& dlist,
int n, std::vector<int>& iperm);
524 template<
class ScalarType,
class MV,
class OP>
528 template<
class ScalarType,
class MV,
class OP>
532 template<
class ScalarType,
class MV,
class OP>
535 template<
class ScalarType,
class MV,
class OP>
538 template<
class ScalarType,
class MV,
class OP>
541 template<
class ScalarType,
class MV,
class OP>
544 template<
class ScalarType,
class MV,
class OP>
547 template<
class ScalarType,
class MV,
class OP>
550 template<
class ScalarType,
class MV,
class OP>
553 template<
class ScalarType,
class MV,
class OP>
556 template<
class ScalarType,
class MV,
class OP>
559 template<
class ScalarType,
class MV,
class OP>
562 template<
class ScalarType,
class MV,
class OP>
565 template<
class ScalarType,
class MV,
class OP>
568 template<
class ScalarType,
class MV,
class OP>
574 template<
class ScalarType,
class MV,
class OP>
584 template<
class ScalarType,
class MV,
class OP>
596 problem == Teuchos::null, std::invalid_argument,
597 "Belos::GCRODRSolMgr constructor: The solver manager's " 598 "constructor needs the linear problem argument 'problem' " 613 template<
class ScalarType,
class MV,
class OP>
615 outputStream_ = outputStream_default_;
616 convTol_ = convTol_default_;
617 orthoKappa_ = orthoKappa_default_;
618 maxRestarts_ = maxRestarts_default_;
619 maxIters_ = maxIters_default_;
620 numBlocks_ = numBlocks_default_;
621 recycledBlocks_ = recycledBlocks_default_;
622 verbosity_ = verbosity_default_;
623 outputStyle_ = outputStyle_default_;
624 outputFreq_ = outputFreq_default_;
625 orthoType_ = orthoType_default_;
626 impResScale_ = impResScale_default_;
627 expResScale_ = expResScale_default_;
628 label_ = label_default_;
630 builtRecycleSpace_ =
false;
646 template<
class ScalarType,
class MV,
class OP>
651 using Teuchos::isParameterType;
652 using Teuchos::getParameter;
655 using Teuchos::parameterList;
658 using Teuchos::rcp_dynamic_cast;
659 using Teuchos::rcpFromRef;
666 RCP<const ParameterList> defaultParams = getValidParameters();
684 if (params_.is_null()) {
685 params_ = parameterList (*defaultParams);
693 if (params_ != params) {
699 params_ = parameterList (*params);
734 params_->validateParametersAndSetDefaults (*defaultParams);
739 maxRestarts_ = params->
get(
"Maximum Restarts", maxRestarts_default_);
742 params_->set (
"Maximum Restarts", maxRestarts_);
747 maxIters_ = params->
get (
"Maximum Iterations", maxIters_default_);
750 params_->set (
"Maximum Iterations", maxIters_);
751 if (! maxIterTest_.is_null())
752 maxIterTest_->setMaxIters (maxIters_);
757 numBlocks_ = params->
get (
"Num Blocks", numBlocks_default_);
759 "Belos::GCRODRSolMgr: The \"Num Blocks\" parameter must " 760 "be strictly positive, but you specified a value of " 761 << numBlocks_ <<
".");
763 params_->set (
"Num Blocks", numBlocks_);
768 recycledBlocks_ = params->
get (
"Num Recycled Blocks",
769 recycledBlocks_default_);
771 "Belos::GCRODRSolMgr: The \"Num Recycled Blocks\" " 772 "parameter must be strictly positive, but you specified " 773 "a value of " << recycledBlocks_ <<
".");
775 "Belos::GCRODRSolMgr: The \"Num Recycled Blocks\" " 776 "parameter must be less than the \"Num Blocks\" " 777 "parameter, but you specified \"Num Recycled Blocks\" " 778 "= " << recycledBlocks_ <<
" and \"Num Blocks\" = " 779 << numBlocks_ <<
".");
781 params_->set(
"Num Recycled Blocks", recycledBlocks_);
788 std::string tempLabel = params->
get (
"Timer Label", label_default_);
791 if (tempLabel != label_) {
793 params_->set (
"Timer Label", label_);
794 std::string solveLabel = label_ +
": GCRODRSolMgr total solve time";
795 #ifdef BELOS_TEUCHOS_TIME_MONITOR 798 if (ortho_ != Teuchos::null) {
799 ortho_->setLabel( label_ );
806 if (isParameterType<int> (*params,
"Verbosity")) {
807 verbosity_ = params->
get (
"Verbosity", verbosity_default_);
809 verbosity_ = (int) getParameter<Belos::MsgType> (*params,
"Verbosity");
812 params_->set (
"Verbosity", verbosity_);
815 if (! printer_.is_null())
816 printer_->setVerbosity (verbosity_);
821 if (isParameterType<int> (*params,
"Output Style")) {
822 outputStyle_ = params->
get (
"Output Style", outputStyle_default_);
824 outputStyle_ = (int) getParameter<OutputType> (*params,
"Output Style");
828 params_->set (
"Output Style", outputStyle_);
848 outputStream_ = getParameter<RCP<std::ostream> > (*params,
"Output Stream");
849 }
catch (InvalidParameter&) {
850 outputStream_ = rcpFromRef (std::cout);
857 if (outputStream_.is_null()) {
861 params_->set (
"Output Stream", outputStream_);
864 if (! printer_.is_null()) {
865 printer_->setOStream (outputStream_);
872 outputFreq_ = params->
get (
"Output Frequency", outputFreq_default_);
876 params_->set(
"Output Frequency", outputFreq_);
877 if (! outputTest_.is_null())
878 outputTest_->setOutputFrequency (outputFreq_);
885 if (printer_.is_null()) {
896 bool changedOrthoType =
false;
898 const std::string& tempOrthoType =
899 params->
get (
"Orthogonalization", orthoType_default_);
902 std::ostringstream os;
903 os <<
"Belos::GCRODRSolMgr: Invalid orthogonalization name \"" 904 << tempOrthoType <<
"\". The following are valid options " 905 <<
"for the \"Orthogonalization\" name parameter: ";
907 throw std::invalid_argument (os.str());
909 if (tempOrthoType != orthoType_) {
910 changedOrthoType =
true;
911 orthoType_ = tempOrthoType;
913 params_->set (
"Orthogonalization", orthoType_);
929 RCP<ParameterList> orthoParams;
932 using Teuchos::sublist;
934 const std::string paramName (
"Orthogonalization Parameters");
937 orthoParams = sublist (params_, paramName,
true);
938 }
catch (InvalidParameter&) {
945 orthoParams = sublist (params_, paramName,
true);
949 "Failed to get orthogonalization parameters. " 950 "Please report this bug to the Belos developers.");
955 if (ortho_.is_null() || changedOrthoType) {
961 label_, orthoParams);
970 RCP<PLA> pla = rcp_dynamic_cast<PLA> (ortho_);
976 label_, orthoParams);
978 pla->setParameterList (orthoParams);
990 if (params->
isParameter (
"Orthogonalization Constant")) {
992 params->
get (
"Orthogonalization Constant", orthoKappa_default_);
993 if (orthoKappa > 0) {
994 orthoKappa_ = orthoKappa;
996 params_->set(
"Orthogonalization Constant", orthoKappa_);
998 if (orthoType_ ==
"DGKS" && ! ortho_.is_null()) {
1005 rcp_dynamic_cast<ortho_man_type>(ortho_)->setDepTol (orthoKappa_);
1015 if (params->
isParameter(
"Convergence Tolerance")) {
1016 convTol_ = params->
get (
"Convergence Tolerance", convTol_default_);
1019 params_->set (
"Convergence Tolerance", convTol_);
1020 if (! impConvTest_.is_null())
1021 impConvTest_->setTolerance (convTol_);
1022 if (! expConvTest_.is_null())
1023 expConvTest_->setTolerance (convTol_);
1027 if (params->
isParameter (
"Implicit Residual Scaling")) {
1028 std::string tempImpResScale =
1029 getParameter<std::string> (*params,
"Implicit Residual Scaling");
1032 if (impResScale_ != tempImpResScale) {
1034 impResScale_ = tempImpResScale;
1037 params_->set(
"Implicit Residual Scaling", impResScale_);
1047 if (! impConvTest_.is_null()) {
1053 impConvTest_ = null;
1060 if (params->
isParameter(
"Explicit Residual Scaling")) {
1061 std::string tempExpResScale =
1062 getParameter<std::string> (*params,
"Explicit Residual Scaling");
1065 if (expResScale_ != tempExpResScale) {
1067 expResScale_ = tempExpResScale;
1070 params_->set(
"Explicit Residual Scaling", expResScale_);
1073 if (! expConvTest_.is_null()) {
1079 expConvTest_ = null;
1090 if (maxIterTest_.is_null())
1095 if (impConvTest_.is_null()) {
1096 impConvTest_ =
rcp (
new StatusTestResNorm_t (convTol_));
1102 if (expConvTest_.is_null()) {
1103 expConvTest_ =
rcp (
new StatusTestResNorm_t (convTol_));
1104 expConvTest_->defineResForm (StatusTestResNorm_t::Explicit,
Belos::TwoNorm);
1110 if (convTest_.is_null()) {
1111 convTest_ =
rcp (
new StatusTestCombo_t (StatusTestCombo_t::SEQ,
1119 sTest_ =
rcp (
new StatusTestCombo_t (StatusTestCombo_t::OR,
1125 outputTest_ = stoFactory.
create (printer_, sTest_, outputFreq_,
1129 std::string solverDesc =
" GCRODR ";
1130 outputTest_->setSolverDesc( solverDesc );
1133 if (timerSolve_.is_null()) {
1134 std::string solveLabel = label_ +
": GCRODRSolMgr total solve time";
1135 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1145 template<
class ScalarType,
class MV,
class OP>
1150 using Teuchos::parameterList;
1153 static RCP<const ParameterList> validPL;
1155 RCP<ParameterList> pl = parameterList ();
1158 pl->set(
"Convergence Tolerance", convTol_default_,
1159 "The relative residual tolerance that needs to be achieved by the\n" 1160 "iterative solver in order for the linear system to be declared converged.");
1161 pl->set(
"Maximum Restarts", maxRestarts_default_,
1162 "The maximum number of cycles allowed for each\n" 1163 "set of RHS solved.");
1164 pl->set(
"Maximum Iterations", maxIters_default_,
1165 "The maximum number of iterations allowed for each\n" 1166 "set of RHS solved.");
1170 pl->set(
"Block Size", blockSize_default_,
1171 "Block Size Parameter -- currently must be 1 for GCRODR");
1172 pl->set(
"Num Blocks", numBlocks_default_,
1173 "The maximum number of vectors allowed in the Krylov subspace\n" 1174 "for each set of RHS solved.");
1175 pl->set(
"Num Recycled Blocks", recycledBlocks_default_,
1176 "The maximum number of vectors in the recycled subspace." );
1177 pl->set(
"Verbosity", verbosity_default_,
1178 "What type(s) of solver information should be outputted\n" 1179 "to the output stream.");
1180 pl->set(
"Output Style", outputStyle_default_,
1181 "What style is used for the solver information outputted\n" 1182 "to the output stream.");
1183 pl->set(
"Output Frequency", outputFreq_default_,
1184 "How often convergence information should be outputted\n" 1185 "to the output stream.");
1186 pl->set(
"Output Stream", outputStream_default_,
1187 "A reference-counted pointer to the output stream where all\n" 1188 "solver output is sent.");
1189 pl->set(
"Implicit Residual Scaling", impResScale_default_,
1190 "The type of scaling used in the implicit residual convergence test.");
1191 pl->set(
"Explicit Residual Scaling", expResScale_default_,
1192 "The type of scaling used in the explicit residual convergence test.");
1193 pl->set(
"Timer Label", label_default_,
1194 "The string to use as a prefix for the timer labels.");
1198 pl->set(
"Orthogonalization", orthoType_default_,
1199 "The type of orthogonalization to use. Valid options: " +
1201 RCP<const ParameterList> orthoParams =
1203 pl->
set (
"Orthogonalization Parameters", *orthoParams,
1204 "Parameters specific to the type of orthogonalization used.");
1206 pl->
set(
"Orthogonalization Constant", orthoKappa_default_,
1207 "When using DGKS orthogonalization: the \"depTol\" constant, used " 1208 "to determine whether another step of classical Gram-Schmidt is " 1209 "necessary. Otherwise ignored.");
1216 template<
class ScalarType,
class MV,
class OP>
1223 if (rhsMV == Teuchos::null) {
1231 "Belos::GCRODRSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
1234 if (U_ == Teuchos::null) {
1235 U_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1239 if (MVT::GetNumberVecs(*U_) < recycledBlocks_+1) {
1241 U_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1246 if (C_ == Teuchos::null) {
1247 C_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1251 if (MVT::GetNumberVecs(*C_) < recycledBlocks_+1) {
1253 C_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1258 if (V_ == Teuchos::null) {
1259 V_ = MVT::Clone( *rhsMV, numBlocks_+1 );
1263 if (MVT::GetNumberVecs(*V_) < numBlocks_+1) {
1265 V_ = MVT::Clone( *tmp, numBlocks_+1 );
1270 if (U1_ == Teuchos::null) {
1271 U1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1275 if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1277 U1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1282 if (C1_ == Teuchos::null) {
1283 C1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1287 if (MVT::GetNumberVecs(*C1_) < recycledBlocks_+1) {
1289 C1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1294 if (r_ == Teuchos::null)
1295 r_ = MVT::Clone( *rhsMV, 1 );
1298 tau_.resize(recycledBlocks_+1);
1301 work_.resize(recycledBlocks_+1);
1304 ipiv_.resize(recycledBlocks_+1);
1307 if (H2_ == Teuchos::null)
1310 if ( (H2_->numRows() != numBlocks_+recycledBlocks_+2) || (H2_->numCols() != numBlocks_+recycledBlocks_+1) )
1311 H2_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
1313 H2_->putScalar(zero);
1316 if (R_ == Teuchos::null)
1319 if ( (R_->numRows() != recycledBlocks_+1) || (R_->numCols() != recycledBlocks_+1) )
1320 R_->reshape( recycledBlocks_+1, recycledBlocks_+1 );
1322 R_->putScalar(zero);
1325 if (PP_ == Teuchos::null)
1328 if ( (PP_->numRows() != numBlocks_+recycledBlocks_+2) || (PP_->numCols() != recycledBlocks_+1) )
1329 PP_->reshape( numBlocks_+recycledBlocks_+2, recycledBlocks_+1 );
1333 if (HP_ == Teuchos::null)
1336 if ( (HP_->numRows() != numBlocks_+recycledBlocks_+2) || (HP_->numCols() != numBlocks_+recycledBlocks_+1) )
1337 HP_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
1345 template<
class ScalarType,
class MV,
class OP>
1353 if (!isSet_) { setParameters( params_ ); }
1357 std::vector<int> index(numBlocks_+1);
1364 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1365 std::vector<int> currIdx(1);
1369 problem_->setLSIndex( currIdx );
1372 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1373 if (static_cast<ptrdiff_t>(numBlocks_) > dim) {
1374 numBlocks_ = Teuchos::as<int>(dim);
1376 "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
1377 " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
1378 params_->set(
"Num Blocks", numBlocks_);
1382 bool isConverged =
true;
1385 initializeStateStorage();
1391 plist.
set(
"Num Blocks",numBlocks_);
1392 plist.
set(
"Recycled Blocks",recycledBlocks_);
1397 RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_iter;
1400 int prime_iterations = 0;
1404 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1408 while ( numRHS2Solve > 0 ) {
1411 builtRecycleSpace_ =
false;
1414 outputTest_->reset();
1422 "Belos::GCRODRSolMgr::solve(): Requested size of recycled subspace is not consistent with the current recycle subspace.");
1424 printer_->stream(
Debug) <<
" Now solving RHS index " << currIdx[0] <<
" using recycled subspace of dimension " << keff << std::endl << std::endl;
1427 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1428 RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1429 RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1430 problem_->apply( *Utmp, *Ctmp );
1432 RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1437 int rank = ortho_->normalize(*Ctmp,
rcp(&Rtmp,
false));
1450 work_.resize(lwork);
1455 MVT::MvTimesMatAddMv( one, *Utmp, Rtmp, zero, *U1tmp );
1460 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1461 Ctmp = MVT::CloneViewNonConst( *C_, index );
1462 Utmp = MVT::CloneView( *U_, index );
1466 problem_->computeCurrPrecResVec( &*r_ );
1467 MVT::MvTransMv( one, *Ctmp, *r_, Ctr );
1470 RCP<MV> update = MVT::Clone( *problem_->getCurrLHSVec(), 1 );
1471 MVT::MvInit( *update, 0.0 );
1472 MVT::MvTimesMatAddMv( one, *Utmp, Ctr, one, *update );
1473 problem_->updateSolution( update,
true );
1476 MVT::MvTimesMatAddMv( -one, *Ctmp, Ctr, one, *r_ );
1479 prime_iterations = 0;
1485 printer_->stream(
Debug) <<
" No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
1490 primeList.
set(
"Num Blocks",numBlocks_);
1491 primeList.
set(
"Recycled Blocks",0);
1494 RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_prime_iter;
1498 problem_->computeCurrPrecResVec( &*r_ );
1499 index.resize( 1 ); index[0] = 0;
1500 RCP<MV> v0 = MVT::CloneViewNonConst( *V_, index );
1501 MVT::SetBlock(*r_,index,*v0);
1505 index.resize( numBlocks_+1 );
1506 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1507 newstate.
V = MVT::CloneViewNonConst( *V_, index );
1508 newstate.
U = Teuchos::null;
1509 newstate.
C = Teuchos::null;
1511 newstate.
B = Teuchos::null;
1513 gcrodr_prime_iter->initialize(newstate);
1516 bool primeConverged =
false;
1518 gcrodr_prime_iter->iterate();
1521 if ( convTest_->getStatus() ==
Passed ) {
1523 primeConverged =
true;
1528 gcrodr_prime_iter->updateLSQR( gcrodr_prime_iter->getCurSubspaceDim() );
1531 sTest_->checkStatus( &*gcrodr_prime_iter );
1532 if (convTest_->getStatus() ==
Passed)
1533 primeConverged =
true;
1535 catch (
const std::exception &e) {
1536 printer_->stream(
Errors) <<
"Error! Caught exception in GCRODRIter::iterate() at iteration " 1537 << gcrodr_prime_iter->getNumIters() << std::endl
1538 << e.what() << std::endl;
1542 prime_iterations = gcrodr_prime_iter->getNumIters();
1545 RCP<MV> update = gcrodr_prime_iter->getCurrentUpdate();
1546 problem_->updateSolution( update,
true );
1549 newstate = gcrodr_prime_iter->getState();
1557 if (recycledBlocks_ < p+1) {
1561 keff = getHarmonicVecs1( p, *newstate.
H, *PPtmp );
1566 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1567 RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1568 RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1569 RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1571 for (
int ii=0; ii < p; ++ii) { index[ii] = ii; }
1572 RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1576 MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
1593 HPtmp.
stride (), &tau_[0], &work_[0], lwork, &info);
1596 " LAPACK's _GEQRF failed to compute a workspace size.");
1605 work_.resize (lwork);
1607 HPtmp.
stride (), &tau_[0], &work_[0], lwork, &info);
1610 " LAPACK's _GEQRF failed to compute a QR factorization.");
1615 for (
int ii = 0; ii < keff; ++ii) {
1616 for (
int jj = ii; jj < keff; ++jj) {
1617 Rtmp(ii,jj) = HPtmp(ii,jj);
1629 "LAPACK's _UNGQR failed to construct the Q factor.");
1634 index.resize (p + 1);
1635 for (
int ii = 0; ii < (p+1); ++ii) {
1638 Vtmp = MVT::CloneView( *V_, index );
1639 MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
1650 "LAPACK's _GETRF failed to compute an LU factorization.");
1660 work_.resize(lwork);
1664 "LAPACK's _GETRI failed to invert triangular matrix.");
1667 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
1669 printer_->stream(
Debug)
1670 <<
" Generated recycled subspace using RHS index " << currIdx[0]
1671 <<
" of dimension " << keff << std::endl << std::endl;
1676 if (primeConverged) {
1678 problem_->setCurrLS();
1682 if (numRHS2Solve > 0) {
1684 problem_->setLSIndex (currIdx);
1687 currIdx.resize (numRHS2Solve);
1697 gcrodr_iter->setSize( keff, numBlocks_ );
1700 gcrodr_iter->resetNumIters(prime_iterations);
1703 outputTest_->resetNumCalls();
1706 problem_->computeCurrPrecResVec( &*r_ );
1707 index.resize( 1 ); index[0] = 0;
1708 RCP<MV> v0 = MVT::CloneViewNonConst( *V_, index );
1709 MVT::SetBlock(*r_,index,*v0);
1713 index.resize( numBlocks_+1 );
1714 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1715 newstate.
V = MVT::CloneViewNonConst( *V_, index );
1716 index.resize( keff );
1717 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1718 newstate.
C = MVT::CloneViewNonConst( *C_, index );
1719 newstate.
U = MVT::CloneViewNonConst( *U_, index );
1723 gcrodr_iter->initialize(newstate);
1726 int numRestarts = 0;
1731 gcrodr_iter->iterate();
1738 if ( convTest_->getStatus() ==
Passed ) {
1747 else if ( maxIterTest_->getStatus() ==
Passed ) {
1749 isConverged =
false;
1757 else if ( gcrodr_iter->getCurSubspaceDim() == gcrodr_iter->getMaxSubspaceDim() ) {
1762 RCP<MV> update = gcrodr_iter->getCurrentUpdate();
1763 problem_->updateSolution( update,
true );
1765 buildRecycleSpace2(gcrodr_iter);
1767 printer_->stream(
Debug)
1768 <<
" Generated new recycled subspace using RHS index " 1769 << currIdx[0] <<
" of dimension " << keff << std::endl
1773 if (numRestarts >= maxRestarts_) {
1774 isConverged =
false;
1779 printer_->stream(
Debug)
1780 <<
" Performing restart number " << numRestarts <<
" of " 1781 << maxRestarts_ << std::endl << std::endl;
1784 problem_->computeCurrPrecResVec( &*r_ );
1785 index.resize( 1 ); index[0] = 0;
1786 RCP<MV> v00 = MVT::CloneViewNonConst( *V_, index );
1787 MVT::SetBlock(*r_,index,*v00);
1791 index.resize( numBlocks_+1 );
1792 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1793 restartState.
V = MVT::CloneViewNonConst( *V_, index );
1794 index.resize( keff );
1795 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1796 restartState.
U = MVT::CloneViewNonConst( *U_, index );
1797 restartState.
C = MVT::CloneViewNonConst( *C_, index );
1801 gcrodr_iter->initialize(restartState);
1815 true, std::logic_error,
"Belos::GCRODRSolMgr::solve: " 1816 "Invalid return from GCRODRIter::iterate().");
1821 gcrodr_iter->updateLSQR( gcrodr_iter->getCurSubspaceDim() );
1824 sTest_->checkStatus( &*gcrodr_iter );
1825 if (convTest_->getStatus() !=
Passed)
1826 isConverged =
false;
1829 catch (
const std::exception& e) {
1831 <<
"Error! Caught exception in GCRODRIter::iterate() at iteration " 1832 << gcrodr_iter->getNumIters() << std::endl << e.what() << std::endl;
1839 RCP<MV> update = gcrodr_iter->getCurrentUpdate();
1840 problem_->updateSolution( update,
true );
1843 problem_->setCurrLS();
1848 if (!builtRecycleSpace_) {
1849 buildRecycleSpace2(gcrodr_iter);
1850 printer_->stream(
Debug)
1851 <<
" Generated new recycled subspace using RHS index " << currIdx[0]
1852 <<
" of dimension " << keff << std::endl << std::endl;
1857 if (numRHS2Solve > 0) {
1859 problem_->setLSIndex (currIdx);
1862 currIdx.resize (numRHS2Solve);
1870 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1876 #endif // BELOS_TEUCHOS_TIME_MONITOR 1879 numIters_ = maxIterTest_->getNumIters ();
1891 const std::vector<MagnitudeType>* pTestValues = expConvTest_->getTestValue();
1892 if (pTestValues == NULL || pTestValues->size() < 1) {
1893 pTestValues = impConvTest_->getTestValue();
1896 "Belos::GCRODRSolMgr::solve(): The implicit convergence test's getTestValue() " 1897 "method returned NULL. Please report this bug to the Belos developers.");
1899 "Belos::GCRODRSolMgr::solve(): The implicit convergence test's getTestValue() " 1900 "method returned a vector of length zero. Please report this bug to the " 1901 "Belos developers.");
1906 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1913 template<
class ScalarType,
class MV,
class OP>
1919 std::vector<MagnitudeType> d(keff);
1920 std::vector<ScalarType> dscalar(keff);
1921 std::vector<int> index(numBlocks_+1);
1933 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1936 dscalar.resize(keff);
1937 MVT::MvNorm( *Utmp, d );
1938 for (
int i=0; i<keff; ++i) {
1940 dscalar[i] = (ScalarType)d[i];
1942 MVT::MvScale( *Utmp, dscalar );
1949 for (
int i=0; i<keff; ++i) {
1950 (*H2tmp)(i,i) = d[i];
1958 keff_new = getHarmonicVecs2( keff, p, *H2tmp, oldState.
V, PPtmp );
1967 index.resize( keff );
1968 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1970 index.resize( keff_new );
1971 for (
int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
1972 U1tmp = MVT::CloneViewNonConst( *U1_, index );
1974 MVT::MvTimesMatAddMv( one, *Utmp, PPtmp, zero, *U1tmp );
1980 for (
int ii=0; ii < p; ii++) { index[ii] = ii; }
1983 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
1994 int info = 0, lwork = -1;
1995 tau_.resize (keff_new);
1997 HPtmp.
stride (), &tau_[0], &work_[0], lwork, &info);
2000 "LAPACK's _GEQRF failed to compute a workspace size.");
2007 work_.resize (lwork);
2009 HPtmp.
stride (), &tau_[0], &work_[0], lwork, &info);
2012 "LAPACK's _GEQRF failed to compute a QR factorization.");
2017 for(
int i=0;i<keff_new;i++) {
for(
int j=i;j<keff_new;j++) Rtmp(i,j) = HPtmp(i,j); }
2028 "LAPACK's _UNGQR failed to construct the Q factor.");
2038 for (
int i=0; i < keff; i++) { index[i] = i; }
2040 index.resize(keff_new);
2041 for (
int i=0; i < keff_new; i++) { index[i] = i; }
2042 C1tmp = MVT::CloneViewNonConst( *C1_, index );
2044 MVT::MvTimesMatAddMv( one, *Ctmp, PPtmp, zero, *C1tmp );
2048 index.resize( p+1 );
2049 for (
int i=0; i < p+1; ++i) { index[i] = i; }
2052 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *C1tmp );
2067 work_.resize(lwork);
2072 index.resize(keff_new);
2073 for (
int i=0; i < keff_new; i++) { index[i] = i; }
2075 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
2079 if (keff != keff_new) {
2081 gcrodr_iter->setSize( keff, numBlocks_ );
2091 template<
class ScalarType,
class MV,
class OP>
2096 bool xtraVec =
false;
2101 std::vector<MagnitudeType> wr(m), wi(m);
2107 std::vector<MagnitudeType> w(m);
2110 std::vector<int> iperm(m);
2114 std::vector<ScalarType> work(lwork);
2115 std::vector<MagnitudeType> rwork(lwork);
2121 builtRecycleSpace_ =
true;
2131 ScalarType d = HH(m, m-1) * HH(m, m-1);
2133 for( i=0; i<m; ++i )
2134 harmHH(i, m-1) += d * e_m[i];
2142 lapack.GEEV(
'N',
'V', m, harmHH.values(), harmHH.stride(), &wr[0], &wi[0],
2143 vl, ldvl, vr.
values(), vr.
stride(), &work[0], lwork, &rwork[0], &info);
2147 for( i=0; i<m; ++i )
2151 this->sort(w, m, iperm);
2156 for( i=0; i<recycledBlocks_; ++i ) {
2157 for( j=0; j<m; j++ ) {
2158 PP(j,i) = vr(j,iperm[i]);
2162 if(!scalarTypeIsComplex) {
2165 if (wi[iperm[recycledBlocks_-1]] != 0.0) {
2167 for ( i=0; i<recycledBlocks_; ++i ) {
2168 if (wi[iperm[i]] != 0.0)
2177 if (wi[iperm[recycledBlocks_-1]] > 0.0) {
2178 for( j=0; j<m; ++j ) {
2179 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
2183 for( j=0; j<m; ++j ) {
2184 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
2193 return recycledBlocks_+1;
2196 return recycledBlocks_;
2202 template<
class ScalarType,
class MV,
class OP>
2209 bool xtraVec =
false;
2212 std::vector<int> index;
2215 std::vector<MagnitudeType> wr(m2), wi(m2);
2218 std::vector<MagnitudeType> w(m2);
2224 std::vector<int> iperm(m2);
2227 builtRecycleSpace_ =
true;
2240 index.resize(keffloc);
2241 for (i=0; i<keffloc; ++i) { index[i] = i; }
2245 MVT::MvTransMv( one, *Ctmp, *Utmp, A11 );
2250 for (i=0; i < m+1; i++) { index[i] = i; }
2252 MVT::MvTransMv( one, *Vp, *Utmp, A21 );
2255 for( i=keffloc; i<keffloc+m; i++ ) {
2269 char balanc=
'P', jobvl=
'N', jobvr=
'V', sense=
'N';
2270 int ld =
A.numRows();
2272 int ldvl = ld, ldvr = ld;
2273 int info = 0,ilo = 0,ihi = 0;
2276 std::vector<ScalarType> beta(ld);
2277 std::vector<ScalarType> work(lwork);
2278 std::vector<MagnitudeType> rwork(lwork);
2279 std::vector<MagnitudeType> lscale(ld), rscale(ld);
2280 std::vector<MagnitudeType> rconde(ld), rcondv(ld);
2281 std::vector<int> iwork(ld+6);
2286 lapack.GGEVX(balanc, jobvl, jobvr, sense, ld,
A.values(), ld,
B.values(), ld, &wr[0], &wi[0],
2287 &beta[0], vl, ldvl, vr.
values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
2288 &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &rwork[0],
2289 &iwork[0], bwork, &info);
2294 for( i=0; i<ld; i++ ) {
2300 this->sort(w,ld,iperm);
2305 for( i=0; i<recycledBlocks_; i++ ) {
2306 for( j=0; j<ld; j++ ) {
2307 PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
2311 if(!scalarTypeIsComplex) {
2314 if (wi[iperm[ld-recycledBlocks_]] != 0.0) {
2316 for ( i=ld-recycledBlocks_; i<ld; i++ ) {
2317 if (wi[iperm[i]] != 0.0)
2326 if (wi[iperm[ld-recycledBlocks_]] > 0.0) {
2327 for( j=0; j<ld; j++ ) {
2328 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
2332 for( j=0; j<ld; j++ ) {
2333 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
2342 return recycledBlocks_+1;
2345 return recycledBlocks_;
2352 template<
class ScalarType,
class MV,
class OP>
2354 int l, r, j, i, flag;
2383 if (dlist[j] > dlist[j - 1]) j = j + 1;
2385 if (dlist[j - 1] > dK) {
2386 dlist[i - 1] = dlist[j - 1];
2387 iperm[i - 1] = iperm[j - 1];
2401 dlist[r] = dlist[0];
2402 iperm[r] = iperm[0];
2417 template<
class ScalarType,
class MV,
class OP>
2419 std::ostringstream out;
2422 out <<
"Ortho Type: \"" << orthoType_ <<
"\"";
2423 out <<
", Num Blocks: " <<numBlocks_;
2424 out <<
", Num Recycle Blocks: " << recycledBlocks_;
2425 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
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 std::string impResScale_default_
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().
static const int outputFreq_default_
T & get(ParameterList &l, const std::string &name)
GCRODRSolMgrLAPACKFailure is thrown when a nonzero value is retuned from an LAPACK call...
static const std::string expResScale_default_
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
MagnitudeType achievedTol() const
Tolerance achieved by the last solve() invocation.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
A factory class for generating StatusTestOutput objects.
static const int outputStyle_default_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H_
static const std::string label_default_
static const MagnitudeType orthoKappa_default_
void reset(const ResetType type)
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
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)
static const int maxRestarts_default_
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 const MagnitudeType convTol_default_
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.
bool isLOADetected() const
Return whether a loss of accuracy was detected by this solver during the most current solve...
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
static const int recycledBlocks_default_
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)
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. ...
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.
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.
static const std::string orthoType_default_
Teuchos::RCP< Teuchos::Time > timerSolve_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B_
Belos concrete class for performing the GCRO-DR iteration.
static const Teuchos::RCP< std::ostream > outputStream_default_
Structure to contain pointers to GCRODRIter state variables.
GCRODRSolMgrRecyclingFailure is thrown when any problem occurs in using/creating the recycling subspa...
int getNumIters() const
Get the iteration count for the most recent call to solve().
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)
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
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.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem)
Set the linear problem that needs to be solved.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
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...
A class for extending the status testing capabilities of Belos via logical combinations.
bool isParameter(const std::string &name) const
Details::SolverManagerRequiresLapack< ScalarType, MV, OP, requiresLapack > base_type
static const int blockSize_default_
static const int verbosity_default_
Class which defines basic traits for the operator type.
static const int numBlocks_default_
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...
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get current linear problem being solved for in this object.
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_
static const int maxIters_default_
Belos concrete class for performing the block, flexible GMRES iteration.