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>
184 #if defined(HAVE_TEUCHOSCORE_CXX11)
185 # if defined(HAVE_TEUCHOS_COMPLEX)
186 #if defined(HAVE_TEUCHOS_LONG_DOUBLE)
187 static_assert (std::is_same<ScalarType, std::complex<float> >::value ||
188 std::is_same<ScalarType, std::complex<double> >::value ||
189 std::is_same<ScalarType, Kokkos::complex<double> >::value ||
190 std::is_same<ScalarType, float>::value ||
191 std::is_same<ScalarType, double>::value ||
192 std::is_same<ScalarType, long double>::value,
193 "Belos::GCRODRSolMgr: ScalarType must be one of the four "
194 "types (S,D,C,Z) supported by LAPACK or long double (largely not impl'd).");
196 static_assert (std::is_same<ScalarType, std::complex<float> >::value ||
197 std::is_same<ScalarType, std::complex<double> >::value ||
198 std::is_same<ScalarType, Kokkos::complex<double> >::value ||
199 std::is_same<ScalarType, float>::value ||
200 std::is_same<ScalarType, double>::value,
201 "Belos::GCRODRSolMgr: ScalarType must be one of the four "
202 "types (S,D,C,Z) supported by LAPACK.");
205 #if defined(HAVE_TEUCHOS_LONG_DOUBLE)
206 static_assert (std::is_same<ScalarType, float>::value ||
207 std::is_same<ScalarType, double>::value ||
208 std::is_same<ScalarType, long double>::value,
209 "Belos::GCRODRSolMgr: ScalarType must be float, double or long double. "
210 "Complex arithmetic support is currently disabled. To "
211 "enable it, set Teuchos_ENABLE_COMPLEX=ON.");
213 static_assert (std::is_same<ScalarType, float>::value ||
214 std::is_same<ScalarType, double>::value,
215 "Belos::GCRODRSolMgr: ScalarType must be float or double. "
216 "Complex arithmetic support is currently disabled. To "
217 "enable it, set Teuchos_ENABLE_COMPLEX=ON.");
219 # endif // defined(HAVE_TEUCHOS_COMPLEX)
220 #endif // defined(HAVE_TEUCHOSCORE_CXX11)
330 return Teuchos::tuple(timerSolve_);
374 bool set = problem_->setProblem();
376 throw "Could not set problem.";
420 std::string description()
const override;
430 void initializeStateStorage();
439 int getHarmonicVecs1(
int m,
448 int getHarmonicVecs2(
int keff,
int m,
454 void sort(std::vector<MagnitudeType>& dlist,
int n, std::vector<int>& iperm);
482 static constexpr
double orthoKappa_default_ = 0.0;
483 static constexpr
int maxRestarts_default_ = 100;
484 static constexpr
int maxIters_default_ = 1000;
485 static constexpr
int numBlocks_default_ = 50;
486 static constexpr
int blockSize_default_ = 1;
487 static constexpr
int recycledBlocks_default_ = 5;
490 static constexpr
int outputFreq_default_ = -1;
491 static constexpr
const char * impResScale_default_ =
"Norm of Preconditioned Initial Residual";
492 static constexpr
const char * expResScale_default_ =
"Norm of Initial Residual";
493 static constexpr
const char * label_default_ =
"Belos";
494 static constexpr
const char * orthoType_default_ =
"ICGS";
549 template<
class ScalarType,
class MV,
class OP>
559 template<
class ScalarType,
class MV,
class OP>
572 "Belos::GCRODRSolMgr constructor: The solver manager's "
573 "constructor needs the linear problem argument 'problem' "
588 template<
class ScalarType,
class MV,
class OP>
590 outputStream_ = Teuchos::rcpFromRef(std::cout);
592 orthoKappa_ = orthoKappa_default_;
593 maxRestarts_ = maxRestarts_default_;
594 maxIters_ = maxIters_default_;
595 numBlocks_ = numBlocks_default_;
596 recycledBlocks_ = recycledBlocks_default_;
597 verbosity_ = verbosity_default_;
598 outputStyle_ = outputStyle_default_;
599 outputFreq_ = outputFreq_default_;
600 orthoType_ = orthoType_default_;
601 impResScale_ = impResScale_default_;
602 expResScale_ = expResScale_default_;
603 label_ = label_default_;
605 builtRecycleSpace_ =
false;
621 template<
class ScalarType,
class MV,
class OP>
626 using Teuchos::isParameterType;
627 using Teuchos::getParameter;
630 using Teuchos::parameterList;
633 using Teuchos::rcp_dynamic_cast;
634 using Teuchos::rcpFromRef;
641 RCP<const ParameterList> defaultParams = getValidParameters();
659 if (params_.is_null()) {
660 params_ = parameterList (*defaultParams);
668 if (params_ != params) {
674 params_ = parameterList (*params);
709 params_->validateParametersAndSetDefaults (*defaultParams);
714 maxRestarts_ = params->
get(
"Maximum Restarts", maxRestarts_default_);
717 params_->set (
"Maximum Restarts", maxRestarts_);
722 maxIters_ = params->
get (
"Maximum Iterations", maxIters_default_);
725 params_->set (
"Maximum Iterations", maxIters_);
726 if (! maxIterTest_.is_null())
727 maxIterTest_->setMaxIters (maxIters_);
732 numBlocks_ = params->
get (
"Num Blocks", numBlocks_default_);
734 "Belos::GCRODRSolMgr: The \"Num Blocks\" parameter must "
735 "be strictly positive, but you specified a value of "
736 << numBlocks_ <<
".");
738 params_->set (
"Num Blocks", numBlocks_);
743 recycledBlocks_ = params->
get (
"Num Recycled Blocks",
744 recycledBlocks_default_);
746 "Belos::GCRODRSolMgr: The \"Num Recycled Blocks\" "
747 "parameter must be strictly positive, but you specified "
748 "a value of " << recycledBlocks_ <<
".");
750 "Belos::GCRODRSolMgr: The \"Num Recycled Blocks\" "
751 "parameter must be less than the \"Num Blocks\" "
752 "parameter, but you specified \"Num Recycled Blocks\" "
753 "= " << recycledBlocks_ <<
" and \"Num Blocks\" = "
754 << numBlocks_ <<
".");
756 params_->set(
"Num Recycled Blocks", recycledBlocks_);
763 std::string tempLabel = params->
get (
"Timer Label", label_default_);
766 if (tempLabel != label_) {
768 params_->set (
"Timer Label", label_);
769 std::string solveLabel = label_ +
": GCRODRSolMgr total solve time";
770 #ifdef BELOS_TEUCHOS_TIME_MONITOR
774 ortho_->setLabel( label_ );
781 if (isParameterType<int> (*params,
"Verbosity")) {
782 verbosity_ = params->
get (
"Verbosity", verbosity_default_);
784 verbosity_ = (int) getParameter<Belos::MsgType> (*params,
"Verbosity");
787 params_->set (
"Verbosity", verbosity_);
790 if (! printer_.is_null())
791 printer_->setVerbosity (verbosity_);
796 if (isParameterType<int> (*params,
"Output Style")) {
797 outputStyle_ = params->
get (
"Output Style", outputStyle_default_);
799 outputStyle_ = (int) getParameter<OutputType> (*params,
"Output Style");
803 params_->set (
"Output Style", outputStyle_);
823 outputStream_ = getParameter<RCP<std::ostream> > (*params,
"Output Stream");
824 }
catch (InvalidParameter&) {
825 outputStream_ = rcpFromRef (std::cout);
832 if (outputStream_.is_null()) {
836 params_->set (
"Output Stream", outputStream_);
839 if (! printer_.is_null()) {
840 printer_->setOStream (outputStream_);
847 outputFreq_ = params->
get (
"Output Frequency", outputFreq_default_);
851 params_->set(
"Output Frequency", outputFreq_);
852 if (! outputTest_.is_null())
853 outputTest_->setOutputFrequency (outputFreq_);
860 if (printer_.is_null()) {
871 bool changedOrthoType =
false;
873 const std::string& tempOrthoType =
874 params->
get (
"Orthogonalization", orthoType_default_);
877 std::ostringstream os;
878 os <<
"Belos::GCRODRSolMgr: Invalid orthogonalization name \""
879 << tempOrthoType <<
"\". The following are valid options "
880 <<
"for the \"Orthogonalization\" name parameter: ";
882 throw std::invalid_argument (os.str());
884 if (tempOrthoType != orthoType_) {
885 changedOrthoType =
true;
886 orthoType_ = tempOrthoType;
888 params_->set (
"Orthogonalization", orthoType_);
904 RCP<ParameterList> orthoParams;
907 using Teuchos::sublist;
909 const std::string paramName (
"Orthogonalization Parameters");
912 orthoParams = sublist (params_, paramName,
true);
913 }
catch (InvalidParameter&) {
920 orthoParams = sublist (params_, paramName,
true);
924 "Failed to get orthogonalization parameters. "
925 "Please report this bug to the Belos developers.");
930 if (ortho_.is_null() || changedOrthoType) {
936 label_, orthoParams);
945 RCP<PLA> pla = rcp_dynamic_cast<PLA> (ortho_);
951 label_, orthoParams);
953 pla->setParameterList (orthoParams);
965 if (params->
isParameter (
"Orthogonalization Constant")) {
968 orthoKappa = params->
get (
"Orthogonalization Constant", orthoKappa);
971 orthoKappa = params->
get (
"Orthogonalization Constant", orthoKappa_default_);
974 if (orthoKappa > 0) {
975 orthoKappa_ = orthoKappa;
977 params_->set(
"Orthogonalization Constant", orthoKappa_);
979 if (orthoType_ ==
"DGKS" && ! ortho_.is_null()) {
986 rcp_dynamic_cast<ortho_man_type>(ortho_)->setDepTol (orthoKappa_);
996 if (params->
isParameter(
"Convergence Tolerance")) {
998 convTol_ = params->
get (
"Convergence Tolerance",
1006 params_->set (
"Convergence Tolerance", convTol_);
1007 if (! impConvTest_.is_null())
1008 impConvTest_->setTolerance (convTol_);
1009 if (! expConvTest_.is_null())
1010 expConvTest_->setTolerance (convTol_);
1014 if (params->
isParameter (
"Implicit Residual Scaling")) {
1015 std::string tempImpResScale =
1016 getParameter<std::string> (*params,
"Implicit Residual Scaling");
1019 if (impResScale_ != tempImpResScale) {
1021 impResScale_ = tempImpResScale;
1024 params_->set(
"Implicit Residual Scaling", impResScale_);
1034 if (! impConvTest_.is_null()) {
1040 impConvTest_ = null;
1047 if (params->
isParameter(
"Explicit Residual Scaling")) {
1048 std::string tempExpResScale =
1049 getParameter<std::string> (*params,
"Explicit Residual Scaling");
1052 if (expResScale_ != tempExpResScale) {
1054 expResScale_ = tempExpResScale;
1057 params_->set(
"Explicit Residual Scaling", expResScale_);
1060 if (! expConvTest_.is_null()) {
1066 expConvTest_ = null;
1077 if (maxIterTest_.is_null())
1082 if (impConvTest_.is_null()) {
1083 impConvTest_ =
rcp (
new StatusTestResNorm_t (convTol_));
1089 if (expConvTest_.is_null()) {
1090 expConvTest_ =
rcp (
new StatusTestResNorm_t (convTol_));
1091 expConvTest_->defineResForm (StatusTestResNorm_t::Explicit,
Belos::TwoNorm);
1097 if (convTest_.is_null()) {
1098 convTest_ =
rcp (
new StatusTestCombo_t (StatusTestCombo_t::SEQ,
1106 sTest_ =
rcp (
new StatusTestCombo_t (StatusTestCombo_t::OR,
1112 outputTest_ = stoFactory.
create (printer_, sTest_, outputFreq_,
1116 std::string solverDesc =
" GCRODR ";
1117 outputTest_->setSolverDesc( solverDesc );
1120 if (timerSolve_.is_null()) {
1121 std::string solveLabel = label_ +
": GCRODRSolMgr total solve time";
1122 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1132 template<
class ScalarType,
class MV,
class OP>
1137 using Teuchos::parameterList;
1140 static RCP<const ParameterList> validPL;
1142 RCP<ParameterList> pl = parameterList ();
1146 "The relative residual tolerance that needs to be achieved by the\n"
1147 "iterative solver in order for the linear system to be declared converged.");
1148 pl->set(
"Maximum Restarts", static_cast<int>(maxRestarts_default_),
1149 "The maximum number of cycles allowed for each\n"
1150 "set of RHS solved.");
1151 pl->set(
"Maximum Iterations", static_cast<int>(maxIters_default_),
1152 "The maximum number of iterations allowed for each\n"
1153 "set of RHS solved.");
1157 pl->set(
"Block Size", static_cast<int>(blockSize_default_),
1158 "Block Size Parameter -- currently must be 1 for GCRODR");
1159 pl->set(
"Num Blocks", static_cast<int>(numBlocks_default_),
1160 "The maximum number of vectors allowed in the Krylov subspace\n"
1161 "for each set of RHS solved.");
1162 pl->set(
"Num Recycled Blocks", static_cast<int>(recycledBlocks_default_),
1163 "The maximum number of vectors in the recycled subspace." );
1164 pl->set(
"Verbosity", static_cast<int>(verbosity_default_),
1165 "What type(s) of solver information should be outputted\n"
1166 "to the output stream.");
1167 pl->set(
"Output Style", static_cast<int>(outputStyle_default_),
1168 "What style is used for the solver information outputted\n"
1169 "to the output stream.");
1170 pl->set(
"Output Frequency", static_cast<int>(outputFreq_default_),
1171 "How often convergence information should be outputted\n"
1172 "to the output stream.");
1173 pl->set(
"Output Stream", Teuchos::rcpFromRef(std::cout),
1174 "A reference-counted pointer to the output stream where all\n"
1175 "solver output is sent.");
1176 pl->set(
"Implicit Residual Scaling", static_cast<const char *>(impResScale_default_),
1177 "The type of scaling used in the implicit residual convergence test.");
1178 pl->set(
"Explicit Residual Scaling", static_cast<const char *>(expResScale_default_),
1179 "The type of scaling used in the explicit residual convergence test.");
1180 pl->set(
"Timer Label", static_cast<const char *>(label_default_),
1181 "The string to use as a prefix for the timer labels.");
1184 pl->set(
"Orthogonalization", static_cast<const char *>(orthoType_default_),
1185 "The type of orthogonalization to use. Valid options: " +
1187 RCP<const ParameterList> orthoParams =
1189 pl->
set (
"Orthogonalization Parameters", *orthoParams,
1190 "Parameters specific to the type of orthogonalization used.");
1192 pl->
set(
"Orthogonalization Constant",static_cast<MagnitudeType>(orthoKappa_default_),
1193 "When using DGKS orthogonalization: the \"depTol\" constant, used "
1194 "to determine whether another step of classical Gram-Schmidt is "
1195 "necessary. Otherwise ignored.");
1202 template<
class ScalarType,
class MV,
class OP>
1217 "Belos::GCRODRSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
1221 U_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1225 if (MVT::GetNumberVecs(*U_) < recycledBlocks_+1) {
1227 U_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1233 C_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1237 if (MVT::GetNumberVecs(*C_) < recycledBlocks_+1) {
1239 C_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1245 V_ = MVT::Clone( *rhsMV, numBlocks_+1 );
1249 if (MVT::GetNumberVecs(*V_) < numBlocks_+1) {
1251 V_ = MVT::Clone( *tmp, numBlocks_+1 );
1257 U1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1261 if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1263 U1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1269 C1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1273 if (MVT::GetNumberVecs(*C1_) < recycledBlocks_+1) {
1275 C1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1281 r_ = MVT::Clone( *rhsMV, 1 );
1284 tau_.resize(recycledBlocks_+1);
1287 work_.resize(recycledBlocks_+1);
1290 ipiv_.resize(recycledBlocks_+1);
1296 if ( (H2_->numRows() != numBlocks_+recycledBlocks_+2) || (H2_->numCols() != numBlocks_+recycledBlocks_+1) )
1297 H2_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
1299 H2_->putScalar(zero);
1305 if ( (R_->numRows() != recycledBlocks_+1) || (R_->numCols() != recycledBlocks_+1) )
1306 R_->reshape( recycledBlocks_+1, recycledBlocks_+1 );
1308 R_->putScalar(zero);
1314 if ( (PP_->numRows() != numBlocks_+recycledBlocks_+2) || (PP_->numCols() != recycledBlocks_+1) )
1315 PP_->reshape( numBlocks_+recycledBlocks_+2, recycledBlocks_+1 );
1322 if ( (HP_->numRows() != numBlocks_+recycledBlocks_+2) || (HP_->numCols() != numBlocks_+recycledBlocks_+1) )
1323 HP_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
1331 template<
class ScalarType,
class MV,
class OP>
1339 if (!isSet_) { setParameters( params_ ); }
1343 std::vector<int> index(numBlocks_+1);
1350 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1351 std::vector<int> currIdx(1);
1355 problem_->setLSIndex( currIdx );
1358 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1359 if (static_cast<ptrdiff_t>(numBlocks_) > dim) {
1360 numBlocks_ = Teuchos::as<int>(dim);
1362 "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
1363 " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
1364 params_->set(
"Num Blocks", numBlocks_);
1368 bool isConverged =
true;
1371 initializeStateStorage();
1377 plist.
set(
"Num Blocks",numBlocks_);
1378 plist.
set(
"Recycled Blocks",recycledBlocks_);
1383 RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_iter;
1386 int prime_iterations = 0;
1390 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1394 while ( numRHS2Solve > 0 ) {
1397 builtRecycleSpace_ =
false;
1400 outputTest_->reset();
1408 "Belos::GCRODRSolMgr::solve(): Requested size of recycled subspace is not consistent with the current recycle subspace.");
1410 printer_->stream(
Debug) <<
" Now solving RHS index " << currIdx[0] <<
" using recycled subspace of dimension " << keff << std::endl << std::endl;
1413 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1414 RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1415 RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1416 problem_->apply( *Utmp, *Ctmp );
1418 RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1423 int rank = ortho_->normalize(*Ctmp,
rcp(&Rtmp,
false));
1436 work_.resize(lwork);
1441 MVT::MvTimesMatAddMv( one, *Utmp, Rtmp, zero, *U1tmp );
1446 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1447 Ctmp = MVT::CloneViewNonConst( *C_, index );
1448 Utmp = MVT::CloneView( *U_, index );
1452 problem_->computeCurrPrecResVec( &*r_ );
1453 MVT::MvTransMv( one, *Ctmp, *r_, Ctr );
1456 RCP<MV> update = MVT::Clone( *problem_->getCurrLHSVec(), 1 );
1457 MVT::MvInit( *update, 0.0 );
1458 MVT::MvTimesMatAddMv( one, *Utmp, Ctr, one, *update );
1459 problem_->updateSolution( update,
true );
1462 MVT::MvTimesMatAddMv( -one, *Ctmp, Ctr, one, *r_ );
1465 prime_iterations = 0;
1471 printer_->stream(
Debug) <<
" No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
1476 primeList.
set(
"Num Blocks",numBlocks_);
1477 primeList.
set(
"Recycled Blocks",0);
1480 RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_prime_iter;
1484 problem_->computeCurrPrecResVec( &*r_ );
1485 index.resize( 1 ); index[0] = 0;
1486 RCP<MV> v0 = MVT::CloneViewNonConst( *V_, index );
1487 MVT::SetBlock(*r_,index,*v0);
1491 index.resize( numBlocks_+1 );
1492 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1493 newstate.
V = MVT::CloneViewNonConst( *V_, index );
1499 gcrodr_prime_iter->initialize(newstate);
1502 bool primeConverged =
false;
1504 gcrodr_prime_iter->iterate();
1507 if ( convTest_->getStatus() ==
Passed ) {
1509 primeConverged =
true;
1514 gcrodr_prime_iter->updateLSQR( gcrodr_prime_iter->getCurSubspaceDim() );
1517 sTest_->checkStatus( &*gcrodr_prime_iter );
1518 if (convTest_->getStatus() ==
Passed)
1519 primeConverged =
true;
1523 achievedTol_ = MT::one();
1525 MVT::MvInit( *X, SCT::zero() );
1526 printer_->stream(
Warnings) <<
"Belos::GCRODRSolMgr::solve(): Warning! NaN has been detected!"
1530 catch (
const std::exception &e) {
1531 printer_->stream(
Errors) <<
"Error! Caught exception in GCRODRIter::iterate() at iteration "
1532 << gcrodr_prime_iter->getNumIters() << std::endl
1533 << e.what() << std::endl;
1537 prime_iterations = gcrodr_prime_iter->getNumIters();
1540 RCP<MV> update = gcrodr_prime_iter->getCurrentUpdate();
1541 problem_->updateSolution( update,
true );
1544 newstate = gcrodr_prime_iter->getState();
1552 if (recycledBlocks_ < p+1) {
1556 keff = getHarmonicVecs1( p, *newstate.
H, *PPtmp );
1561 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1562 RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1563 RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1564 RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1566 for (
int ii=0; ii < p; ++ii) { index[ii] = ii; }
1567 RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1571 MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
1588 HPtmp.
stride (), &tau_[0], &work_[0], lwork, &info);
1591 " LAPACK's _GEQRF failed to compute a workspace size.");
1600 work_.resize (lwork);
1602 HPtmp.
stride (), &tau_[0], &work_[0], lwork, &info);
1605 " LAPACK's _GEQRF failed to compute a QR factorization.");
1610 for (
int ii = 0; ii < keff; ++ii) {
1611 for (
int jj = ii; jj < keff; ++jj) {
1612 Rtmp(ii,jj) = HPtmp(ii,jj);
1624 "LAPACK's _UNGQR failed to construct the Q factor.");
1629 index.resize (p + 1);
1630 for (
int ii = 0; ii < (p+1); ++ii) {
1633 Vtmp = MVT::CloneView( *V_, index );
1634 MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
1645 "LAPACK's _GETRF failed to compute an LU factorization.");
1655 work_.resize(lwork);
1659 "LAPACK's _GETRI failed to invert triangular matrix.");
1662 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
1664 printer_->stream(
Debug)
1665 <<
" Generated recycled subspace using RHS index " << currIdx[0]
1666 <<
" of dimension " << keff << std::endl << std::endl;
1671 if (primeConverged) {
1673 problem_->setCurrLS();
1677 if (numRHS2Solve > 0) {
1679 problem_->setLSIndex (currIdx);
1682 currIdx.resize (numRHS2Solve);
1692 gcrodr_iter->setSize( keff, numBlocks_ );
1695 gcrodr_iter->resetNumIters(prime_iterations);
1698 outputTest_->resetNumCalls();
1701 problem_->computeCurrPrecResVec( &*r_ );
1702 index.resize( 1 ); index[0] = 0;
1703 RCP<MV> v0 = MVT::CloneViewNonConst( *V_, index );
1704 MVT::SetBlock(*r_,index,*v0);
1708 index.resize( numBlocks_+1 );
1709 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1710 newstate.
V = MVT::CloneViewNonConst( *V_, index );
1711 index.resize( keff );
1712 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1713 newstate.
C = MVT::CloneViewNonConst( *C_, index );
1714 newstate.
U = MVT::CloneViewNonConst( *U_, index );
1718 gcrodr_iter->initialize(newstate);
1721 int numRestarts = 0;
1726 gcrodr_iter->iterate();
1733 if ( convTest_->getStatus() ==
Passed ) {
1742 else if ( maxIterTest_->getStatus() ==
Passed ) {
1744 isConverged =
false;
1752 else if ( gcrodr_iter->getCurSubspaceDim() == gcrodr_iter->getMaxSubspaceDim() ) {
1757 RCP<MV> update = gcrodr_iter->getCurrentUpdate();
1758 problem_->updateSolution( update,
true );
1760 buildRecycleSpace2(gcrodr_iter);
1762 printer_->stream(
Debug)
1763 <<
" Generated new recycled subspace using RHS index "
1764 << currIdx[0] <<
" of dimension " << keff << std::endl
1768 if (numRestarts >= maxRestarts_) {
1769 isConverged =
false;
1774 printer_->stream(
Debug)
1775 <<
" Performing restart number " << numRestarts <<
" of "
1776 << maxRestarts_ << std::endl << std::endl;
1779 problem_->computeCurrPrecResVec( &*r_ );
1780 index.resize( 1 ); index[0] = 0;
1781 RCP<MV> v00 = MVT::CloneViewNonConst( *V_, index );
1782 MVT::SetBlock(*r_,index,*v00);
1786 index.resize( numBlocks_+1 );
1787 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1788 restartState.
V = MVT::CloneViewNonConst( *V_, index );
1789 index.resize( keff );
1790 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1791 restartState.
U = MVT::CloneViewNonConst( *U_, index );
1792 restartState.
C = MVT::CloneViewNonConst( *C_, index );
1796 gcrodr_iter->initialize(restartState);
1810 true, std::logic_error,
"Belos::GCRODRSolMgr::solve: "
1811 "Invalid return from GCRODRIter::iterate().");
1816 gcrodr_iter->updateLSQR( gcrodr_iter->getCurSubspaceDim() );
1819 sTest_->checkStatus( &*gcrodr_iter );
1820 if (convTest_->getStatus() !=
Passed)
1821 isConverged =
false;
1824 catch (
const std::exception& e) {
1826 <<
"Error! Caught exception in GCRODRIter::iterate() at iteration "
1827 << gcrodr_iter->getNumIters() << std::endl << e.what() << std::endl;
1834 RCP<MV> update = gcrodr_iter->getCurrentUpdate();
1835 problem_->updateSolution( update,
true );
1838 problem_->setCurrLS();
1843 if (!builtRecycleSpace_) {
1844 buildRecycleSpace2(gcrodr_iter);
1845 printer_->stream(
Debug)
1846 <<
" Generated new recycled subspace using RHS index " << currIdx[0]
1847 <<
" of dimension " << keff << std::endl << std::endl;
1852 if (numRHS2Solve > 0) {
1854 problem_->setLSIndex (currIdx);
1857 currIdx.resize (numRHS2Solve);
1865 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1871 #endif // BELOS_TEUCHOS_TIME_MONITOR
1874 numIters_ = maxIterTest_->getNumIters ();
1886 const std::vector<MagnitudeType>* pTestValues = expConvTest_->getTestValue();
1887 if (pTestValues == NULL || pTestValues->size() < 1) {
1888 pTestValues = impConvTest_->getTestValue();
1891 "Belos::GCRODRSolMgr::solve(): The implicit convergence test's getTestValue() "
1892 "method returned NULL. Please report this bug to the Belos developers.");
1894 "Belos::GCRODRSolMgr::solve(): The implicit convergence test's getTestValue() "
1895 "method returned a vector of length zero. Please report this bug to the "
1896 "Belos developers.");
1901 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1908 template<
class ScalarType,
class MV,
class OP>
1914 std::vector<MagnitudeType> d(keff);
1915 std::vector<ScalarType> dscalar(keff);
1916 std::vector<int> index(numBlocks_+1);
1928 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1931 dscalar.resize(keff);
1932 MVT::MvNorm( *Utmp, d );
1933 for (
int i=0; i<keff; ++i) {
1935 dscalar[i] = (ScalarType)d[i];
1937 MVT::MvScale( *Utmp, dscalar );
1944 for (
int i=0; i<keff; ++i) {
1945 (*H2tmp)(i,i) = d[i];
1953 keff_new = getHarmonicVecs2( keff, p, *H2tmp, oldState.
V, PPtmp );
1962 index.resize( keff );
1963 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1965 index.resize( keff_new );
1966 for (
int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
1967 U1tmp = MVT::CloneViewNonConst( *U1_, index );
1969 MVT::MvTimesMatAddMv( one, *Utmp, PPtmp, zero, *U1tmp );
1975 for (
int ii=0; ii < p; ii++) { index[ii] = ii; }
1978 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
1989 int info = 0, lwork = -1;
1990 tau_.resize (keff_new);
1992 HPtmp.
stride (), &tau_[0], &work_[0], lwork, &info);
1995 "LAPACK's _GEQRF failed to compute a workspace size.");
2002 work_.resize (lwork);
2004 HPtmp.
stride (), &tau_[0], &work_[0], lwork, &info);
2007 "LAPACK's _GEQRF failed to compute a QR factorization.");
2012 for(
int i=0;i<keff_new;i++) {
for(
int j=i;j<keff_new;j++) Rtmp(i,j) = HPtmp(i,j); }
2023 "LAPACK's _UNGQR failed to construct the Q factor.");
2033 for (
int i=0; i < keff; i++) { index[i] = i; }
2035 index.resize(keff_new);
2036 for (
int i=0; i < keff_new; i++) { index[i] = i; }
2037 C1tmp = MVT::CloneViewNonConst( *C1_, index );
2039 MVT::MvTimesMatAddMv( one, *Ctmp, PPtmp, zero, *C1tmp );
2043 index.resize( p+1 );
2044 for (
int i=0; i < p+1; ++i) { index[i] = i; }
2047 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *C1tmp );
2062 work_.resize(lwork);
2067 index.resize(keff_new);
2068 for (
int i=0; i < keff_new; i++) { index[i] = i; }
2070 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
2074 if (keff != keff_new) {
2076 gcrodr_iter->setSize( keff, numBlocks_ );
2086 template<
class ScalarType,
class MV,
class OP>
2091 bool xtraVec =
false;
2095 std::vector<MagnitudeType> wr(m), wi(m);
2101 std::vector<MagnitudeType> w(m);
2104 std::vector<int> iperm(m);
2110 builtRecycleSpace_ =
true;
2120 ScalarType d = HH(m, m-1) * HH(m, m-1);
2122 for( i=0; i<m; ++i )
2123 harmHH(i, m-1) += d * e_m[i];
2132 std::vector<ScalarType> work(1);
2133 std::vector<MagnitudeType> rwork(2*m);
2136 lapack.GEEV(
'N',
'V', m, harmHH.
values(), harmHH.
stride(), &wr[0], &wi[0],
2137 vl, ldvl, vr.
values(), vr.
stride(), &work[0], lwork, &rwork[0], &info);
2140 work.resize( lwork );
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.
ScalarType * values() const
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.
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_
Teuchos::RCP< const Teuchos::ParameterList > getDefaultParameters(const std::string &name) const
Default parameters for the given MatOrthoManager subclass.
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_
static const double convTol
Default convergence tolerance.
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.
std::string validNamesString() const
List (as a string) of recognized MatOrthoManager names.
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
bool isParameter(const std::string &name) const
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. ...
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)
bool isValidName(const std::string &name) const
Whether this factory recognizes the MatOrthoManager with the given name.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > convTest_
Teuchos::RCP< MV > U
The recycled subspace and its projection.
std::ostream & printValidNames(std::ostream &out) const
Print all recognized MatOrthoManager names to the given ostream.
std::vector< ScalarType > work_
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
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...
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.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
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.
OrdinalType numCols() const
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
bool isType(const std::string &name) const
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.
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
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.
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::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > R_
OrdinalType stride() const
OrdinalType numRows() const
Belos concrete class for performing the block, flexible GMRES iteration.
Teuchos::RCP< Belos::MatOrthoManager< Scalar, MV, OP > > makeMatOrthoManager(const std::string &ortho, const Teuchos::RCP< const OP > &M, const Teuchos::RCP< OutputManager< Scalar > > &, const std::string &label, const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Return an instance of the specified MatOrthoManager subclass.