42 #ifndef BELOS_GCRODR_SOLMGR_HPP 43 #define BELOS_GCRODR_SOLMGR_HPP 62 #include "Teuchos_BLAS.hpp" 63 #include "Teuchos_LAPACK.hpp" 64 #include "Teuchos_as.hpp" 65 #ifdef BELOS_TEUCHOS_TIME_MONITOR 66 # include "Teuchos_TimeMonitor.hpp" 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 =
160 static const bool requiresLapack =
170 const Teuchos::RCP<Teuchos::ParameterList>& pl) :
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, float>::value ||
190 std::is_same<ScalarType, double>::value ||
191 std::is_same<ScalarType, long double>::value,
192 "Belos::GCRODRSolMgr: ScalarType must be one of the four " 193 "types (S,D,C,Z) supported by LAPACK or long double (largely not impl'd).");
195 static_assert (std::is_same<ScalarType, std::complex<float> >::value ||
196 std::is_same<ScalarType, std::complex<double> >::value ||
197 std::is_same<ScalarType, float>::value ||
198 std::is_same<ScalarType, double>::value,
199 "Belos::GCRODRSolMgr: ScalarType must be one of the four " 200 "types (S,D,C,Z) supported by LAPACK.");
203 #if defined(HAVE_TEUCHOS_LONG_DOUBLE) 204 static_assert (std::is_same<ScalarType, float>::value ||
205 std::is_same<ScalarType, double>::value ||
206 std::is_same<ScalarType, long double>::value,
207 "Belos::GCRODRSolMgr: ScalarType must be float, double or long double. " 208 "Complex arithmetic support is currently disabled. To " 209 "enable it, set Teuchos_ENABLE_COMPLEX=ON.");
211 static_assert (std::is_same<ScalarType, float>::value ||
212 std::is_same<ScalarType, double>::value,
213 "Belos::GCRODRSolMgr: ScalarType must be float or double. " 214 "Complex arithmetic support is currently disabled. To " 215 "enable it, set Teuchos_ENABLE_COMPLEX=ON.");
217 # endif // defined(HAVE_TEUCHOS_COMPLEX) 218 #endif // defined(HAVE_TEUCHOSCORE_CXX11) 223 typedef Teuchos::ScalarTraits<ScalarType> SCT;
224 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
225 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
292 const Teuchos::RCP<Teuchos::ParameterList> &pl);
298 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
314 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters()
const override;
327 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
328 return Teuchos::tuple(timerSolve_);
360 void setParameters(
const Teuchos::RCP<Teuchos::ParameterList> ¶ms )
override;
372 bool set = problem_->setProblem();
374 throw "Could not set problem.";
418 std::string description()
const override;
428 void initializeStateStorage();
437 int getHarmonicVecs1(
int m,
438 const Teuchos::SerialDenseMatrix<int,ScalarType>& HH,
439 Teuchos::SerialDenseMatrix<int,ScalarType>& PP);
446 int getHarmonicVecs2(
int keff,
int m,
447 const Teuchos::SerialDenseMatrix<int,ScalarType>& HH,
448 const Teuchos::RCP<const MV>& VV,
449 Teuchos::SerialDenseMatrix<int,ScalarType>& PP);
452 void sort(std::vector<MagnitudeType>& dlist,
int n, std::vector<int>& iperm);
455 Teuchos::LAPACK<int,ScalarType> lapack;
458 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
461 Teuchos::RCP<OutputManager<ScalarType> > printer_;
462 Teuchos::RCP<std::ostream> outputStream_;
465 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
466 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
467 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
468 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
469 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
474 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
477 Teuchos::RCP<Teuchos::ParameterList> params_;
480 static constexpr
double orthoKappa_default_ = 0.0;
481 static constexpr
int maxRestarts_default_ = 100;
482 static constexpr
int maxIters_default_ = 1000;
483 static constexpr
int numBlocks_default_ = 50;
484 static constexpr
int blockSize_default_ = 1;
485 static constexpr
int recycledBlocks_default_ = 5;
488 static constexpr
int outputFreq_default_ = -1;
489 static constexpr
const char * impResScale_default_ =
"Norm of Preconditioned Initial Residual";
490 static constexpr
const char * expResScale_default_ =
"Norm of Initial Residual";
491 static constexpr
const char * label_default_ =
"Belos";
492 static constexpr
const char * orthoType_default_ =
"ICGS";
494 #if defined(_WIN32) && defined(__clang__) 495 static constexpr std::ostream * outputStream_default_ =
496 __builtin_constant_p(reinterpret_cast<const std::ostream*>(&std::cout));
498 static constexpr std::ostream * outputStream_default_ = &std::cout;
502 MagnitudeType convTol_, orthoKappa_, achievedTol_;
503 int maxRestarts_, maxIters_, numIters_;
504 int verbosity_, outputStyle_, outputFreq_;
505 std::string orthoType_;
506 std::string impResScale_, expResScale_;
513 int numBlocks_, recycledBlocks_;
524 Teuchos::RCP<MV> U_, C_;
527 Teuchos::RCP<MV> U1_, C1_;
530 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H2_;
531 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_;
532 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B_;
533 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PP_;
534 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > HP_;
535 std::vector<ScalarType> tau_;
536 std::vector<ScalarType> work_;
537 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > R_;
538 std::vector<int> ipiv_;
543 Teuchos::RCP<Teuchos::Time> timerSolve_;
549 bool builtRecycleSpace_;
554 template<
class ScalarType,
class MV,
class OP>
564 template<
class ScalarType,
class MV,
class OP>
567 const Teuchos::RCP<Teuchos::ParameterList>& pl):
575 TEUCHOS_TEST_FOR_EXCEPTION(
576 problem == Teuchos::null, std::invalid_argument,
577 "Belos::GCRODRSolMgr constructor: The solver manager's " 578 "constructor needs the linear problem argument 'problem' " 587 if (! pl.is_null ()) {
593 template<
class ScalarType,
class MV,
class OP>
595 outputStream_ = Teuchos::rcp(outputStream_default_,
false);
597 orthoKappa_ = orthoKappa_default_;
598 maxRestarts_ = maxRestarts_default_;
599 maxIters_ = maxIters_default_;
600 numBlocks_ = numBlocks_default_;
601 recycledBlocks_ = recycledBlocks_default_;
602 verbosity_ = verbosity_default_;
603 outputStyle_ = outputStyle_default_;
604 outputFreq_ = outputFreq_default_;
605 orthoType_ = orthoType_default_;
606 impResScale_ = impResScale_default_;
607 expResScale_ = expResScale_default_;
608 label_ = label_default_;
610 builtRecycleSpace_ =
false;
626 template<
class ScalarType,
class MV,
class OP>
631 using Teuchos::isParameterType;
632 using Teuchos::getParameter;
634 using Teuchos::ParameterList;
635 using Teuchos::parameterList;
638 using Teuchos::rcp_dynamic_cast;
639 using Teuchos::rcpFromRef;
640 using Teuchos::Exceptions::InvalidParameter;
641 using Teuchos::Exceptions::InvalidParameterName;
642 using Teuchos::Exceptions::InvalidParameterType;
646 RCP<const ParameterList> defaultParams = getValidParameters();
664 if (params_.is_null()) {
665 params_ = parameterList (*defaultParams);
673 if (params_ != params) {
679 params_ = parameterList (*params);
714 params_->validateParametersAndSetDefaults (*defaultParams);
718 if (params->isParameter (
"Maximum Restarts")) {
719 maxRestarts_ = params->get(
"Maximum Restarts", maxRestarts_default_);
722 params_->set (
"Maximum Restarts", maxRestarts_);
726 if (params->isParameter (
"Maximum Iterations")) {
727 maxIters_ = params->get (
"Maximum Iterations", maxIters_default_);
730 params_->set (
"Maximum Iterations", maxIters_);
731 if (! maxIterTest_.is_null())
732 maxIterTest_->setMaxIters (maxIters_);
736 if (params->isParameter (
"Num Blocks")) {
737 numBlocks_ = params->get (
"Num Blocks", numBlocks_default_);
738 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
739 "Belos::GCRODRSolMgr: The \"Num Blocks\" parameter must " 740 "be strictly positive, but you specified a value of " 741 << numBlocks_ <<
".");
743 params_->set (
"Num Blocks", numBlocks_);
747 if (params->isParameter (
"Num Recycled Blocks")) {
748 recycledBlocks_ = params->get (
"Num Recycled Blocks",
749 recycledBlocks_default_);
750 TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks_ <= 0, std::invalid_argument,
751 "Belos::GCRODRSolMgr: The \"Num Recycled Blocks\" " 752 "parameter must be strictly positive, but you specified " 753 "a value of " << recycledBlocks_ <<
".");
754 TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks_ >= numBlocks_, std::invalid_argument,
755 "Belos::GCRODRSolMgr: The \"Num Recycled Blocks\" " 756 "parameter must be less than the \"Num Blocks\" " 757 "parameter, but you specified \"Num Recycled Blocks\" " 758 "= " << recycledBlocks_ <<
" and \"Num Blocks\" = " 759 << numBlocks_ <<
".");
761 params_->set(
"Num Recycled Blocks", recycledBlocks_);
767 if (params->isParameter (
"Timer Label")) {
768 std::string tempLabel = params->get (
"Timer Label", label_default_);
771 if (tempLabel != label_) {
773 params_->set (
"Timer Label", label_);
774 std::string solveLabel = label_ +
": GCRODRSolMgr total solve time";
775 #ifdef BELOS_TEUCHOS_TIME_MONITOR 776 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
778 if (ortho_ != Teuchos::null) {
779 ortho_->setLabel( label_ );
785 if (params->isParameter (
"Verbosity")) {
786 if (isParameterType<int> (*params,
"Verbosity")) {
787 verbosity_ = params->get (
"Verbosity", verbosity_default_);
789 verbosity_ = (int) getParameter<Belos::MsgType> (*params,
"Verbosity");
792 params_->set (
"Verbosity", verbosity_);
795 if (! printer_.is_null())
796 printer_->setVerbosity (verbosity_);
800 if (params->isParameter (
"Output Style")) {
801 if (isParameterType<int> (*params,
"Output Style")) {
802 outputStyle_ = params->get (
"Output Style", outputStyle_default_);
804 outputStyle_ = (int) getParameter<OutputType> (*params,
"Output Style");
808 params_->set (
"Output Style", outputStyle_);
826 if (params->isParameter (
"Output Stream")) {
828 outputStream_ = getParameter<RCP<std::ostream> > (*params,
"Output Stream");
829 }
catch (InvalidParameter&) {
830 outputStream_ = rcpFromRef (std::cout);
837 if (outputStream_.is_null()) {
838 outputStream_ = rcp (
new Teuchos::oblackholestream);
841 params_->set (
"Output Stream", outputStream_);
844 if (! printer_.is_null()) {
845 printer_->setOStream (outputStream_);
851 if (params->isParameter (
"Output Frequency")) {
852 outputFreq_ = params->get (
"Output Frequency", outputFreq_default_);
856 params_->set(
"Output Frequency", outputFreq_);
857 if (! outputTest_.is_null())
858 outputTest_->setOutputFrequency (outputFreq_);
865 if (printer_.is_null()) {
876 bool changedOrthoType =
false;
877 if (params->isParameter (
"Orthogonalization")) {
878 const std::string& tempOrthoType =
879 params->get (
"Orthogonalization", orthoType_default_);
882 std::ostringstream os;
883 os <<
"Belos::GCRODRSolMgr: Invalid orthogonalization name \"" 884 << tempOrthoType <<
"\". The following are valid options " 885 <<
"for the \"Orthogonalization\" name parameter: ";
887 throw std::invalid_argument (os.str());
889 if (tempOrthoType != orthoType_) {
890 changedOrthoType =
true;
891 orthoType_ = tempOrthoType;
893 params_->set (
"Orthogonalization", orthoType_);
909 RCP<ParameterList> orthoParams;
912 using Teuchos::sublist;
914 const std::string paramName (
"Orthogonalization Parameters");
917 orthoParams = sublist (params_, paramName,
true);
918 }
catch (InvalidParameter&) {
925 orthoParams = sublist (params_, paramName,
true);
928 TEUCHOS_TEST_FOR_EXCEPTION(orthoParams.is_null(), std::logic_error,
929 "Failed to get orthogonalization parameters. " 930 "Please report this bug to the Belos developers.");
935 if (ortho_.is_null() || changedOrthoType) {
941 label_, orthoParams);
949 typedef Teuchos::ParameterListAcceptor PLA;
950 RCP<PLA> pla = rcp_dynamic_cast<PLA> (ortho_);
956 label_, orthoParams);
958 pla->setParameterList (orthoParams);
970 if (params->isParameter (
"Orthogonalization Constant")) {
971 MagnitudeType orthoKappa = orthoKappa_default_;
972 if (params->isType<MagnitudeType> (
"Orthogonalization Constant")) {
973 orthoKappa = params->get (
"Orthogonalization Constant", orthoKappa);
976 orthoKappa = params->get (
"Orthogonalization Constant", orthoKappa_default_);
979 if (orthoKappa > 0) {
980 orthoKappa_ = orthoKappa;
982 params_->set(
"Orthogonalization Constant", orthoKappa_);
984 if (orthoType_ ==
"DGKS" && ! ortho_.is_null()) {
991 rcp_dynamic_cast<ortho_man_type>(ortho_)->setDepTol (orthoKappa_);
1001 if (params->isParameter(
"Convergence Tolerance")) {
1002 if (params->isType<MagnitudeType> (
"Convergence Tolerance")) {
1003 convTol_ = params->get (
"Convergence Tolerance",
1011 params_->set (
"Convergence Tolerance", convTol_);
1012 if (! impConvTest_.is_null())
1014 if (! expConvTest_.is_null())
1015 expConvTest_->setTolerance (convTol_);
1019 if (params->isParameter (
"Implicit Residual Scaling")) {
1020 std::string tempImpResScale =
1021 getParameter<std::string> (*params,
"Implicit Residual Scaling");
1024 if (impResScale_ != tempImpResScale) {
1026 impResScale_ = tempImpResScale;
1029 params_->set(
"Implicit Residual Scaling", impResScale_);
1039 if (! impConvTest_.is_null()) {
1045 impConvTest_ = null;
1052 if (params->isParameter(
"Explicit Residual Scaling")) {
1053 std::string tempExpResScale =
1054 getParameter<std::string> (*params,
"Explicit Residual Scaling");
1057 if (expResScale_ != tempExpResScale) {
1059 expResScale_ = tempExpResScale;
1062 params_->set(
"Explicit Residual Scaling", expResScale_);
1065 if (! expConvTest_.is_null()) {
1071 expConvTest_ = null;
1082 if (maxIterTest_.is_null())
1087 if (impConvTest_.is_null()) {
1088 impConvTest_ = rcp (
new StatusTestResNorm_t (convTol_));
1094 if (expConvTest_.is_null()) {
1095 expConvTest_ = rcp (
new StatusTestResNorm_t (convTol_));
1096 expConvTest_->defineResForm (StatusTestResNorm_t::Explicit,
Belos::TwoNorm);
1102 if (convTest_.is_null()) {
1103 convTest_ = rcp (
new StatusTestCombo_t (StatusTestCombo_t::SEQ,
1111 sTest_ = rcp (
new StatusTestCombo_t (StatusTestCombo_t::OR,
1117 outputTest_ = stoFactory.
create (printer_, sTest_, outputFreq_,
1121 std::string solverDesc =
" GCRODR ";
1122 outputTest_->setSolverDesc( solverDesc );
1125 if (timerSolve_.is_null()) {
1126 std::string solveLabel = label_ +
": GCRODRSolMgr total solve time";
1127 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1128 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
1137 template<
class ScalarType,
class MV,
class OP>
1138 Teuchos::RCP<const Teuchos::ParameterList>
1141 using Teuchos::ParameterList;
1142 using Teuchos::parameterList;
1145 static RCP<const ParameterList> validPL;
1146 if (is_null(validPL)) {
1147 RCP<ParameterList> pl = parameterList ();
1151 "The relative residual tolerance that needs to be achieved by the\n" 1152 "iterative solver in order for the linear system to be declared converged.");
1153 pl->set(
"Maximum Restarts", static_cast<int>(maxRestarts_default_),
1154 "The maximum number of cycles allowed for each\n" 1155 "set of RHS solved.");
1156 pl->set(
"Maximum Iterations", static_cast<int>(maxIters_default_),
1157 "The maximum number of iterations allowed for each\n" 1158 "set of RHS solved.");
1162 pl->set(
"Block Size", static_cast<int>(blockSize_default_),
1163 "Block Size Parameter -- currently must be 1 for GCRODR");
1164 pl->set(
"Num Blocks", static_cast<int>(numBlocks_default_),
1165 "The maximum number of vectors allowed in the Krylov subspace\n" 1166 "for each set of RHS solved.");
1167 pl->set(
"Num Recycled Blocks", static_cast<int>(recycledBlocks_default_),
1168 "The maximum number of vectors in the recycled subspace." );
1169 pl->set(
"Verbosity", static_cast<int>(verbosity_default_),
1170 "What type(s) of solver information should be outputted\n" 1171 "to the output stream.");
1172 pl->set(
"Output Style", static_cast<int>(outputStyle_default_),
1173 "What style is used for the solver information outputted\n" 1174 "to the output stream.");
1175 pl->set(
"Output Frequency", static_cast<int>(outputFreq_default_),
1176 "How often convergence information should be outputted\n" 1177 "to the output stream.");
1178 pl->set(
"Output Stream", Teuchos::rcp(outputStream_default_,
false),
1179 "A reference-counted pointer to the output stream where all\n" 1180 "solver output is sent.");
1181 pl->set(
"Implicit Residual Scaling", static_cast<const char *>(impResScale_default_),
1182 "The type of scaling used in the implicit residual convergence test.");
1183 pl->set(
"Explicit Residual Scaling", static_cast<const char *>(expResScale_default_),
1184 "The type of scaling used in the explicit residual convergence test.");
1185 pl->set(
"Timer Label", static_cast<const char *>(label_default_),
1186 "The string to use as a prefix for the timer labels.");
1189 pl->set(
"Orthogonalization", static_cast<const char *>(orthoType_default_),
1190 "The type of orthogonalization to use. Valid options: " +
1192 RCP<const ParameterList> orthoParams =
1194 pl->set (
"Orthogonalization Parameters", *orthoParams,
1195 "Parameters specific to the type of orthogonalization used.");
1197 pl->set(
"Orthogonalization Constant",static_cast<MagnitudeType>(orthoKappa_default_),
1198 "When using DGKS orthogonalization: the \"depTol\" constant, used " 1199 "to determine whether another step of classical Gram-Schmidt is " 1200 "necessary. Otherwise ignored.");
1207 template<
class ScalarType,
class MV,
class OP>
1210 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1213 Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
1214 if (rhsMV == Teuchos::null) {
1221 TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(numBlocks_) > MVT::GetGlobalLength(*rhsMV),std::invalid_argument,
1222 "Belos::GCRODRSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
1225 if (U_ == Teuchos::null) {
1226 U_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1230 if (MVT::GetNumberVecs(*U_) < recycledBlocks_+1) {
1231 Teuchos::RCP<const MV> tmp = U_;
1232 U_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1237 if (C_ == Teuchos::null) {
1238 C_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1242 if (MVT::GetNumberVecs(*C_) < recycledBlocks_+1) {
1243 Teuchos::RCP<const MV> tmp = C_;
1244 C_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1249 if (V_ == Teuchos::null) {
1250 V_ = MVT::Clone( *rhsMV, numBlocks_+1 );
1254 if (MVT::GetNumberVecs(*V_) < numBlocks_+1) {
1255 Teuchos::RCP<const MV> tmp = V_;
1256 V_ = MVT::Clone( *tmp, numBlocks_+1 );
1261 if (U1_ == Teuchos::null) {
1262 U1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1266 if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1267 Teuchos::RCP<const MV> tmp = U1_;
1268 U1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1273 if (C1_ == Teuchos::null) {
1274 C1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1278 if (MVT::GetNumberVecs(*C1_) < recycledBlocks_+1) {
1279 Teuchos::RCP<const MV> tmp = C1_;
1280 C1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1285 if (r_ == Teuchos::null)
1286 r_ = MVT::Clone( *rhsMV, 1 );
1289 tau_.resize(recycledBlocks_+1);
1292 work_.resize(recycledBlocks_+1);
1295 ipiv_.resize(recycledBlocks_+1);
1298 if (H2_ == Teuchos::null)
1299 H2_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 ) );
1301 if ( (H2_->numRows() != numBlocks_+recycledBlocks_+2) || (H2_->numCols() != numBlocks_+recycledBlocks_+1) )
1302 H2_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
1304 H2_->putScalar(zero);
1307 if (R_ == Teuchos::null)
1308 R_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycledBlocks_+1, recycledBlocks_+1 ) );
1310 if ( (R_->numRows() != recycledBlocks_+1) || (R_->numCols() != recycledBlocks_+1) )
1311 R_->reshape( recycledBlocks_+1, recycledBlocks_+1 );
1313 R_->putScalar(zero);
1316 if (PP_ == Teuchos::null)
1317 PP_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, recycledBlocks_+1 ) );
1319 if ( (PP_->numRows() != numBlocks_+recycledBlocks_+2) || (PP_->numCols() != recycledBlocks_+1) )
1320 PP_->reshape( numBlocks_+recycledBlocks_+2, recycledBlocks_+1 );
1324 if (HP_ == Teuchos::null)
1325 HP_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 ) );
1327 if ( (HP_->numRows() != numBlocks_+recycledBlocks_+2) || (HP_->numCols() != numBlocks_+recycledBlocks_+1) )
1328 HP_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
1336 template<
class ScalarType,
class MV,
class OP>
1344 if (!isSet_) { setParameters( params_ ); }
1346 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1347 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1348 std::vector<int> index(numBlocks_+1);
1350 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,
GCRODRSolMgrLinearProblemFailure,
"Belos::GCRODRSolMgr::solve(): Linear problem is not a valid object.");
1352 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),
GCRODRSolMgrLinearProblemFailure,
"Belos::GCRODRSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
1355 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1356 std::vector<int> currIdx(1);
1360 problem_->setLSIndex( currIdx );
1363 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1364 if (static_cast<ptrdiff_t>(numBlocks_) > dim) {
1365 numBlocks_ = Teuchos::as<int>(dim);
1367 "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
1368 " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
1369 params_->set(
"Num Blocks", numBlocks_);
1373 bool isConverged =
true;
1376 initializeStateStorage();
1380 Teuchos::ParameterList plist;
1382 plist.set(
"Num Blocks",numBlocks_);
1383 plist.set(
"Recycled Blocks",recycledBlocks_);
1388 RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_iter;
1391 int prime_iterations = 0;
1395 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1396 Teuchos::TimeMonitor slvtimer(*timerSolve_);
1399 while ( numRHS2Solve > 0 ) {
1402 builtRecycleSpace_ =
false;
1405 outputTest_->reset();
1413 "Belos::GCRODRSolMgr::solve(): Requested size of recycled subspace is not consistent with the current recycle subspace.");
1415 printer_->stream(
Debug) <<
" Now solving RHS index " << currIdx[0] <<
" using recycled subspace of dimension " << keff << std::endl << std::endl;
1418 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1419 RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1420 RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1421 problem_->apply( *Utmp, *Ctmp );
1423 RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1427 Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff, keff );
1428 int rank = ortho_->normalize(*Ctmp, rcp(&Rtmp,
false));
1430 TEUCHOS_TEST_FOR_EXCEPTION(rank != keff,
GCRODRSolMgrOrthoFailure,
"Belos::GCRODRSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
1435 ipiv_.resize(Rtmp.numRows());
1436 lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
1437 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1440 int lwork = Rtmp.numRows();
1441 work_.resize(lwork);
1442 lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1443 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
1446 MVT::MvTimesMatAddMv( one, *Utmp, Rtmp, zero, *U1tmp );
1451 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1452 Ctmp = MVT::CloneViewNonConst( *C_, index );
1453 Utmp = MVT::CloneView( *U_, index );
1456 Teuchos::SerialDenseMatrix<int,ScalarType> Ctr(keff,1);
1457 problem_->computeCurrPrecResVec( &*r_ );
1458 MVT::MvTransMv( one, *Ctmp, *r_, Ctr );
1461 RCP<MV> update = MVT::Clone( *problem_->getCurrLHSVec(), 1 );
1462 MVT::MvInit( *update, 0.0 );
1463 MVT::MvTimesMatAddMv( one, *Utmp, Ctr, one, *update );
1464 problem_->updateSolution( update,
true );
1467 MVT::MvTimesMatAddMv( -one, *Ctmp, Ctr, one, *r_ );
1470 prime_iterations = 0;
1476 printer_->stream(
Debug) <<
" No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
1478 Teuchos::ParameterList primeList;
1481 primeList.set(
"Num Blocks",numBlocks_);
1482 primeList.set(
"Recycled Blocks",0);
1485 RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_prime_iter;
1489 problem_->computeCurrPrecResVec( &*r_ );
1490 index.resize( 1 ); index[0] = 0;
1491 RCP<MV> v0 = MVT::CloneViewNonConst( *V_, index );
1492 MVT::SetBlock(*r_,index,*v0);
1496 index.resize( numBlocks_+1 );
1497 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1498 newstate.
V = MVT::CloneViewNonConst( *V_, index );
1499 newstate.
U = Teuchos::null;
1500 newstate.
C = Teuchos::null;
1501 newstate.
H = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, recycledBlocks_+1, recycledBlocks_+1 ) );
1502 newstate.
B = Teuchos::null;
1504 gcrodr_prime_iter->initialize(newstate);
1507 bool primeConverged =
false;
1509 gcrodr_prime_iter->iterate();
1512 if ( convTest_->getStatus() ==
Passed ) {
1514 primeConverged =
true;
1519 gcrodr_prime_iter->updateLSQR( gcrodr_prime_iter->getCurSubspaceDim() );
1522 sTest_->checkStatus( &*gcrodr_prime_iter );
1523 if (convTest_->getStatus() ==
Passed)
1524 primeConverged =
true;
1526 catch (
const std::exception &e) {
1527 printer_->stream(
Errors) <<
"Error! Caught exception in GCRODRIter::iterate() at iteration " 1528 << gcrodr_prime_iter->getNumIters() << std::endl
1529 << e.what() << std::endl;
1533 prime_iterations = gcrodr_prime_iter->getNumIters();
1536 RCP<MV> update = gcrodr_prime_iter->getCurrentUpdate();
1537 problem_->updateSolution( update,
true );
1540 newstate = gcrodr_prime_iter->getState();
1548 if (recycledBlocks_ < p+1) {
1550 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PPtmp = rcp (
new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, p, recycledBlocks_+1 ) );
1552 keff = getHarmonicVecs1( p, *newstate.
H, *PPtmp );
1554 PPtmp = rcp (
new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, p, keff ) );
1557 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1558 RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1559 RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1560 RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1562 for (
int ii=0; ii < p; ++ii) { index[ii] = ii; }
1563 RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1567 MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
1574 Teuchos::SerialDenseMatrix<int,ScalarType> Htmp( Teuchos::View, *H2_, p+1, p, recycledBlocks_+1,recycledBlocks_+1);
1575 Teuchos::SerialDenseMatrix<int,ScalarType> HPtmp( Teuchos::View, *HP_, p+1, keff );
1576 HPtmp.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, Htmp, *PPtmp, zero );
1583 lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
1584 HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
1585 TEUCHOS_TEST_FOR_EXCEPTION(
1587 " LAPACK's _GEQRF failed to compute a workspace size.");
1595 lwork = std::abs (static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work_[0])));
1596 work_.resize (lwork);
1597 lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
1598 HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
1599 TEUCHOS_TEST_FOR_EXCEPTION(
1601 " LAPACK's _GEQRF failed to compute a QR factorization.");
1605 Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff, keff );
1606 for (
int ii = 0; ii < keff; ++ii) {
1607 for (
int jj = ii; jj < keff; ++jj) {
1608 Rtmp(ii,jj) = HPtmp(ii,jj);
1615 lapack.UNGQR (HPtmp.numRows (), HPtmp.numCols (), HPtmp.numCols (),
1616 HPtmp.values (), HPtmp.stride (), &tau_[0], &work_[0],
1618 TEUCHOS_TEST_FOR_EXCEPTION(
1620 "LAPACK's _UNGQR failed to construct the Q factor.");
1625 index.resize (p + 1);
1626 for (
int ii = 0; ii < (p+1); ++ii) {
1629 Vtmp = MVT::CloneView( *V_, index );
1630 MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
1637 ipiv_.resize(Rtmp.numRows());
1638 lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
1639 TEUCHOS_TEST_FOR_EXCEPTION(
1641 "LAPACK's _GETRF failed to compute an LU factorization.");
1650 lwork = Rtmp.numRows();
1651 work_.resize(lwork);
1652 lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1653 TEUCHOS_TEST_FOR_EXCEPTION(
1655 "LAPACK's _GETRI failed to invert triangular matrix.");
1658 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
1660 printer_->stream(
Debug)
1661 <<
" Generated recycled subspace using RHS index " << currIdx[0]
1662 <<
" of dimension " << keff << std::endl << std::endl;
1667 if (primeConverged) {
1669 problem_->setCurrLS();
1673 if (numRHS2Solve > 0) {
1675 problem_->setLSIndex (currIdx);
1678 currIdx.resize (numRHS2Solve);
1688 gcrodr_iter->setSize( keff, numBlocks_ );
1691 gcrodr_iter->resetNumIters(prime_iterations);
1694 outputTest_->resetNumCalls();
1697 problem_->computeCurrPrecResVec( &*r_ );
1698 index.resize( 1 ); index[0] = 0;
1699 RCP<MV> v0 = MVT::CloneViewNonConst( *V_, index );
1700 MVT::SetBlock(*r_,index,*v0);
1704 index.resize( numBlocks_+1 );
1705 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1706 newstate.
V = MVT::CloneViewNonConst( *V_, index );
1707 index.resize( keff );
1708 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1709 newstate.
C = MVT::CloneViewNonConst( *C_, index );
1710 newstate.
U = MVT::CloneViewNonConst( *U_, index );
1711 newstate.
B = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, keff, numBlocks_, 0, keff ) );
1712 newstate.
H = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, keff, keff ) );
1714 gcrodr_iter->initialize(newstate);
1717 int numRestarts = 0;
1722 gcrodr_iter->iterate();
1729 if ( convTest_->getStatus() ==
Passed ) {
1738 else if ( maxIterTest_->getStatus() ==
Passed ) {
1740 isConverged =
false;
1748 else if ( gcrodr_iter->getCurSubspaceDim() == gcrodr_iter->getMaxSubspaceDim() ) {
1753 RCP<MV> update = gcrodr_iter->getCurrentUpdate();
1754 problem_->updateSolution( update,
true );
1756 buildRecycleSpace2(gcrodr_iter);
1758 printer_->stream(
Debug)
1759 <<
" Generated new recycled subspace using RHS index " 1760 << currIdx[0] <<
" of dimension " << keff << std::endl
1764 if (numRestarts >= maxRestarts_) {
1765 isConverged =
false;
1770 printer_->stream(
Debug)
1771 <<
" Performing restart number " << numRestarts <<
" of " 1772 << maxRestarts_ << std::endl << std::endl;
1775 problem_->computeCurrPrecResVec( &*r_ );
1776 index.resize( 1 ); index[0] = 0;
1777 RCP<MV> v00 = MVT::CloneViewNonConst( *V_, index );
1778 MVT::SetBlock(*r_,index,*v00);
1782 index.resize( numBlocks_+1 );
1783 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1784 restartState.
V = MVT::CloneViewNonConst( *V_, index );
1785 index.resize( keff );
1786 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1787 restartState.
U = MVT::CloneViewNonConst( *U_, index );
1788 restartState.
C = MVT::CloneViewNonConst( *C_, index );
1789 restartState.
B = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, keff, numBlocks_, 0, keff ) );
1790 restartState.
H = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, keff, keff ) );
1792 gcrodr_iter->initialize(restartState);
1805 TEUCHOS_TEST_FOR_EXCEPTION(
1806 true, std::logic_error,
"Belos::GCRODRSolMgr::solve: " 1807 "Invalid return from GCRODRIter::iterate().");
1812 gcrodr_iter->updateLSQR( gcrodr_iter->getCurSubspaceDim() );
1815 sTest_->checkStatus( &*gcrodr_iter );
1816 if (convTest_->getStatus() !=
Passed)
1817 isConverged =
false;
1820 catch (
const std::exception& e) {
1822 <<
"Error! Caught exception in GCRODRIter::iterate() at iteration " 1823 << gcrodr_iter->getNumIters() << std::endl << e.what() << std::endl;
1830 RCP<MV> update = gcrodr_iter->getCurrentUpdate();
1831 problem_->updateSolution( update,
true );
1834 problem_->setCurrLS();
1839 if (!builtRecycleSpace_) {
1840 buildRecycleSpace2(gcrodr_iter);
1841 printer_->stream(
Debug)
1842 <<
" Generated new recycled subspace using RHS index " << currIdx[0]
1843 <<
" of dimension " << keff << std::endl << std::endl;
1848 if (numRHS2Solve > 0) {
1850 problem_->setLSIndex (currIdx);
1853 currIdx.resize (numRHS2Solve);
1861 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1866 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails) );
1867 #endif // BELOS_TEUCHOS_TIME_MONITOR 1870 numIters_ = maxIterTest_->getNumIters ();
1882 const std::vector<MagnitudeType>* pTestValues = expConvTest_->getTestValue();
1883 if (pTestValues == NULL || pTestValues->size() < 1) {
1884 pTestValues = impConvTest_->getTestValue();
1886 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1887 "Belos::GCRODRSolMgr::solve(): The implicit convergence test's getTestValue() " 1888 "method returned NULL. Please report this bug to the Belos developers.");
1889 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1890 "Belos::GCRODRSolMgr::solve(): The implicit convergence test's getTestValue() " 1891 "method returned a vector of length zero. Please report this bug to the " 1892 "Belos developers.");
1897 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1904 template<
class ScalarType,
class MV,
class OP>
1907 MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
1908 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1910 std::vector<MagnitudeType> d(keff);
1911 std::vector<ScalarType> dscalar(keff);
1912 std::vector<int> index(numBlocks_+1);
1924 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1925 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1927 dscalar.resize(keff);
1928 MVT::MvNorm( *Utmp, d );
1929 for (
int i=0; i<keff; ++i) {
1931 dscalar[i] = (ScalarType)d[i];
1933 MVT::MvScale( *Utmp, dscalar );
1937 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H2tmp = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, p+keff+1, p+keff ) );
1940 for (
int i=0; i<keff; ++i) {
1941 (*H2tmp)(i,i) = d[i];
1948 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p+keff, recycledBlocks_+1 );
1949 keff_new = getHarmonicVecs2( keff, p, *H2tmp, oldState.
V, PPtmp );
1956 Teuchos::RCP<MV> U1tmp;
1958 index.resize( keff );
1959 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1960 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1961 index.resize( keff_new );
1962 for (
int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
1963 U1tmp = MVT::CloneViewNonConst( *U1_, index );
1964 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, keff, keff_new );
1965 MVT::MvTimesMatAddMv( one, *Utmp, PPtmp, zero, *U1tmp );
1971 for (
int ii=0; ii < p; ii++) { index[ii] = ii; }
1972 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1973 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p, keff_new, keff );
1974 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
1978 Teuchos::SerialDenseMatrix<int,ScalarType> HPtmp( Teuchos::View, *HP_, p+keff+1, keff_new );
1980 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p+keff, keff_new );
1981 HPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*H2tmp,PPtmp,zero);
1985 int info = 0, lwork = -1;
1986 tau_.resize (keff_new);
1987 lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
1988 HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
1989 TEUCHOS_TEST_FOR_EXCEPTION(
1990 info != 0, GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve: " 1991 "LAPACK's _GEQRF failed to compute a workspace size.");
1997 lwork = std::abs (static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work_[0])));
1998 work_.resize (lwork);
1999 lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
2000 HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
2001 TEUCHOS_TEST_FOR_EXCEPTION(
2002 info != 0, GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve: " 2003 "LAPACK's _GEQRF failed to compute a QR factorization.");
2007 Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff_new, keff_new );
2008 for(
int i=0;i<keff_new;i++) {
for(
int j=i;j<keff_new;j++) Rtmp(i,j) = HPtmp(i,j); }
2014 lapack.UNGQR (HPtmp.numRows (), HPtmp.numCols (), HPtmp.numCols (),
2015 HPtmp.values (), HPtmp.stride (), &tau_[0], &work_[0],
2017 TEUCHOS_TEST_FOR_EXCEPTION(
2018 info != 0, GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve: " 2019 "LAPACK's _UNGQR failed to construct the Q factor.");
2026 Teuchos::RCP<MV> C1tmp;
2029 for (
int i=0; i < keff; i++) { index[i] = i; }
2030 Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
2031 index.resize(keff_new);
2032 for (
int i=0; i < keff_new; i++) { index[i] = i; }
2033 C1tmp = MVT::CloneViewNonConst( *C1_, index );
2034 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *HP_, keff, keff_new );
2035 MVT::MvTimesMatAddMv( one, *Ctmp, PPtmp, zero, *C1tmp );
2039 index.resize( p+1 );
2040 for (
int i=0; i < p+1; ++i) { index[i] = i; }
2041 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
2042 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *HP_, p+1, keff_new, keff, 0 );
2043 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *C1tmp );
2052 ipiv_.resize(Rtmp.numRows());
2053 lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
2054 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
2057 lwork = Rtmp.numRows();
2058 work_.resize(lwork);
2059 lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
2060 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to compute an LU factorization.");
2063 index.resize(keff_new);
2064 for (
int i=0; i < keff_new; i++) { index[i] = i; }
2065 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
2066 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
2070 if (keff != keff_new) {
2072 gcrodr_iter->setSize( keff, numBlocks_ );
2074 Teuchos::SerialDenseMatrix<int,ScalarType> b1( Teuchos::View, *H2_, recycledBlocks_+2, 1, 0, recycledBlocks_ );
2082 template<
class ScalarType,
class MV,
class OP>
2083 int GCRODRSolMgr<ScalarType,MV,OP,true>::getHarmonicVecs1(
int m,
2084 const Teuchos::SerialDenseMatrix<int,ScalarType>& HH,
2085 Teuchos::SerialDenseMatrix<int,ScalarType>& PP) {
2087 bool xtraVec =
false;
2088 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
2091 std::vector<MagnitudeType> wr(m), wi(m);
2094 Teuchos::SerialDenseMatrix<int,ScalarType> vr(m,m,
false);
2097 std::vector<MagnitudeType> w(m);
2100 std::vector<int> iperm(m);
2106 builtRecycleSpace_ =
true;
2109 Teuchos::SerialDenseMatrix<int, ScalarType> HHt( HH,
Teuchos::TRANS );
2110 Teuchos::SerialDenseVector<int, ScalarType> e_m( m );
2112 lapack.GESV(m, 1, HHt.values(), HHt.stride(), &iperm[0], e_m.values(), e_m.stride(), &info);
2113 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK GESV failed to compute a solution.");
2116 ScalarType d = HH(m, m-1) * HH(m, m-1);
2117 Teuchos::SerialDenseMatrix<int, ScalarType> harmHH( Teuchos::Copy, HH, m, m );
2118 for( i=0; i<m; ++i )
2119 harmHH(i, m-1) += d * e_m[i];
2128 std::vector<ScalarType> work(1);
2129 std::vector<MagnitudeType> rwork(2*m);
2132 lapack.GEEV(
'N',
'V', m, harmHH.values(), harmHH.stride(), &wr[0], &wi[0],
2133 vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
2135 lwork = std::abs (static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work[0])));
2136 work.resize( lwork );
2138 lapack.GEEV(
'N',
'V', m, harmHH.values(), harmHH.stride(), &wr[0], &wi[0],
2139 vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
2140 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK GEEV failed to compute eigensolutions.");
2143 for( i=0; i<m; ++i )
2144 w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] );
2147 this->sort(w, m, iperm);
2149 const bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
2152 for( i=0; i<recycledBlocks_; ++i ) {
2153 for( j=0; j<m; j++ ) {
2154 PP(j,i) = vr(j,iperm[i]);
2158 if(!scalarTypeIsComplex) {
2161 if (wi[iperm[recycledBlocks_-1]] != 0.0) {
2163 for ( i=0; i<recycledBlocks_; ++i ) {
2164 if (wi[iperm[i]] != 0.0)
2173 if (wi[iperm[recycledBlocks_-1]] > 0.0) {
2174 for( j=0; j<m; ++j ) {
2175 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
2179 for( j=0; j<m; ++j ) {
2180 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
2189 return recycledBlocks_+1;
2192 return recycledBlocks_;
2198 template<
class ScalarType,
class MV,
class OP>
2199 int GCRODRSolMgr<ScalarType,MV,OP,true>::getHarmonicVecs2(
int keffloc,
int m,
2200 const Teuchos::SerialDenseMatrix<int,ScalarType>& HH,
2201 const Teuchos::RCP<const MV>& VV,
2202 Teuchos::SerialDenseMatrix<int,ScalarType>& PP) {
2204 int m2 = HH.numCols();
2205 bool xtraVec =
false;
2206 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
2207 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
2208 std::vector<int> index;
2211 std::vector<MagnitudeType> wr(m2), wi(m2);
2214 std::vector<MagnitudeType> w(m2);
2217 Teuchos::SerialDenseMatrix<int,ScalarType> vr(m2,m2,
false);
2220 std::vector<int> iperm(m2);
2223 builtRecycleSpace_ =
true;
2228 Teuchos::SerialDenseMatrix<int,ScalarType> B(m2,m2,
false);
2233 Teuchos::SerialDenseMatrix<int,ScalarType> A_tmp( keffloc+m+1, keffloc+m );
2236 index.resize(keffloc);
2237 for (i=0; i<keffloc; ++i) { index[i] = i; }
2238 Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
2239 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
2240 Teuchos::SerialDenseMatrix<int,ScalarType> A11( Teuchos::View, A_tmp, keffloc, keffloc );
2241 MVT::MvTransMv( one, *Ctmp, *Utmp, A11 );
2244 Teuchos::SerialDenseMatrix<int,ScalarType> A21( Teuchos::View, A_tmp, m+1, keffloc, keffloc );
2246 for (i=0; i < m+1; i++) { index[i] = i; }
2247 Teuchos::RCP<const MV> Vp = MVT::CloneView( *VV, index );
2248 MVT::MvTransMv( one, *Vp, *Utmp, A21 );
2251 for( i=keffloc; i<keffloc+m; i++ ) {
2256 Teuchos::SerialDenseMatrix<int,ScalarType> A( m2, A_tmp.numCols() );
2257 A.multiply(
Teuchos::TRANS, Teuchos::NO_TRANS, one, HH, A_tmp, zero );
2265 char balanc=
'P', jobvl=
'N', jobvr=
'V', sense=
'N';
2266 int ld = A.numRows();
2268 int ldvl = ld, ldvr = ld;
2269 int info = 0,ilo = 0,ihi = 0;
2270 MagnitudeType abnrm = 0.0, bbnrm = 0.0;
2272 std::vector<ScalarType> beta(ld);
2273 std::vector<ScalarType> work(lwork);
2274 std::vector<MagnitudeType> rwork(lwork);
2275 std::vector<MagnitudeType> lscale(ld), rscale(ld);
2276 std::vector<MagnitudeType> rconde(ld), rcondv(ld);
2277 std::vector<int> iwork(ld+6);
2282 lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0],
2283 &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
2284 &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &rwork[0],
2285 &iwork[0], bwork, &info);
2286 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK GGEVX failed to compute eigensolutions.");
2290 for( i=0; i<ld; i++ ) {
2291 w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot (wr[i]*wr[i] + wi[i]*wi[i]) /
2292 Teuchos::ScalarTraits<ScalarType>::magnitude (beta[i]);
2296 this->sort(w,ld,iperm);
2298 const bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
2301 for( i=0; i<recycledBlocks_; i++ ) {
2302 for( j=0; j<ld; j++ ) {
2303 PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
2307 if(!scalarTypeIsComplex) {
2310 if (wi[iperm[ld-recycledBlocks_]] != 0.0) {
2312 for ( i=ld-recycledBlocks_; i<ld; i++ ) {
2313 if (wi[iperm[i]] != 0.0)
2322 if (wi[iperm[ld-recycledBlocks_]] > 0.0) {
2323 for( j=0; j<ld; j++ ) {
2324 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
2328 for( j=0; j<ld; j++ ) {
2329 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
2338 return recycledBlocks_+1;
2341 return recycledBlocks_;
2348 template<
class ScalarType,
class MV,
class OP>
2349 void GCRODRSolMgr<ScalarType,MV,OP,true>::sort(std::vector<MagnitudeType>& dlist,
int n, std::vector<int>& iperm) {
2350 int l, r, j, i, flag;
2352 MagnitudeType dRR, dK;
2379 if (dlist[j] > dlist[j - 1]) j = j + 1;
2381 if (dlist[j - 1] > dK) {
2382 dlist[i - 1] = dlist[j - 1];
2383 iperm[i - 1] = iperm[j - 1];
2397 dlist[r] = dlist[0];
2398 iperm[r] = iperm[0];
2413 template<
class ScalarType,
class MV,
class OP>
2415 std::ostringstream out;
2416 out <<
"Belos::GCRODRSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<
">";
2418 out <<
"Ortho Type: \"" << orthoType_ <<
"\"";
2419 out <<
", Num Blocks: " <<numBlocks_;
2420 out <<
", Num Recycle Blocks: " << recycledBlocks_;
2421 out <<
", Max Restarts: " << maxRestarts_;
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Collection of types and exceptions used within the Belos solvers.
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...
GCRODRSolMgrLAPACKFailure(const std::string &what_arg)
Class which manages the output and verbosity of the Belos solvers.
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().
GCRODRSolMgrLAPACKFailure is thrown when a nonzero value is retuned from an LAPACK call...
A factory class for generating StatusTestOutput objects.
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::ostream & printValidNames(std::ostream &out) const
Print all recognized MatOrthoManager names to the given ostream.
static const double convTol
Default convergence tolerance.
Belos::StatusTest class for specifying a maximum number of iterations.
int curDim
The current dimension of the reduction.
A factory class for generating StatusTestOutput objects.
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.
Pure virtual base class which describes the basic interface for a solver manager. ...
Teuchos::RCP< MV > V
The current Krylov basis.
Teuchos::RCP< MV > U
The recycled subspace and its projection.
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.
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.
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 setParameters(const Teuchos::RCP< Teuchos::ParameterList > ¶ms) override
Set the parameters the solver manager should use to solve the linear problem.
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...
virtual ~GCRODRSolMgr()
Destructor.
int setTolerance(MagnitudeType tolerance)
Set the value of the tolerance.
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.
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().
Class which defines basic traits for the operator type.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
Parent class to all Belos exceptions.
Base class for Belos::SolverManager subclasses which normally can only compile with ScalarType types ...
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::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
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.