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"
68#if defined(HAVE_TEUCHOSCORE_CXX11)
69# include <type_traits>
154 template<
class ScalarType,
class MV,
class OP,
155 const bool lapackSupportsScalarType =
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.");
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);
458 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >
problem_;
465 Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
sTest_;
468 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> >
expConvTest_, impConvTest_;
474 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> >
ortho_;
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";
517 Teuchos::RCP<MV> U_,
C_;
520 Teuchos::RCP<MV> U1_,
C1_;
523 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
H2_;
524 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
H_;
525 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
B_;
526 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
PP_;
527 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
HP_;
530 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
R_;
547template<
class ScalarType,
class MV,
class OP>
557template<
class ScalarType,
class MV,
class OP>
560 const Teuchos::RCP<Teuchos::ParameterList>& pl):
568 TEUCHOS_TEST_FOR_EXCEPTION(
569 problem == Teuchos::null, std::invalid_argument,
570 "Belos::GCRODRSolMgr constructor: The solver manager's "
571 "constructor needs the linear problem argument 'problem' "
580 if (! pl.is_null ()) {
586template<
class ScalarType,
class MV,
class OP>
588 outputStream_ = Teuchos::rcpFromRef(std::cout);
590 orthoKappa_ = orthoKappa_default_;
591 maxRestarts_ = maxRestarts_default_;
592 maxIters_ = maxIters_default_;
593 numBlocks_ = numBlocks_default_;
594 recycledBlocks_ = recycledBlocks_default_;
595 verbosity_ = verbosity_default_;
596 outputStyle_ = outputStyle_default_;
597 outputFreq_ = outputFreq_default_;
598 orthoType_ = orthoType_default_;
599 impResScale_ = impResScale_default_;
600 expResScale_ = expResScale_default_;
601 label_ = label_default_;
603 builtRecycleSpace_ =
false;
619template<
class ScalarType,
class MV,
class OP>
622setParameters (
const Teuchos::RCP<Teuchos::ParameterList> ¶ms)
624 using Teuchos::isParameterType;
625 using Teuchos::getParameter;
627 using Teuchos::ParameterList;
628 using Teuchos::parameterList;
631 using Teuchos::rcp_dynamic_cast;
632 using Teuchos::rcpFromRef;
633 using Teuchos::Exceptions::InvalidParameter;
634 using Teuchos::Exceptions::InvalidParameterName;
635 using Teuchos::Exceptions::InvalidParameterType;
639 RCP<const ParameterList> defaultParams = getValidParameters();
657 if (params_.is_null()) {
658 params_ = parameterList (*defaultParams);
666 if (params_ != params) {
672 params_ = parameterList (*params);
707 params_->validateParametersAndSetDefaults (*defaultParams);
711 if (params->isParameter (
"Maximum Restarts")) {
712 maxRestarts_ = params->get(
"Maximum Restarts", maxRestarts_default_);
715 params_->set (
"Maximum Restarts", maxRestarts_);
719 if (params->isParameter (
"Maximum Iterations")) {
720 maxIters_ = params->get (
"Maximum Iterations", maxIters_default_);
723 params_->set (
"Maximum Iterations", maxIters_);
724 if (! maxIterTest_.is_null())
725 maxIterTest_->setMaxIters (maxIters_);
729 if (params->isParameter (
"Num Blocks")) {
730 numBlocks_ = params->get (
"Num Blocks", numBlocks_default_);
731 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
732 "Belos::GCRODRSolMgr: The \"Num Blocks\" parameter must "
733 "be strictly positive, but you specified a value of "
734 << numBlocks_ <<
".");
736 params_->set (
"Num Blocks", numBlocks_);
740 if (params->isParameter (
"Num Recycled Blocks")) {
741 recycledBlocks_ = params->get (
"Num Recycled Blocks",
742 recycledBlocks_default_);
743 TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks_ <= 0, std::invalid_argument,
744 "Belos::GCRODRSolMgr: The \"Num Recycled Blocks\" "
745 "parameter must be strictly positive, but you specified "
746 "a value of " << recycledBlocks_ <<
".");
747 TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks_ >= numBlocks_, std::invalid_argument,
748 "Belos::GCRODRSolMgr: The \"Num Recycled Blocks\" "
749 "parameter must be less than the \"Num Blocks\" "
750 "parameter, but you specified \"Num Recycled Blocks\" "
751 "= " << recycledBlocks_ <<
" and \"Num Blocks\" = "
752 << numBlocks_ <<
".");
754 params_->set(
"Num Recycled Blocks", recycledBlocks_);
760 if (params->isParameter (
"Timer Label")) {
761 std::string tempLabel = params->get (
"Timer Label", label_default_);
764 if (tempLabel != label_) {
766 params_->set (
"Timer Label", label_);
767 std::string solveLabel = label_ +
": GCRODRSolMgr total solve time";
768#ifdef BELOS_TEUCHOS_TIME_MONITOR
769 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
771 if (ortho_ != Teuchos::null) {
772 ortho_->setLabel( label_ );
778 if (params->isParameter (
"Verbosity")) {
779 if (isParameterType<int> (*params,
"Verbosity")) {
780 verbosity_ = params->get (
"Verbosity", verbosity_default_);
782 verbosity_ = (int) getParameter<Belos::MsgType> (*params,
"Verbosity");
785 params_->set (
"Verbosity", verbosity_);
788 if (! printer_.is_null())
789 printer_->setVerbosity (verbosity_);
793 if (params->isParameter (
"Output Style")) {
794 if (isParameterType<int> (*params,
"Output Style")) {
795 outputStyle_ = params->get (
"Output Style", outputStyle_default_);
797 outputStyle_ = (int) getParameter<OutputType> (*params,
"Output Style");
801 params_->set (
"Output Style", outputStyle_);
819 if (params->isParameter (
"Output Stream")) {
821 outputStream_ = getParameter<RCP<std::ostream> > (*params,
"Output Stream");
822 }
catch (InvalidParameter&) {
823 outputStream_ = rcpFromRef (std::cout);
830 if (outputStream_.is_null()) {
831 outputStream_ = rcp (
new Teuchos::oblackholestream);
834 params_->set (
"Output Stream", outputStream_);
837 if (! printer_.is_null()) {
838 printer_->setOStream (outputStream_);
844 if (params->isParameter (
"Output Frequency")) {
845 outputFreq_ = params->get (
"Output Frequency", outputFreq_default_);
849 params_->set(
"Output Frequency", outputFreq_);
850 if (! outputTest_.is_null())
851 outputTest_->setOutputFrequency (outputFreq_);
858 if (printer_.is_null()) {
869 bool changedOrthoType =
false;
870 if (params->isParameter (
"Orthogonalization")) {
871 const std::string& tempOrthoType =
872 params->get (
"Orthogonalization", orthoType_default_);
875 std::ostringstream os;
876 os <<
"Belos::GCRODRSolMgr: Invalid orthogonalization name \""
877 << tempOrthoType <<
"\". The following are valid options "
878 <<
"for the \"Orthogonalization\" name parameter: ";
880 throw std::invalid_argument (os.str());
882 if (tempOrthoType != orthoType_) {
883 changedOrthoType =
true;
884 orthoType_ = tempOrthoType;
886 params_->set (
"Orthogonalization", orthoType_);
902 RCP<ParameterList> orthoParams;
905 using Teuchos::sublist;
907 const std::string paramName (
"Orthogonalization Parameters");
910 orthoParams = sublist (params_, paramName,
true);
911 }
catch (InvalidParameter&) {
918 orthoParams = sublist (params_, paramName,
true);
921 TEUCHOS_TEST_FOR_EXCEPTION(orthoParams.is_null(), std::logic_error,
922 "Failed to get orthogonalization parameters. "
923 "Please report this bug to the Belos developers.");
928 if (ortho_.is_null() || changedOrthoType) {
934 label_, orthoParams);
942 typedef Teuchos::ParameterListAcceptor PLA;
943 RCP<PLA> pla = rcp_dynamic_cast<PLA> (ortho_);
949 label_, orthoParams);
951 pla->setParameterList (orthoParams);
963 if (params->isParameter (
"Orthogonalization Constant")) {
965 if (params->isType<
MagnitudeType> (
"Orthogonalization Constant")) {
966 orthoKappa = params->get (
"Orthogonalization Constant", orthoKappa);
969 orthoKappa = params->get (
"Orthogonalization Constant", orthoKappa_default_);
972 if (orthoKappa > 0) {
973 orthoKappa_ = orthoKappa;
975 params_->set(
"Orthogonalization Constant", orthoKappa_);
977 if (orthoType_ ==
"DGKS" && ! ortho_.is_null()) {
984 rcp_dynamic_cast<ortho_man_type>(ortho_)->
setDepTol (orthoKappa_);
994 if (params->isParameter(
"Convergence Tolerance")) {
995 if (params->isType<
MagnitudeType> (
"Convergence Tolerance")) {
996 convTol_ = params->get (
"Convergence Tolerance",
1004 params_->set (
"Convergence Tolerance", convTol_);
1005 if (! impConvTest_.is_null())
1007 if (! expConvTest_.is_null())
1008 expConvTest_->setTolerance (convTol_);
1012 if (params->isParameter (
"Implicit Residual Scaling")) {
1013 std::string tempImpResScale =
1014 getParameter<std::string> (*params,
"Implicit Residual Scaling");
1017 if (impResScale_ != tempImpResScale) {
1019 impResScale_ = tempImpResScale;
1022 params_->set(
"Implicit Residual Scaling", impResScale_);
1032 if (! impConvTest_.is_null()) {
1038 impConvTest_ = null;
1045 if (params->isParameter(
"Explicit Residual Scaling")) {
1046 std::string tempExpResScale =
1047 getParameter<std::string> (*params,
"Explicit Residual Scaling");
1050 if (expResScale_ != tempExpResScale) {
1052 expResScale_ = tempExpResScale;
1055 params_->set(
"Explicit Residual Scaling", expResScale_);
1058 if (! expConvTest_.is_null()) {
1064 expConvTest_ = null;
1075 if (maxIterTest_.is_null())
1080 if (impConvTest_.is_null()) {
1081 impConvTest_ = rcp (
new StatusTestResNorm_t (convTol_));
1087 if (expConvTest_.is_null()) {
1088 expConvTest_ = rcp (
new StatusTestResNorm_t (convTol_));
1089 expConvTest_->defineResForm (StatusTestResNorm_t::Explicit,
Belos::TwoNorm);
1095 if (convTest_.is_null()) {
1096 convTest_ = rcp (
new StatusTestCombo_t (StatusTestCombo_t::SEQ,
1104 sTest_ = rcp (
new StatusTestCombo_t (StatusTestCombo_t::OR,
1110 outputTest_ = stoFactory.
create (printer_, sTest_, outputFreq_,
1114 std::string solverDesc =
" GCRODR ";
1115 outputTest_->setSolverDesc( solverDesc );
1118 if (timerSolve_.is_null()) {
1119 std::string solveLabel = label_ +
": GCRODRSolMgr total solve time";
1120#ifdef BELOS_TEUCHOS_TIME_MONITOR
1121 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
1130template<
class ScalarType,
class MV,
class OP>
1131Teuchos::RCP<const Teuchos::ParameterList>
1134 using Teuchos::ParameterList;
1135 using Teuchos::parameterList;
1138 static RCP<const ParameterList> validPL;
1139 if (is_null(validPL)) {
1140 RCP<ParameterList> pl = parameterList ();
1144 "The relative residual tolerance that needs to be achieved by the\n"
1145 "iterative solver in order for the linear system to be declared converged.");
1146 pl->set(
"Maximum Restarts",
static_cast<int>(maxRestarts_default_),
1147 "The maximum number of cycles allowed for each\n"
1148 "set of RHS solved.");
1149 pl->set(
"Maximum Iterations",
static_cast<int>(maxIters_default_),
1150 "The maximum number of iterations allowed for each\n"
1151 "set of RHS solved.");
1155 pl->set(
"Block Size",
static_cast<int>(blockSize_default_),
1156 "Block Size Parameter -- currently must be 1 for GCRODR");
1157 pl->set(
"Num Blocks",
static_cast<int>(numBlocks_default_),
1158 "The maximum number of vectors allowed in the Krylov subspace\n"
1159 "for each set of RHS solved.");
1160 pl->set(
"Num Recycled Blocks",
static_cast<int>(recycledBlocks_default_),
1161 "The maximum number of vectors in the recycled subspace." );
1162 pl->set(
"Verbosity",
static_cast<int>(verbosity_default_),
1163 "What type(s) of solver information should be outputted\n"
1164 "to the output stream.");
1165 pl->set(
"Output Style",
static_cast<int>(outputStyle_default_),
1166 "What style is used for the solver information outputted\n"
1167 "to the output stream.");
1168 pl->set(
"Output Frequency",
static_cast<int>(outputFreq_default_),
1169 "How often convergence information should be outputted\n"
1170 "to the output stream.");
1171 pl->set(
"Output Stream", Teuchos::rcpFromRef(std::cout),
1172 "A reference-counted pointer to the output stream where all\n"
1173 "solver output is sent.");
1174 pl->set(
"Implicit Residual Scaling",
static_cast<const char *
>(impResScale_default_),
1175 "The type of scaling used in the implicit residual convergence test.");
1176 pl->set(
"Explicit Residual Scaling",
static_cast<const char *
>(expResScale_default_),
1177 "The type of scaling used in the explicit residual convergence test.");
1178 pl->set(
"Timer Label",
static_cast<const char *
>(label_default_),
1179 "The string to use as a prefix for the timer labels.");
1182 pl->set(
"Orthogonalization",
static_cast<const char *
>(orthoType_default_),
1183 "The type of orthogonalization to use. Valid options: " +
1185 RCP<const ParameterList> orthoParams =
1187 pl->set (
"Orthogonalization Parameters", *orthoParams,
1188 "Parameters specific to the type of orthogonalization used.");
1190 pl->set(
"Orthogonalization Constant",
static_cast<MagnitudeType>(orthoKappa_default_),
1191 "When using DGKS orthogonalization: the \"depTol\" constant, used "
1192 "to determine whether another step of classical Gram-Schmidt is "
1193 "necessary. Otherwise ignored.");
1200template<
class ScalarType,
class MV,
class OP>
1203 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1206 Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
1207 if (rhsMV == Teuchos::null) {
1214 TEUCHOS_TEST_FOR_EXCEPTION(
static_cast<ptrdiff_t
>(numBlocks_) > MVT::GetGlobalLength(*rhsMV),std::invalid_argument,
1215 "Belos::GCRODRSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
1218 if (U_ == Teuchos::null) {
1219 U_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1223 if (MVT::GetNumberVecs(*U_) < recycledBlocks_+1) {
1224 Teuchos::RCP<const MV> tmp = U_;
1225 U_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1230 if (C_ == Teuchos::null) {
1231 C_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1235 if (MVT::GetNumberVecs(*C_) < recycledBlocks_+1) {
1236 Teuchos::RCP<const MV> tmp = C_;
1237 C_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1242 if (V_ == Teuchos::null) {
1243 V_ = MVT::Clone( *rhsMV, numBlocks_+1 );
1247 if (MVT::GetNumberVecs(*V_) < numBlocks_+1) {
1248 Teuchos::RCP<const MV> tmp = V_;
1249 V_ = MVT::Clone( *tmp, numBlocks_+1 );
1254 if (U1_ == Teuchos::null) {
1255 U1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1259 if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1260 Teuchos::RCP<const MV> tmp = U1_;
1261 U1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1266 if (C1_ == Teuchos::null) {
1267 C1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1271 if (MVT::GetNumberVecs(*C1_) < recycledBlocks_+1) {
1272 Teuchos::RCP<const MV> tmp = C1_;
1273 C1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1278 if (r_ == Teuchos::null)
1279 r_ = MVT::Clone( *rhsMV, 1 );
1282 tau_.resize(recycledBlocks_+1);
1285 work_.resize(recycledBlocks_+1);
1288 ipiv_.resize(recycledBlocks_+1);
1291 if (H2_ == Teuchos::null)
1292 H2_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 ) );
1294 if ( (H2_->numRows() != numBlocks_+recycledBlocks_+2) || (H2_->numCols() != numBlocks_+recycledBlocks_+1) )
1295 H2_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
1297 H2_->putScalar(zero);
1300 if (R_ == Teuchos::null)
1301 R_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycledBlocks_+1, recycledBlocks_+1 ) );
1303 if ( (R_->numRows() != recycledBlocks_+1) || (R_->numCols() != recycledBlocks_+1) )
1304 R_->reshape( recycledBlocks_+1, recycledBlocks_+1 );
1306 R_->putScalar(zero);
1309 if (PP_ == Teuchos::null)
1310 PP_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, recycledBlocks_+1 ) );
1312 if ( (PP_->numRows() != numBlocks_+recycledBlocks_+2) || (PP_->numCols() != recycledBlocks_+1) )
1313 PP_->reshape( numBlocks_+recycledBlocks_+2, recycledBlocks_+1 );
1317 if (HP_ == Teuchos::null)
1318 HP_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 ) );
1320 if ( (HP_->numRows() != numBlocks_+recycledBlocks_+2) || (HP_->numCols() != numBlocks_+recycledBlocks_+1) )
1321 HP_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
1329template<
class ScalarType,
class MV,
class OP>
1337 if (!isSet_) { setParameters( params_ ); }
1339 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1340 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1341 std::vector<int> index(numBlocks_+1);
1343 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,
GCRODRSolMgrLinearProblemFailure,
"Belos::GCRODRSolMgr::solve(): Linear problem is not a valid object.");
1345 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),
GCRODRSolMgrLinearProblemFailure,
"Belos::GCRODRSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
1348 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1349 std::vector<int> currIdx(1);
1353 problem_->setLSIndex( currIdx );
1356 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1357 if (
static_cast<ptrdiff_t
>(numBlocks_) > dim) {
1358 numBlocks_ = Teuchos::as<int>(dim);
1360 "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
1361 " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
1362 params_->set(
"Num Blocks", numBlocks_);
1366 bool isConverged =
true;
1369 initializeStateStorage();
1373 Teuchos::ParameterList plist;
1375 plist.set(
"Num Blocks",numBlocks_);
1376 plist.set(
"Recycled Blocks",recycledBlocks_);
1381 RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_iter;
1384 int prime_iterations = 0;
1388#ifdef BELOS_TEUCHOS_TIME_MONITOR
1389 Teuchos::TimeMonitor slvtimer(*timerSolve_);
1392 while ( numRHS2Solve > 0 ) {
1395 builtRecycleSpace_ =
false;
1398 outputTest_->reset();
1406 "Belos::GCRODRSolMgr::solve(): Requested size of recycled subspace is not consistent with the current recycle subspace.");
1408 printer_->stream(
Debug) <<
" Now solving RHS index " << currIdx[0] <<
" using recycled subspace of dimension " << keff << std::endl << std::endl;
1411 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1412 RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1413 RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1414 problem_->apply( *Utmp, *Ctmp );
1416 RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1420 Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff, keff );
1421 int rank = ortho_->normalize(*Ctmp, rcp(&Rtmp,
false));
1423 TEUCHOS_TEST_FOR_EXCEPTION(rank != keff,
GCRODRSolMgrOrthoFailure,
"Belos::GCRODRSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
1428 ipiv_.resize(Rtmp.numRows());
1429 lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
1430 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1433 int lwork = Rtmp.numRows();
1434 work_.resize(lwork);
1435 lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1436 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
1439 MVT::MvTimesMatAddMv( one, *Utmp, Rtmp, zero, *U1tmp );
1444 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1445 Ctmp = MVT::CloneViewNonConst( *C_, index );
1446 Utmp = MVT::CloneView( *U_, index );
1449 Teuchos::SerialDenseMatrix<int,ScalarType> Ctr(keff,1);
1450 problem_->computeCurrPrecResVec( &*r_ );
1451 MVT::MvTransMv( one, *Ctmp, *r_, Ctr );
1454 RCP<MV> update = MVT::Clone( *problem_->getCurrLHSVec(), 1 );
1455 MVT::MvInit( *update, 0.0 );
1456 MVT::MvTimesMatAddMv( one, *Utmp, Ctr, one, *update );
1457 problem_->updateSolution( update,
true );
1460 MVT::MvTimesMatAddMv( -one, *Ctmp, Ctr, one, *r_ );
1463 prime_iterations = 0;
1469 printer_->stream(
Debug) <<
" No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
1471 Teuchos::ParameterList primeList;
1474 primeList.set(
"Num Blocks",numBlocks_);
1475 primeList.set(
"Recycled Blocks",0);
1478 RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_prime_iter;
1482 problem_->computeCurrPrecResVec( &*r_ );
1483 index.resize( 1 ); index[0] = 0;
1484 RCP<MV> v0 = MVT::CloneViewNonConst( *V_, index );
1485 MVT::SetBlock(*r_,index,*v0);
1489 index.resize( numBlocks_+1 );
1490 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1491 newstate.
V = MVT::CloneViewNonConst( *V_, index );
1492 newstate.
U = Teuchos::null;
1493 newstate.
C = Teuchos::null;
1494 newstate.
H = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, recycledBlocks_+1, recycledBlocks_+1 ) );
1495 newstate.
B = Teuchos::null;
1497 gcrodr_prime_iter->initialize(newstate);
1500 bool primeConverged =
false;
1502 gcrodr_prime_iter->iterate();
1505 if ( convTest_->getStatus() ==
Passed ) {
1507 primeConverged =
true;
1512 gcrodr_prime_iter->updateLSQR( gcrodr_prime_iter->getCurSubspaceDim() );
1515 sTest_->checkStatus( &*gcrodr_prime_iter );
1516 if (convTest_->getStatus() ==
Passed)
1517 primeConverged =
true;
1519 catch (
const std::exception &e) {
1520 printer_->stream(
Errors) <<
"Error! Caught exception in GCRODRIter::iterate() at iteration "
1521 << gcrodr_prime_iter->getNumIters() << std::endl
1522 << e.what() << std::endl;
1526 prime_iterations = gcrodr_prime_iter->getNumIters();
1529 RCP<MV> update = gcrodr_prime_iter->getCurrentUpdate();
1530 problem_->updateSolution( update,
true );
1533 newstate = gcrodr_prime_iter->getState();
1541 if (recycledBlocks_ < p+1) {
1543 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PPtmp = rcp (
new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, p, recycledBlocks_+1 ) );
1545 keff = getHarmonicVecs1( p, *newstate.
H, *PPtmp );
1547 PPtmp = rcp (
new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, p, keff ) );
1550 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1551 RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1552 RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1553 RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1555 for (
int ii=0; ii < p; ++ii) { index[ii] = ii; }
1556 RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1560 MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
1567 Teuchos::SerialDenseMatrix<int,ScalarType> Htmp( Teuchos::View, *H2_, p+1, p, recycledBlocks_+1,recycledBlocks_+1);
1568 Teuchos::SerialDenseMatrix<int,ScalarType> HPtmp( Teuchos::View, *HP_, p+1, keff );
1569 HPtmp.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, Htmp, *PPtmp, zero );
1576 lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
1577 HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
1578 TEUCHOS_TEST_FOR_EXCEPTION(
1580 " LAPACK's _GEQRF failed to compute a workspace size.");
1588 lwork = std::abs (
static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work_[0])));
1589 work_.resize (lwork);
1590 lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
1591 HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
1592 TEUCHOS_TEST_FOR_EXCEPTION(
1594 " LAPACK's _GEQRF failed to compute a QR factorization.");
1598 Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff, keff );
1599 for (
int ii = 0; ii < keff; ++ii) {
1600 for (
int jj = ii; jj < keff; ++jj) {
1601 Rtmp(ii,jj) = HPtmp(ii,jj);
1608 lapack.UNGQR (HPtmp.numRows (), HPtmp.numCols (), HPtmp.numCols (),
1609 HPtmp.values (), HPtmp.stride (), &tau_[0], &work_[0],
1611 TEUCHOS_TEST_FOR_EXCEPTION(
1613 "LAPACK's _UNGQR failed to construct the Q factor.");
1618 index.resize (p + 1);
1619 for (
int ii = 0; ii < (p+1); ++ii) {
1622 Vtmp = MVT::CloneView( *V_, index );
1623 MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
1630 ipiv_.resize(Rtmp.numRows());
1631 lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
1632 TEUCHOS_TEST_FOR_EXCEPTION(
1634 "LAPACK's _GETRF failed to compute an LU factorization.");
1643 lwork = Rtmp.numRows();
1644 work_.resize(lwork);
1645 lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1646 TEUCHOS_TEST_FOR_EXCEPTION(
1648 "LAPACK's _GETRI failed to invert triangular matrix.");
1651 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
1653 printer_->stream(
Debug)
1654 <<
" Generated recycled subspace using RHS index " << currIdx[0]
1655 <<
" of dimension " << keff << std::endl << std::endl;
1660 if (primeConverged) {
1662 problem_->setCurrLS();
1666 if (numRHS2Solve > 0) {
1668 problem_->setLSIndex (currIdx);
1671 currIdx.resize (numRHS2Solve);
1681 gcrodr_iter->setSize( keff, numBlocks_ );
1684 gcrodr_iter->resetNumIters(prime_iterations);
1687 outputTest_->resetNumCalls();
1690 problem_->computeCurrPrecResVec( &*r_ );
1691 index.resize( 1 ); index[0] = 0;
1692 RCP<MV> v0 = MVT::CloneViewNonConst( *V_, index );
1693 MVT::SetBlock(*r_,index,*v0);
1697 index.resize( numBlocks_+1 );
1698 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1699 newstate.
V = MVT::CloneViewNonConst( *V_, index );
1700 index.resize( keff );
1701 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1702 newstate.
C = MVT::CloneViewNonConst( *C_, index );
1703 newstate.
U = MVT::CloneViewNonConst( *U_, index );
1704 newstate.
B = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, keff, numBlocks_, 0, keff ) );
1705 newstate.
H = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, keff, keff ) );
1707 gcrodr_iter->initialize(newstate);
1710 int numRestarts = 0;
1715 gcrodr_iter->iterate();
1722 if ( convTest_->getStatus() ==
Passed ) {
1731 else if ( maxIterTest_->getStatus() ==
Passed ) {
1733 isConverged =
false;
1741 else if ( gcrodr_iter->getCurSubspaceDim() == gcrodr_iter->getMaxSubspaceDim() ) {
1746 RCP<MV> update = gcrodr_iter->getCurrentUpdate();
1747 problem_->updateSolution( update,
true );
1749 buildRecycleSpace2(gcrodr_iter);
1751 printer_->stream(
Debug)
1752 <<
" Generated new recycled subspace using RHS index "
1753 << currIdx[0] <<
" of dimension " << keff << std::endl
1757 if (numRestarts >= maxRestarts_) {
1758 isConverged =
false;
1763 printer_->stream(
Debug)
1764 <<
" Performing restart number " << numRestarts <<
" of "
1765 << maxRestarts_ << std::endl << std::endl;
1768 problem_->computeCurrPrecResVec( &*r_ );
1769 index.resize( 1 ); index[0] = 0;
1770 RCP<MV> v00 = MVT::CloneViewNonConst( *V_, index );
1771 MVT::SetBlock(*r_,index,*v00);
1775 index.resize( numBlocks_+1 );
1776 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1777 restartState.
V = MVT::CloneViewNonConst( *V_, index );
1778 index.resize( keff );
1779 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1780 restartState.
U = MVT::CloneViewNonConst( *U_, index );
1781 restartState.
C = MVT::CloneViewNonConst( *C_, index );
1782 restartState.
B = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, keff, numBlocks_, 0, keff ) );
1783 restartState.
H = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, keff, keff ) );
1785 gcrodr_iter->initialize(restartState);
1798 TEUCHOS_TEST_FOR_EXCEPTION(
1799 true, std::logic_error,
"Belos::GCRODRSolMgr::solve: "
1800 "Invalid return from GCRODRIter::iterate().");
1805 gcrodr_iter->updateLSQR( gcrodr_iter->getCurSubspaceDim() );
1808 sTest_->checkStatus( &*gcrodr_iter );
1809 if (convTest_->getStatus() !=
Passed)
1810 isConverged =
false;
1813 catch (
const std::exception& e) {
1815 <<
"Error! Caught exception in GCRODRIter::iterate() at iteration "
1816 << gcrodr_iter->getNumIters() << std::endl << e.what() << std::endl;
1823 RCP<MV> update = gcrodr_iter->getCurrentUpdate();
1824 problem_->updateSolution( update,
true );
1827 problem_->setCurrLS();
1832 if (!builtRecycleSpace_) {
1833 buildRecycleSpace2(gcrodr_iter);
1834 printer_->stream(
Debug)
1835 <<
" Generated new recycled subspace using RHS index " << currIdx[0]
1836 <<
" of dimension " << keff << std::endl << std::endl;
1841 if (numRHS2Solve > 0) {
1843 problem_->setLSIndex (currIdx);
1846 currIdx.resize (numRHS2Solve);
1854#ifdef BELOS_TEUCHOS_TIME_MONITOR
1859 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails) );
1863 numIters_ = maxIterTest_->getNumIters ();
1875 const std::vector<MagnitudeType>* pTestValues = expConvTest_->getTestValue();
1876 if (pTestValues == NULL || pTestValues->size() < 1) {
1877 pTestValues = impConvTest_->getTestValue();
1879 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1880 "Belos::GCRODRSolMgr::solve(): The implicit convergence test's getTestValue() "
1881 "method returned NULL. Please report this bug to the Belos developers.");
1882 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1883 "Belos::GCRODRSolMgr::solve(): The implicit convergence test's getTestValue() "
1884 "method returned a vector of length zero. Please report this bug to the "
1885 "Belos developers.");
1890 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1897template<
class ScalarType,
class MV,
class OP>
1900 MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
1901 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1903 std::vector<MagnitudeType> d(keff);
1904 std::vector<ScalarType> dscalar(keff);
1905 std::vector<int> index(numBlocks_+1);
1917 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1918 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1920 dscalar.resize(keff);
1921 MVT::MvNorm( *Utmp, d );
1922 for (
int i=0; i<keff; ++i) {
1924 dscalar[i] = (ScalarType)d[i];
1926 MVT::MvScale( *Utmp, dscalar );
1930 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H2tmp = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, p+keff+1, p+keff ) );
1933 for (
int i=0; i<keff; ++i) {
1934 (*H2tmp)(i,i) = d[i];
1941 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p+keff, recycledBlocks_+1 );
1942 keff_new = getHarmonicVecs2( keff, p, *H2tmp, oldState.
V, PPtmp );
1949 Teuchos::RCP<MV> U1tmp;
1951 index.resize( keff );
1952 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1953 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1954 index.resize( keff_new );
1955 for (
int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
1956 U1tmp = MVT::CloneViewNonConst( *U1_, index );
1957 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, keff, keff_new );
1958 MVT::MvTimesMatAddMv( one, *Utmp, PPtmp, zero, *U1tmp );
1964 for (
int ii=0; ii < p; ii++) { index[ii] = ii; }
1965 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1966 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p, keff_new, keff );
1967 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
1971 Teuchos::SerialDenseMatrix<int,ScalarType> HPtmp( Teuchos::View, *HP_, p+keff+1, keff_new );
1973 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p+keff, keff_new );
1974 HPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*H2tmp,PPtmp,zero);
1978 int info = 0, lwork = -1;
1979 tau_.resize (keff_new);
1980 lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
1981 HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
1982 TEUCHOS_TEST_FOR_EXCEPTION(
1984 "LAPACK's _GEQRF failed to compute a workspace size.");
1990 lwork = std::abs (
static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work_[0])));
1991 work_.resize (lwork);
1992 lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
1993 HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
1994 TEUCHOS_TEST_FOR_EXCEPTION(
1996 "LAPACK's _GEQRF failed to compute a QR factorization.");
2000 Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff_new, keff_new );
2001 for(
int i=0;i<keff_new;i++) {
for(
int j=i;j<keff_new;j++) Rtmp(i,j) = HPtmp(i,j); }
2007 lapack.UNGQR (HPtmp.numRows (), HPtmp.numCols (), HPtmp.numCols (),
2008 HPtmp.values (), HPtmp.stride (), &tau_[0], &work_[0],
2010 TEUCHOS_TEST_FOR_EXCEPTION(
2012 "LAPACK's _UNGQR failed to construct the Q factor.");
2019 Teuchos::RCP<MV> C1tmp;
2022 for (
int i=0; i < keff; i++) { index[i] = i; }
2023 Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
2024 index.resize(keff_new);
2025 for (
int i=0; i < keff_new; i++) { index[i] = i; }
2026 C1tmp = MVT::CloneViewNonConst( *C1_, index );
2027 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *HP_, keff, keff_new );
2028 MVT::MvTimesMatAddMv( one, *Ctmp, PPtmp, zero, *C1tmp );
2032 index.resize( p+1 );
2033 for (
int i=0; i < p+1; ++i) { index[i] = i; }
2034 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
2035 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *HP_, p+1, keff_new, keff, 0 );
2036 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *C1tmp );
2045 ipiv_.resize(Rtmp.numRows());
2046 lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
2047 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
2050 lwork = Rtmp.numRows();
2051 work_.resize(lwork);
2052 lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
2053 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to compute an LU factorization.");
2056 index.resize(keff_new);
2057 for (
int i=0; i < keff_new; i++) { index[i] = i; }
2058 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
2059 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
2063 if (keff != keff_new) {
2065 gcrodr_iter->setSize( keff, numBlocks_ );
2067 Teuchos::SerialDenseMatrix<int,ScalarType> b1( Teuchos::View, *H2_, recycledBlocks_+2, 1, 0, recycledBlocks_ );
2075template<
class ScalarType,
class MV,
class OP>
2077 const Teuchos::SerialDenseMatrix<int,ScalarType>& HH,
2078 Teuchos::SerialDenseMatrix<int,ScalarType>& PP) {
2080 bool xtraVec =
false;
2081 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
2084 std::vector<MagnitudeType> wr(m), wi(m);
2087 Teuchos::SerialDenseMatrix<int,ScalarType> vr(m,m,
false);
2090 std::vector<MagnitudeType> w(m);
2093 std::vector<int> iperm(m);
2099 builtRecycleSpace_ =
true;
2102 Teuchos::SerialDenseMatrix<int, ScalarType> HHt( HH, Teuchos::TRANS );
2103 Teuchos::SerialDenseVector<int, ScalarType> e_m( m );
2105 lapack.GESV(m, 1, HHt.values(), HHt.stride(), &iperm[0], e_m.values(), e_m.stride(), &info);
2106 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK GESV failed to compute a solution.");
2109 ScalarType d = HH(m, m-1) * HH(m, m-1);
2110 Teuchos::SerialDenseMatrix<int, ScalarType> harmHH( Teuchos::Copy, HH, m, m );
2111 for( i=0; i<m; ++i )
2112 harmHH(i, m-1) += d * e_m[i];
2121 std::vector<ScalarType> work(1);
2122 std::vector<MagnitudeType> rwork(2*m);
2125 lapack.GEEV(
'N',
'V', m, harmHH.values(), harmHH.stride(), &wr[0], &wi[0],
2126 vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
2128 lwork = std::abs (
static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work[0])));
2129 work.resize( lwork );
2131 lapack.GEEV(
'N',
'V', m, harmHH.values(), harmHH.stride(), &wr[0], &wi[0],
2132 vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
2133 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK GEEV failed to compute eigensolutions.");
2136 for( i=0; i<m; ++i )
2137 w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] );
2140 this->sort(w, m, iperm);
2142 const bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
2145 for( i=0; i<recycledBlocks_; ++i ) {
2146 for( j=0; j<m; j++ ) {
2147 PP(j,i) = vr(j,iperm[i]);
2151 if(!scalarTypeIsComplex) {
2154 if (wi[iperm[recycledBlocks_-1]] != 0.0) {
2156 for ( i=0; i<recycledBlocks_; ++i ) {
2157 if (wi[iperm[i]] != 0.0)
2166 if (wi[iperm[recycledBlocks_-1]] > 0.0) {
2167 for( j=0; j<m; ++j ) {
2168 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
2172 for( j=0; j<m; ++j ) {
2173 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
2182 return recycledBlocks_+1;
2185 return recycledBlocks_;
2191template<
class ScalarType,
class MV,
class OP>
2193 const Teuchos::SerialDenseMatrix<int,ScalarType>& HH,
2194 const Teuchos::RCP<const MV>& VV,
2195 Teuchos::SerialDenseMatrix<int,ScalarType>& PP) {
2197 int m2 = HH.numCols();
2198 bool xtraVec =
false;
2199 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
2200 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
2201 std::vector<int> index;
2204 std::vector<MagnitudeType> wr(m2), wi(m2);
2207 std::vector<MagnitudeType> w(m2);
2210 Teuchos::SerialDenseMatrix<int,ScalarType> vr(m2,m2,
false);
2213 std::vector<int> iperm(m2);
2216 builtRecycleSpace_ =
true;
2221 Teuchos::SerialDenseMatrix<int,ScalarType> B(m2,m2,
false);
2222 B.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,HH,HH,zero);
2226 Teuchos::SerialDenseMatrix<int,ScalarType> A_tmp( keffloc+m+1, keffloc+m );
2229 index.resize(keffloc);
2230 for (i=0; i<keffloc; ++i) { index[i] = i; }
2231 Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
2232 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
2233 Teuchos::SerialDenseMatrix<int,ScalarType> A11( Teuchos::View, A_tmp, keffloc, keffloc );
2234 MVT::MvTransMv( one, *Ctmp, *Utmp, A11 );
2237 Teuchos::SerialDenseMatrix<int,ScalarType> A21( Teuchos::View, A_tmp, m+1, keffloc, keffloc );
2239 for (i=0; i < m+1; i++) { index[i] = i; }
2240 Teuchos::RCP<const MV> Vp = MVT::CloneView( *VV, index );
2241 MVT::MvTransMv( one, *Vp, *Utmp, A21 );
2244 for( i=keffloc; i<keffloc+m; i++ ) {
2249 Teuchos::SerialDenseMatrix<int,ScalarType> A( m2, A_tmp.numCols() );
2250 A.multiply( Teuchos::TRANS, Teuchos::NO_TRANS, one, HH, A_tmp, zero );
2258 char balanc=
'P', jobvl=
'N', jobvr=
'V', sense=
'N';
2259 int ld = A.numRows();
2261 int ldvl = ld, ldvr = ld;
2262 int info = 0,ilo = 0,ihi = 0;
2265 std::vector<ScalarType> beta(ld);
2266 std::vector<ScalarType> work(lwork);
2267 std::vector<MagnitudeType> rwork(lwork);
2268 std::vector<MagnitudeType> lscale(ld), rscale(ld);
2269 std::vector<MagnitudeType> rconde(ld), rcondv(ld);
2270 std::vector<int> iwork(ld+6);
2275 lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0],
2276 &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
2277 &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &rwork[0],
2278 &iwork[0], bwork, &info);
2279 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK GGEVX failed to compute eigensolutions.");
2283 for( i=0; i<ld; i++ ) {
2284 w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot (wr[i]*wr[i] + wi[i]*wi[i]) /
2285 Teuchos::ScalarTraits<ScalarType>::magnitude (beta[i]);
2289 this->sort(w,ld,iperm);
2291 const bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
2294 for( i=0; i<recycledBlocks_; i++ ) {
2295 for( j=0; j<ld; j++ ) {
2296 PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
2300 if(!scalarTypeIsComplex) {
2303 if (wi[iperm[ld-recycledBlocks_]] != 0.0) {
2305 for ( i=ld-recycledBlocks_; i<ld; i++ ) {
2306 if (wi[iperm[i]] != 0.0)
2315 if (wi[iperm[ld-recycledBlocks_]] > 0.0) {
2316 for( j=0; j<ld; j++ ) {
2317 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
2321 for( j=0; j<ld; j++ ) {
2322 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
2331 return recycledBlocks_+1;
2334 return recycledBlocks_;
2341template<
class ScalarType,
class MV,
class OP>
2343 int l, r, j, i, flag;
2372 if (dlist[j] > dlist[j - 1]) j = j + 1;
2374 if (dlist[j - 1] > dK) {
2375 dlist[i - 1] = dlist[j - 1];
2376 iperm[i - 1] = iperm[j - 1];
2390 dlist[r] = dlist[0];
2391 iperm[r] = iperm[0];
2406template<
class ScalarType,
class MV,
class OP>
2408 std::ostringstream out;
2409 out <<
"Belos::GCRODRSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<
">";
2411 out <<
"Ortho Type: \"" << orthoType_ <<
"\"";
2412 out <<
", Num Blocks: " <<numBlocks_;
2413 out <<
", Num Recycle Blocks: " << recycledBlocks_;
2414 out <<
", Max Restarts: " << maxRestarts_;
Belos concrete class for performing the block, flexible GMRES iteration.
Belos header file which uses auto-configuration information to include necessary C++ headers.
Belos concrete class for performing the GCRO-DR iteration.
Class which describes the linear problem to be solved by the iterative solver.
Class which manages the output and verbosity of the Belos solvers.
Pure virtual base class which describes the basic interface for a solver manager.
Belos::StatusTest for logically combining several status tests.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Belos::StatusTest class for specifying a maximum number of iterations.
A factory class for generating StatusTestOutput objects.
Collection of types and exceptions used within the Belos solvers.
Parent class to all Belos exceptions.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
void setDepTol(const MagnitudeType dep_tol)
Set parameter for re-orthogonalization threshhold.
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
Base class for Belos::SolverManager subclasses which normally can only compile with ScalarType types ...
GCRODRIterOrthoFailure is thrown when the GCRODRIter object is unable to compute independent directio...
This class implements the GCRODR iteration, where a single-std::vector Krylov subspace is constructed...
GCRODRSolMgrLAPACKFailure is thrown when a nonzero value is retuned from an LAPACK call.
GCRODRSolMgrLAPACKFailure(const std::string &what_arg)
GCRODRSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i....
GCRODRSolMgrLinearProblemFailure(const std::string &what_arg)
GCRODRSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate orthonorm...
GCRODRSolMgrOrthoFailure(const std::string &what_arg)
GCRODRSolMgrRecyclingFailure is thrown when any problem occurs in using/creating the recycling subspa...
GCRODRSolMgrRecyclingFailure(const std::string &what_arg)
std::vector< ScalarType > work_
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > convTest_
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Get current linear problem being solved for in 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...
virtual ~GCRODRSolMgr()
Destructor.
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
std::vector< ScalarType > tau_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H_
MagnitudeType achievedTol_
Teuchos::RCP< OutputManager< ScalarType > > printer_
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
Teuchos::RCP< std::ostream > outputStream_
Teuchos::LAPACK< int, ScalarType > lapack
OrthoManagerFactory< ScalarType, MV, OP > ortho_factory_type
Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > ortho_
Orthogonalization manager.
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > PP_
int getNumIters() const override
Get the iteration count for the most recent call to solve().
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
Teuchos::ScalarTraits< MagnitudeType > MT
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
Teuchos::RCP< Teuchos::Time > timerSolve_
MultiVecTraits< ScalarType, MV > MVT
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > expConvTest_
Teuchos::RCP< Teuchos::ParameterList > params_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H2_
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.
Teuchos::ScalarTraits< ScalarType > SCT
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve.
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > R_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > HP_
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
Implementation of the GCRODR (Recycling GMRES) iterative linear solver.
static const bool requiresLapack
Details::SolverManagerRequiresLapack< ScalarType, MV, OP, requiresLapack > base_type
GCRODRSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
A linear system to solve, and its associated information.
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
Enumeration of all valid Belos (Mat)OrthoManager classes.
std::ostream & printValidNames(std::ostream &out) const
Print all recognized MatOrthoManager names to the given ostream.
bool isValidName(const std::string &name) const
Whether this factory recognizes the MatOrthoManager with the given name.
std::string validNamesString() const
List (as a string) of recognized MatOrthoManager names.
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.
Teuchos::RCP< const Teuchos::ParameterList > getDefaultParameters(const std::string &name) const
Default parameters for the given MatOrthoManager subclass.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
A class for extending the status testing capabilities of Belos via logical combinations.
Exception thrown to signal error in a status test during Belos::StatusTest::checkStatus().
An implementation of StatusTestResNorm using a family of residual norms.
int setTolerance(MagnitudeType tolerance)
Set the value of the tolerance.
A Belos::StatusTest class for specifying a maximum number of iterations.
A factory class for generating StatusTestOutput objects.
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.
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
ReturnType
Whether the Belos solve converged for all linear systems.
ScaleType
The type of scaling to use on the residual norm value.
ResetType
How to reset the solver.
static const double convTol
Default convergence tolerance.
Structure to contain pointers to GCRODRIter state variables.
Teuchos::RCP< MV > U
The recycled subspace and its projection.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B
The projection of the Krylov subspace against the recycled subspace.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
Teuchos::RCP< MV > V
The current Krylov basis.
int curDim
The current dimension of the reduction.