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>
183 #if defined(HAVE_TEUCHOSCORE_CXX11) 184 # if defined(HAVE_TEUCHOS_COMPLEX) 185 static_assert (std::is_same<ScalarType, std::complex<float> >::value ||
186 std::is_same<ScalarType, std::complex<double> >::value ||
187 std::is_same<ScalarType, float>::value ||
188 std::is_same<ScalarType, double>::value,
189 "Belos::GCRODRSolMgr: ScalarType must be one of the four " 190 "types (S,D,C,Z) supported by LAPACK.");
192 static_assert (std::is_same<ScalarType, float>::value ||
193 std::is_same<ScalarType, double>::value,
194 "Belos::GCRODRSolMgr: ScalarType must be float or double. " 195 "Complex arithmetic support is currently disabled. To " 196 "enable it, set Teuchos_ENABLE_COMPLEX=ON.");
197 # endif // defined(HAVE_TEUCHOS_COMPLEX) 198 #endif // defined(HAVE_TEUCHOSCORE_CXX11) 203 typedef Teuchos::ScalarTraits<ScalarType> SCT;
204 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
205 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
272 const Teuchos::RCP<Teuchos::ParameterList> &pl);
289 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters()
const;
302 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
303 return Teuchos::tuple(timerSolve_);
335 void setParameters(
const Teuchos::RCP<Teuchos::ParameterList> ¶ms );
347 bool set = problem_->setProblem();
349 throw "Could not set problem.";
393 std::string description()
const;
403 void initializeStateStorage();
412 int getHarmonicVecs1(
int m,
413 const Teuchos::SerialDenseMatrix<int,ScalarType>& HH,
414 Teuchos::SerialDenseMatrix<int,ScalarType>& PP);
421 int getHarmonicVecs2(
int keff,
int m,
422 const Teuchos::SerialDenseMatrix<int,ScalarType>& HH,
423 const Teuchos::RCP<const MV>& VV,
424 Teuchos::SerialDenseMatrix<int,ScalarType>& PP);
427 void sort(std::vector<MagnitudeType>& dlist,
int n, std::vector<int>& iperm);
430 Teuchos::LAPACK<int,ScalarType> lapack;
433 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
436 Teuchos::RCP<OutputManager<ScalarType> > printer_;
437 Teuchos::RCP<std::ostream> outputStream_;
440 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
441 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
442 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
443 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
444 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
449 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
452 Teuchos::RCP<Teuchos::ParameterList> params_;
455 static const MagnitudeType convTol_default_;
456 static const MagnitudeType orthoKappa_default_;
457 static const int maxRestarts_default_;
458 static const int maxIters_default_;
459 static const int numBlocks_default_;
460 static const int blockSize_default_;
461 static const int recycledBlocks_default_;
462 static const int verbosity_default_;
463 static const int outputStyle_default_;
464 static const int outputFreq_default_;
465 static const std::string impResScale_default_;
466 static const std::string expResScale_default_;
467 static const std::string label_default_;
468 static const std::string orthoType_default_;
469 static const Teuchos::RCP<std::ostream> outputStream_default_;
472 MagnitudeType convTol_, orthoKappa_, achievedTol_;
473 int maxRestarts_, maxIters_, numIters_;
474 int verbosity_, outputStyle_, outputFreq_;
475 std::string orthoType_;
476 std::string impResScale_, expResScale_;
483 int numBlocks_, recycledBlocks_;
494 Teuchos::RCP<MV> U_, C_;
497 Teuchos::RCP<MV> U1_, C1_;
500 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H2_;
501 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_;
502 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B_;
503 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PP_;
504 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > HP_;
505 std::vector<ScalarType> tau_;
506 std::vector<ScalarType> work_;
507 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > R_;
508 std::vector<int> ipiv_;
513 Teuchos::RCP<Teuchos::Time> timerSolve_;
519 bool builtRecycleSpace_;
524 template<
class ScalarType,
class MV,
class OP>
528 template<
class ScalarType,
class MV,
class OP>
532 template<
class ScalarType,
class MV,
class OP>
535 template<
class ScalarType,
class MV,
class OP>
538 template<
class ScalarType,
class MV,
class OP>
541 template<
class ScalarType,
class MV,
class OP>
544 template<
class ScalarType,
class MV,
class OP>
547 template<
class ScalarType,
class MV,
class OP>
550 template<
class ScalarType,
class MV,
class OP>
553 template<
class ScalarType,
class MV,
class OP>
556 template<
class ScalarType,
class MV,
class OP>
559 template<
class ScalarType,
class MV,
class OP>
562 template<
class ScalarType,
class MV,
class OP>
565 template<
class ScalarType,
class MV,
class OP>
568 template<
class ScalarType,
class MV,
class OP>
569 const Teuchos::RCP<std::ostream>
574 template<
class ScalarType,
class MV,
class OP>
584 template<
class ScalarType,
class MV,
class OP>
587 const Teuchos::RCP<Teuchos::ParameterList>& pl):
595 TEUCHOS_TEST_FOR_EXCEPTION(
596 problem == Teuchos::null, std::invalid_argument,
597 "Belos::GCRODRSolMgr constructor: The solver manager's " 598 "constructor needs the linear problem argument 'problem' " 607 if (! pl.is_null ()) {
613 template<
class ScalarType,
class MV,
class OP>
615 outputStream_ = outputStream_default_;
616 convTol_ = convTol_default_;
617 orthoKappa_ = orthoKappa_default_;
618 maxRestarts_ = maxRestarts_default_;
619 maxIters_ = maxIters_default_;
620 numBlocks_ = numBlocks_default_;
621 recycledBlocks_ = recycledBlocks_default_;
622 verbosity_ = verbosity_default_;
623 outputStyle_ = outputStyle_default_;
624 outputFreq_ = outputFreq_default_;
625 orthoType_ = orthoType_default_;
626 impResScale_ = impResScale_default_;
627 expResScale_ = expResScale_default_;
628 label_ = label_default_;
630 builtRecycleSpace_ =
false;
646 template<
class ScalarType,
class MV,
class OP>
651 using Teuchos::isParameterType;
652 using Teuchos::getParameter;
654 using Teuchos::ParameterList;
655 using Teuchos::parameterList;
658 using Teuchos::rcp_dynamic_cast;
659 using Teuchos::rcpFromRef;
660 using Teuchos::Exceptions::InvalidParameter;
661 using Teuchos::Exceptions::InvalidParameterName;
662 using Teuchos::Exceptions::InvalidParameterType;
666 RCP<const ParameterList> defaultParams = getValidParameters();
684 if (params_.is_null()) {
685 params_ = parameterList (*defaultParams);
693 if (params_ != params) {
699 params_ = parameterList (*params);
734 params_->validateParametersAndSetDefaults (*defaultParams);
738 if (params->isParameter (
"Maximum Restarts")) {
739 maxRestarts_ = params->get(
"Maximum Restarts", maxRestarts_default_);
742 params_->set (
"Maximum Restarts", maxRestarts_);
746 if (params->isParameter (
"Maximum Iterations")) {
747 maxIters_ = params->get (
"Maximum Iterations", maxIters_default_);
750 params_->set (
"Maximum Iterations", maxIters_);
751 if (! maxIterTest_.is_null())
752 maxIterTest_->setMaxIters (maxIters_);
756 if (params->isParameter (
"Num Blocks")) {
757 numBlocks_ = params->get (
"Num Blocks", numBlocks_default_);
758 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
759 "Belos::GCRODRSolMgr: The \"Num Blocks\" parameter must " 760 "be strictly positive, but you specified a value of " 761 << numBlocks_ <<
".");
763 params_->set (
"Num Blocks", numBlocks_);
767 if (params->isParameter (
"Num Recycled Blocks")) {
768 recycledBlocks_ = params->get (
"Num Recycled Blocks",
769 recycledBlocks_default_);
770 TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks_ <= 0, std::invalid_argument,
771 "Belos::GCRODRSolMgr: The \"Num Recycled Blocks\" " 772 "parameter must be strictly positive, but you specified " 773 "a value of " << recycledBlocks_ <<
".");
774 TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks_ >= numBlocks_, std::invalid_argument,
775 "Belos::GCRODRSolMgr: The \"Num Recycled Blocks\" " 776 "parameter must be less than the \"Num Blocks\" " 777 "parameter, but you specified \"Num Recycled Blocks\" " 778 "= " << recycledBlocks_ <<
" and \"Num Blocks\" = " 779 << numBlocks_ <<
".");
781 params_->set(
"Num Recycled Blocks", recycledBlocks_);
787 if (params->isParameter (
"Timer Label")) {
788 std::string tempLabel = params->get (
"Timer Label", label_default_);
791 if (tempLabel != label_) {
793 params_->set (
"Timer Label", label_);
794 std::string solveLabel = label_ +
": GCRODRSolMgr total solve time";
795 #ifdef BELOS_TEUCHOS_TIME_MONITOR 796 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
798 if (ortho_ != Teuchos::null) {
799 ortho_->setLabel( label_ );
805 if (params->isParameter (
"Verbosity")) {
806 if (isParameterType<int> (*params,
"Verbosity")) {
807 verbosity_ = params->get (
"Verbosity", verbosity_default_);
809 verbosity_ = (int) getParameter<Belos::MsgType> (*params,
"Verbosity");
812 params_->set (
"Verbosity", verbosity_);
815 if (! printer_.is_null())
816 printer_->setVerbosity (verbosity_);
820 if (params->isParameter (
"Output Style")) {
821 if (isParameterType<int> (*params,
"Output Style")) {
822 outputStyle_ = params->get (
"Output Style", outputStyle_default_);
824 outputStyle_ = (int) getParameter<OutputType> (*params,
"Output Style");
828 params_->set (
"Output Style", outputStyle_);
846 if (params->isParameter (
"Output Stream")) {
848 outputStream_ = getParameter<RCP<std::ostream> > (*params,
"Output Stream");
849 }
catch (InvalidParameter&) {
850 outputStream_ = rcpFromRef (std::cout);
857 if (outputStream_.is_null()) {
858 outputStream_ = rcp (
new Teuchos::oblackholestream);
861 params_->set (
"Output Stream", outputStream_);
864 if (! printer_.is_null()) {
865 printer_->setOStream (outputStream_);
871 if (params->isParameter (
"Output Frequency")) {
872 outputFreq_ = params->get (
"Output Frequency", outputFreq_default_);
876 params_->set(
"Output Frequency", outputFreq_);
877 if (! outputTest_.is_null())
878 outputTest_->setOutputFrequency (outputFreq_);
885 if (printer_.is_null()) {
896 bool changedOrthoType =
false;
897 if (params->isParameter (
"Orthogonalization")) {
898 const std::string& tempOrthoType =
899 params->get (
"Orthogonalization", orthoType_default_);
902 std::ostringstream os;
903 os <<
"Belos::GCRODRSolMgr: Invalid orthogonalization name \"" 904 << tempOrthoType <<
"\". The following are valid options " 905 <<
"for the \"Orthogonalization\" name parameter: ";
907 throw std::invalid_argument (os.str());
909 if (tempOrthoType != orthoType_) {
910 changedOrthoType =
true;
911 orthoType_ = tempOrthoType;
913 params_->set (
"Orthogonalization", orthoType_);
929 RCP<ParameterList> orthoParams;
932 using Teuchos::sublist;
934 const std::string paramName (
"Orthogonalization Parameters");
937 orthoParams = sublist (params_, paramName,
true);
938 }
catch (InvalidParameter&) {
945 orthoParams = sublist (params_, paramName,
true);
948 TEUCHOS_TEST_FOR_EXCEPTION(orthoParams.is_null(), std::logic_error,
949 "Failed to get orthogonalization parameters. " 950 "Please report this bug to the Belos developers.");
955 if (ortho_.is_null() || changedOrthoType) {
961 label_, orthoParams);
969 typedef Teuchos::ParameterListAcceptor PLA;
970 RCP<PLA> pla = rcp_dynamic_cast<PLA> (ortho_);
976 label_, orthoParams);
978 pla->setParameterList (orthoParams);
990 if (params->isParameter (
"Orthogonalization Constant")) {
991 const MagnitudeType orthoKappa =
992 params->get (
"Orthogonalization Constant", orthoKappa_default_);
993 if (orthoKappa > 0) {
994 orthoKappa_ = orthoKappa;
996 params_->set(
"Orthogonalization Constant", orthoKappa_);
998 if (orthoType_ ==
"DGKS" && ! ortho_.is_null()) {
1005 rcp_dynamic_cast<ortho_man_type>(ortho_)->setDepTol (orthoKappa_);
1015 if (params->isParameter(
"Convergence Tolerance")) {
1016 convTol_ = params->get (
"Convergence Tolerance", convTol_default_);
1019 params_->set (
"Convergence Tolerance", convTol_);
1020 if (! impConvTest_.is_null())
1022 if (! expConvTest_.is_null())
1023 expConvTest_->setTolerance (convTol_);
1027 if (params->isParameter (
"Implicit Residual Scaling")) {
1028 std::string tempImpResScale =
1029 getParameter<std::string> (*params,
"Implicit Residual Scaling");
1032 if (impResScale_ != tempImpResScale) {
1034 impResScale_ = tempImpResScale;
1037 params_->set(
"Implicit Residual Scaling", impResScale_);
1047 if (! impConvTest_.is_null()) {
1053 impConvTest_ = null;
1060 if (params->isParameter(
"Explicit Residual Scaling")) {
1061 std::string tempExpResScale =
1062 getParameter<std::string> (*params,
"Explicit Residual Scaling");
1065 if (expResScale_ != tempExpResScale) {
1067 expResScale_ = tempExpResScale;
1070 params_->set(
"Explicit Residual Scaling", expResScale_);
1073 if (! expConvTest_.is_null()) {
1079 expConvTest_ = null;
1090 if (maxIterTest_.is_null())
1095 if (impConvTest_.is_null()) {
1096 impConvTest_ = rcp (
new StatusTestResNorm_t (convTol_));
1102 if (expConvTest_.is_null()) {
1103 expConvTest_ = rcp (
new StatusTestResNorm_t (convTol_));
1104 expConvTest_->defineResForm (StatusTestResNorm_t::Explicit,
Belos::TwoNorm);
1110 if (convTest_.is_null()) {
1111 convTest_ = rcp (
new StatusTestCombo_t (StatusTestCombo_t::SEQ,
1119 sTest_ = rcp (
new StatusTestCombo_t (StatusTestCombo_t::OR,
1125 outputTest_ = stoFactory.
create (printer_, sTest_, outputFreq_,
1129 std::string solverDesc =
" GCRODR ";
1130 outputTest_->setSolverDesc( solverDesc );
1133 if (timerSolve_.is_null()) {
1134 std::string solveLabel = label_ +
": GCRODRSolMgr total solve time";
1135 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1136 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
1145 template<
class ScalarType,
class MV,
class OP>
1146 Teuchos::RCP<const Teuchos::ParameterList>
1149 using Teuchos::ParameterList;
1150 using Teuchos::parameterList;
1153 static RCP<const ParameterList> validPL;
1154 if (is_null(validPL)) {
1155 RCP<ParameterList> pl = parameterList ();
1158 pl->set(
"Convergence Tolerance", convTol_default_,
1159 "The relative residual tolerance that needs to be achieved by the\n" 1160 "iterative solver in order for the linear system to be declared converged.");
1161 pl->set(
"Maximum Restarts", maxRestarts_default_,
1162 "The maximum number of cycles allowed for each\n" 1163 "set of RHS solved.");
1164 pl->set(
"Maximum Iterations", maxIters_default_,
1165 "The maximum number of iterations allowed for each\n" 1166 "set of RHS solved.");
1170 pl->set(
"Block Size", blockSize_default_,
1171 "Block Size Parameter -- currently must be 1 for GCRODR");
1172 pl->set(
"Num Blocks", numBlocks_default_,
1173 "The maximum number of vectors allowed in the Krylov subspace\n" 1174 "for each set of RHS solved.");
1175 pl->set(
"Num Recycled Blocks", recycledBlocks_default_,
1176 "The maximum number of vectors in the recycled subspace." );
1177 pl->set(
"Verbosity", verbosity_default_,
1178 "What type(s) of solver information should be outputted\n" 1179 "to the output stream.");
1180 pl->set(
"Output Style", outputStyle_default_,
1181 "What style is used for the solver information outputted\n" 1182 "to the output stream.");
1183 pl->set(
"Output Frequency", outputFreq_default_,
1184 "How often convergence information should be outputted\n" 1185 "to the output stream.");
1186 pl->set(
"Output Stream", outputStream_default_,
1187 "A reference-counted pointer to the output stream where all\n" 1188 "solver output is sent.");
1189 pl->set(
"Implicit Residual Scaling", impResScale_default_,
1190 "The type of scaling used in the implicit residual convergence test.");
1191 pl->set(
"Explicit Residual Scaling", expResScale_default_,
1192 "The type of scaling used in the explicit residual convergence test.");
1193 pl->set(
"Timer Label", label_default_,
1194 "The string to use as a prefix for the timer labels.");
1198 pl->set(
"Orthogonalization", orthoType_default_,
1199 "The type of orthogonalization to use. Valid options: " +
1201 RCP<const ParameterList> orthoParams =
1203 pl->set (
"Orthogonalization Parameters", *orthoParams,
1204 "Parameters specific to the type of orthogonalization used.");
1206 pl->set(
"Orthogonalization Constant", orthoKappa_default_,
1207 "When using DGKS orthogonalization: the \"depTol\" constant, used " 1208 "to determine whether another step of classical Gram-Schmidt is " 1209 "necessary. Otherwise ignored.");
1216 template<
class ScalarType,
class MV,
class OP>
1219 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1222 Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
1223 if (rhsMV == Teuchos::null) {
1230 TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(numBlocks_) > MVT::GetGlobalLength(*rhsMV),std::invalid_argument,
1231 "Belos::GCRODRSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
1234 if (U_ == Teuchos::null) {
1235 U_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1239 if (MVT::GetNumberVecs(*U_) < recycledBlocks_+1) {
1240 Teuchos::RCP<const MV> tmp = U_;
1241 U_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1246 if (C_ == Teuchos::null) {
1247 C_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1251 if (MVT::GetNumberVecs(*C_) < recycledBlocks_+1) {
1252 Teuchos::RCP<const MV> tmp = C_;
1253 C_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1258 if (V_ == Teuchos::null) {
1259 V_ = MVT::Clone( *rhsMV, numBlocks_+1 );
1263 if (MVT::GetNumberVecs(*V_) < numBlocks_+1) {
1264 Teuchos::RCP<const MV> tmp = V_;
1265 V_ = MVT::Clone( *tmp, numBlocks_+1 );
1270 if (U1_ == Teuchos::null) {
1271 U1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1275 if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1276 Teuchos::RCP<const MV> tmp = U1_;
1277 U1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1282 if (C1_ == Teuchos::null) {
1283 C1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1287 if (MVT::GetNumberVecs(*C1_) < recycledBlocks_+1) {
1288 Teuchos::RCP<const MV> tmp = C1_;
1289 C1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1294 if (r_ == Teuchos::null)
1295 r_ = MVT::Clone( *rhsMV, 1 );
1298 tau_.resize(recycledBlocks_+1);
1301 work_.resize(recycledBlocks_+1);
1304 ipiv_.resize(recycledBlocks_+1);
1307 if (H2_ == Teuchos::null)
1308 H2_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 ) );
1310 if ( (H2_->numRows() != numBlocks_+recycledBlocks_+2) || (H2_->numCols() != numBlocks_+recycledBlocks_+1) )
1311 H2_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
1313 H2_->putScalar(zero);
1316 if (R_ == Teuchos::null)
1317 R_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycledBlocks_+1, recycledBlocks_+1 ) );
1319 if ( (R_->numRows() != recycledBlocks_+1) || (R_->numCols() != recycledBlocks_+1) )
1320 R_->reshape( recycledBlocks_+1, recycledBlocks_+1 );
1322 R_->putScalar(zero);
1325 if (PP_ == Teuchos::null)
1326 PP_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, recycledBlocks_+1 ) );
1328 if ( (PP_->numRows() != numBlocks_+recycledBlocks_+2) || (PP_->numCols() != recycledBlocks_+1) )
1329 PP_->reshape( numBlocks_+recycledBlocks_+2, recycledBlocks_+1 );
1333 if (HP_ == Teuchos::null)
1334 HP_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 ) );
1336 if ( (HP_->numRows() != numBlocks_+recycledBlocks_+2) || (HP_->numCols() != numBlocks_+recycledBlocks_+1) )
1337 HP_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
1345 template<
class ScalarType,
class MV,
class OP>
1353 if (!isSet_) { setParameters( params_ ); }
1355 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1356 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1357 std::vector<int> index(numBlocks_+1);
1359 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,
GCRODRSolMgrLinearProblemFailure,
"Belos::GCRODRSolMgr::solve(): Linear problem is not a valid object.");
1361 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),
GCRODRSolMgrLinearProblemFailure,
"Belos::GCRODRSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
1364 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1365 std::vector<int> currIdx(1);
1369 problem_->setLSIndex( currIdx );
1372 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1373 if (static_cast<ptrdiff_t>(numBlocks_) > dim) {
1374 numBlocks_ = Teuchos::as<int>(dim);
1376 "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
1377 " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
1378 params_->set(
"Num Blocks", numBlocks_);
1382 bool isConverged =
true;
1385 initializeStateStorage();
1389 Teuchos::ParameterList plist;
1391 plist.set(
"Num Blocks",numBlocks_);
1392 plist.set(
"Recycled Blocks",recycledBlocks_);
1397 RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_iter;
1400 int prime_iterations = 0;
1404 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1405 Teuchos::TimeMonitor slvtimer(*timerSolve_);
1408 while ( numRHS2Solve > 0 ) {
1411 builtRecycleSpace_ =
false;
1414 outputTest_->reset();
1422 "Belos::GCRODRSolMgr::solve(): Requested size of recycled subspace is not consistent with the current recycle subspace.");
1424 printer_->stream(
Debug) <<
" Now solving RHS index " << currIdx[0] <<
" using recycled subspace of dimension " << keff << std::endl << std::endl;
1427 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1428 RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1429 RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1430 problem_->apply( *Utmp, *Ctmp );
1432 RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1436 Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff, keff );
1437 int rank = ortho_->normalize(*Ctmp, rcp(&Rtmp,
false));
1439 TEUCHOS_TEST_FOR_EXCEPTION(rank != keff,
GCRODRSolMgrOrthoFailure,
"Belos::GCRODRSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
1444 ipiv_.resize(Rtmp.numRows());
1445 lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
1446 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1449 int lwork = Rtmp.numRows();
1450 work_.resize(lwork);
1451 lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1452 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
1455 MVT::MvTimesMatAddMv( one, *Utmp, Rtmp, zero, *U1tmp );
1460 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1461 Ctmp = MVT::CloneViewNonConst( *C_, index );
1462 Utmp = MVT::CloneView( *U_, index );
1465 Teuchos::SerialDenseMatrix<int,ScalarType> Ctr(keff,1);
1466 problem_->computeCurrPrecResVec( &*r_ );
1467 MVT::MvTransMv( one, *Ctmp, *r_, Ctr );
1470 RCP<MV> update = MVT::Clone( *problem_->getCurrLHSVec(), 1 );
1471 MVT::MvInit( *update, 0.0 );
1472 MVT::MvTimesMatAddMv( one, *Utmp, Ctr, one, *update );
1473 problem_->updateSolution( update,
true );
1476 MVT::MvTimesMatAddMv( -one, *Ctmp, Ctr, one, *r_ );
1479 prime_iterations = 0;
1485 printer_->stream(
Debug) <<
" No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
1487 Teuchos::ParameterList primeList;
1490 primeList.set(
"Num Blocks",numBlocks_);
1491 primeList.set(
"Recycled Blocks",0);
1494 RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_prime_iter;
1498 problem_->computeCurrPrecResVec( &*r_ );
1499 index.resize( 1 ); index[0] = 0;
1500 RCP<MV> v0 = MVT::CloneViewNonConst( *V_, index );
1501 MVT::SetBlock(*r_,index,*v0);
1505 index.resize( numBlocks_+1 );
1506 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1507 newstate.
V = MVT::CloneViewNonConst( *V_, index );
1508 newstate.
U = Teuchos::null;
1509 newstate.
C = Teuchos::null;
1510 newstate.
H = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, recycledBlocks_+1, recycledBlocks_+1 ) );
1511 newstate.
B = Teuchos::null;
1513 gcrodr_prime_iter->initialize(newstate);
1516 bool primeConverged =
false;
1518 gcrodr_prime_iter->iterate();
1521 if ( convTest_->getStatus() ==
Passed ) {
1523 primeConverged =
true;
1528 gcrodr_prime_iter->updateLSQR( gcrodr_prime_iter->getCurSubspaceDim() );
1531 sTest_->checkStatus( &*gcrodr_prime_iter );
1532 if (convTest_->getStatus() ==
Passed)
1533 primeConverged =
true;
1535 catch (
const std::exception &e) {
1536 printer_->stream(
Errors) <<
"Error! Caught exception in GCRODRIter::iterate() at iteration " 1537 << gcrodr_prime_iter->getNumIters() << std::endl
1538 << e.what() << std::endl;
1542 prime_iterations = gcrodr_prime_iter->getNumIters();
1545 RCP<MV> update = gcrodr_prime_iter->getCurrentUpdate();
1546 problem_->updateSolution( update,
true );
1549 newstate = gcrodr_prime_iter->getState();
1557 if (recycledBlocks_ < p+1) {
1559 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PPtmp = rcp (
new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, p, recycledBlocks_+1 ) );
1561 keff = getHarmonicVecs1( p, *newstate.
H, *PPtmp );
1563 PPtmp = rcp (
new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, p, keff ) );
1566 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1567 RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1568 RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1569 RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1571 for (
int ii=0; ii < p; ++ii) { index[ii] = ii; }
1572 RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1576 MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
1583 Teuchos::SerialDenseMatrix<int,ScalarType> Htmp( Teuchos::View, *H2_, p+1, p, recycledBlocks_+1,recycledBlocks_+1);
1584 Teuchos::SerialDenseMatrix<int,ScalarType> HPtmp( Teuchos::View, *HP_, p+1, keff );
1585 HPtmp.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, Htmp, *PPtmp, zero );
1592 lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
1593 HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
1594 TEUCHOS_TEST_FOR_EXCEPTION(
1596 " LAPACK's _GEQRF failed to compute a workspace size.");
1604 lwork = std::abs (static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work_[0])));
1605 work_.resize (lwork);
1606 lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
1607 HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
1608 TEUCHOS_TEST_FOR_EXCEPTION(
1610 " LAPACK's _GEQRF failed to compute a QR factorization.");
1614 Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff, keff );
1615 for (
int ii = 0; ii < keff; ++ii) {
1616 for (
int jj = ii; jj < keff; ++jj) {
1617 Rtmp(ii,jj) = HPtmp(ii,jj);
1624 lapack.UNGQR (HPtmp.numRows (), HPtmp.numCols (), HPtmp.numCols (),
1625 HPtmp.values (), HPtmp.stride (), &tau_[0], &work_[0],
1627 TEUCHOS_TEST_FOR_EXCEPTION(
1629 "LAPACK's _UNGQR failed to construct the Q factor.");
1634 index.resize (p + 1);
1635 for (
int ii = 0; ii < (p+1); ++ii) {
1638 Vtmp = MVT::CloneView( *V_, index );
1639 MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
1646 ipiv_.resize(Rtmp.numRows());
1647 lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
1648 TEUCHOS_TEST_FOR_EXCEPTION(
1650 "LAPACK's _GETRF failed to compute an LU factorization.");
1659 lwork = Rtmp.numRows();
1660 work_.resize(lwork);
1661 lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1662 TEUCHOS_TEST_FOR_EXCEPTION(
1664 "LAPACK's _GETRI failed to invert triangular matrix.");
1667 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
1669 printer_->stream(
Debug)
1670 <<
" Generated recycled subspace using RHS index " << currIdx[0]
1671 <<
" of dimension " << keff << std::endl << std::endl;
1676 if (primeConverged) {
1678 problem_->setCurrLS();
1682 if (numRHS2Solve > 0) {
1684 problem_->setLSIndex (currIdx);
1687 currIdx.resize (numRHS2Solve);
1697 gcrodr_iter->setSize( keff, numBlocks_ );
1700 gcrodr_iter->resetNumIters(prime_iterations);
1703 outputTest_->resetNumCalls();
1706 problem_->computeCurrPrecResVec( &*r_ );
1707 index.resize( 1 ); index[0] = 0;
1708 RCP<MV> v0 = MVT::CloneViewNonConst( *V_, index );
1709 MVT::SetBlock(*r_,index,*v0);
1713 index.resize( numBlocks_+1 );
1714 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1715 newstate.
V = MVT::CloneViewNonConst( *V_, index );
1716 index.resize( keff );
1717 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1718 newstate.
C = MVT::CloneViewNonConst( *C_, index );
1719 newstate.
U = MVT::CloneViewNonConst( *U_, index );
1720 newstate.
B = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, keff, numBlocks_, 0, keff ) );
1721 newstate.
H = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, keff, keff ) );
1723 gcrodr_iter->initialize(newstate);
1726 int numRestarts = 0;
1731 gcrodr_iter->iterate();
1738 if ( convTest_->getStatus() ==
Passed ) {
1747 else if ( maxIterTest_->getStatus() ==
Passed ) {
1749 isConverged =
false;
1757 else if ( gcrodr_iter->getCurSubspaceDim() == gcrodr_iter->getMaxSubspaceDim() ) {
1762 RCP<MV> update = gcrodr_iter->getCurrentUpdate();
1763 problem_->updateSolution( update,
true );
1765 buildRecycleSpace2(gcrodr_iter);
1767 printer_->stream(
Debug)
1768 <<
" Generated new recycled subspace using RHS index " 1769 << currIdx[0] <<
" of dimension " << keff << std::endl
1773 if (numRestarts >= maxRestarts_) {
1774 isConverged =
false;
1779 printer_->stream(
Debug)
1780 <<
" Performing restart number " << numRestarts <<
" of " 1781 << maxRestarts_ << std::endl << std::endl;
1784 problem_->computeCurrPrecResVec( &*r_ );
1785 index.resize( 1 ); index[0] = 0;
1786 RCP<MV> v00 = MVT::CloneViewNonConst( *V_, index );
1787 MVT::SetBlock(*r_,index,*v00);
1791 index.resize( numBlocks_+1 );
1792 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1793 restartState.
V = MVT::CloneViewNonConst( *V_, index );
1794 index.resize( keff );
1795 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1796 restartState.
U = MVT::CloneViewNonConst( *U_, index );
1797 restartState.
C = MVT::CloneViewNonConst( *C_, index );
1798 restartState.
B = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, keff, numBlocks_, 0, keff ) );
1799 restartState.
H = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, keff, keff ) );
1801 gcrodr_iter->initialize(restartState);
1814 TEUCHOS_TEST_FOR_EXCEPTION(
1815 true, std::logic_error,
"Belos::GCRODRSolMgr::solve: " 1816 "Invalid return from GCRODRIter::iterate().");
1821 gcrodr_iter->updateLSQR( gcrodr_iter->getCurSubspaceDim() );
1824 sTest_->checkStatus( &*gcrodr_iter );
1825 if (convTest_->getStatus() !=
Passed)
1826 isConverged =
false;
1829 catch (
const std::exception& e) {
1831 <<
"Error! Caught exception in GCRODRIter::iterate() at iteration " 1832 << gcrodr_iter->getNumIters() << std::endl << e.what() << std::endl;
1839 RCP<MV> update = gcrodr_iter->getCurrentUpdate();
1840 problem_->updateSolution( update,
true );
1843 problem_->setCurrLS();
1848 if (!builtRecycleSpace_) {
1849 buildRecycleSpace2(gcrodr_iter);
1850 printer_->stream(
Debug)
1851 <<
" Generated new recycled subspace using RHS index " << currIdx[0]
1852 <<
" of dimension " << keff << std::endl << std::endl;
1857 if (numRHS2Solve > 0) {
1859 problem_->setLSIndex (currIdx);
1862 currIdx.resize (numRHS2Solve);
1870 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1875 Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
1876 #endif // BELOS_TEUCHOS_TIME_MONITOR 1879 numIters_ = maxIterTest_->getNumIters ();
1891 const std::vector<MagnitudeType>* pTestValues = expConvTest_->getTestValue();
1892 if (pTestValues == NULL || pTestValues->size() < 1) {
1893 pTestValues = impConvTest_->getTestValue();
1895 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1896 "Belos::GCRODRSolMgr::solve(): The implicit convergence test's getTestValue() " 1897 "method returned NULL. Please report this bug to the Belos developers.");
1898 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1899 "Belos::GCRODRSolMgr::solve(): The implicit convergence test's getTestValue() " 1900 "method returned a vector of length zero. Please report this bug to the " 1901 "Belos developers.");
1906 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1913 template<
class ScalarType,
class MV,
class OP>
1916 MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
1917 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1919 std::vector<MagnitudeType> d(keff);
1920 std::vector<ScalarType> dscalar(keff);
1921 std::vector<int> index(numBlocks_+1);
1933 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1934 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1936 dscalar.resize(keff);
1937 MVT::MvNorm( *Utmp, d );
1938 for (
int i=0; i<keff; ++i) {
1940 dscalar[i] = (ScalarType)d[i];
1942 MVT::MvScale( *Utmp, dscalar );
1946 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H2tmp = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, p+keff+1, p+keff ) );
1949 for (
int i=0; i<keff; ++i) {
1950 (*H2tmp)(i,i) = d[i];
1957 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p+keff, recycledBlocks_+1 );
1958 keff_new = getHarmonicVecs2( keff, p, *H2tmp, oldState.
V, PPtmp );
1965 Teuchos::RCP<MV> U1tmp;
1967 index.resize( keff );
1968 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1969 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1970 index.resize( keff_new );
1971 for (
int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
1972 U1tmp = MVT::CloneViewNonConst( *U1_, index );
1973 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, keff, keff_new );
1974 MVT::MvTimesMatAddMv( one, *Utmp, PPtmp, zero, *U1tmp );
1980 for (
int ii=0; ii < p; ii++) { index[ii] = ii; }
1981 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1982 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p, keff_new, keff );
1983 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
1987 Teuchos::SerialDenseMatrix<int,ScalarType> HPtmp( Teuchos::View, *HP_, p+keff+1, keff_new );
1989 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p+keff, keff_new );
1990 HPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*H2tmp,PPtmp,zero);
1994 int info = 0, lwork = -1;
1995 tau_.resize (keff_new);
1996 lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
1997 HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
1998 TEUCHOS_TEST_FOR_EXCEPTION(
2000 "LAPACK's _GEQRF failed to compute a workspace size.");
2006 lwork = std::abs (static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work_[0])));
2007 work_.resize (lwork);
2008 lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
2009 HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
2010 TEUCHOS_TEST_FOR_EXCEPTION(
2012 "LAPACK's _GEQRF failed to compute a QR factorization.");
2016 Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff_new, keff_new );
2017 for(
int i=0;i<keff_new;i++) {
for(
int j=i;j<keff_new;j++) Rtmp(i,j) = HPtmp(i,j); }
2023 lapack.UNGQR (HPtmp.numRows (), HPtmp.numCols (), HPtmp.numCols (),
2024 HPtmp.values (), HPtmp.stride (), &tau_[0], &work_[0],
2026 TEUCHOS_TEST_FOR_EXCEPTION(
2028 "LAPACK's _UNGQR failed to construct the Q factor.");
2035 Teuchos::RCP<MV> C1tmp;
2038 for (
int i=0; i < keff; i++) { index[i] = i; }
2039 Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
2040 index.resize(keff_new);
2041 for (
int i=0; i < keff_new; i++) { index[i] = i; }
2042 C1tmp = MVT::CloneViewNonConst( *C1_, index );
2043 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *HP_, keff, keff_new );
2044 MVT::MvTimesMatAddMv( one, *Ctmp, PPtmp, zero, *C1tmp );
2048 index.resize( p+1 );
2049 for (
int i=0; i < p+1; ++i) { index[i] = i; }
2050 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
2051 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *HP_, p+1, keff_new, keff, 0 );
2052 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *C1tmp );
2061 ipiv_.resize(Rtmp.numRows());
2062 lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
2063 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
2066 lwork = Rtmp.numRows();
2067 work_.resize(lwork);
2068 lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
2069 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to compute an LU factorization.");
2072 index.resize(keff_new);
2073 for (
int i=0; i < keff_new; i++) { index[i] = i; }
2074 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
2075 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
2079 if (keff != keff_new) {
2081 gcrodr_iter->setSize( keff, numBlocks_ );
2083 Teuchos::SerialDenseMatrix<int,ScalarType> b1( Teuchos::View, *H2_, recycledBlocks_+2, 1, 0, recycledBlocks_ );
2091 template<
class ScalarType,
class MV,
class OP>
2093 const Teuchos::SerialDenseMatrix<int,ScalarType>& HH,
2094 Teuchos::SerialDenseMatrix<int,ScalarType>& PP) {
2096 bool xtraVec =
false;
2097 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
2101 std::vector<MagnitudeType> wr(m), wi(m);
2104 Teuchos::SerialDenseMatrix<int,ScalarType> vr(m,m,
false);
2107 std::vector<MagnitudeType> w(m);
2110 std::vector<int> iperm(m);
2114 std::vector<ScalarType> work(lwork);
2115 std::vector<MagnitudeType> rwork(lwork);
2121 builtRecycleSpace_ =
true;
2124 Teuchos::SerialDenseMatrix<int, ScalarType> HHt( HH,
Teuchos::TRANS );
2125 Teuchos::SerialDenseVector<int, ScalarType> e_m( m );
2127 lapack.GESV(m, 1, HHt.values(), HHt.stride(), &iperm[0], e_m.values(), e_m.stride(), &info);
2128 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK GESV failed to compute a solution.");
2131 ScalarType d = HH(m, m-1) * HH(m, m-1);
2132 Teuchos::SerialDenseMatrix<int, ScalarType> harmHH( Teuchos::Copy, HH, HH.numRows(), HH.numCols() );
2133 for( i=0; i<m; ++i )
2134 harmHH(i, m-1) += d * e_m[i];
2142 lapack.GEEV(
'N',
'V', m, harmHH.values(), harmHH.stride(), &wr[0], &wi[0],
2143 vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
2144 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK GEEV failed to compute eigensolutions.");
2147 for( i=0; i<m; ++i )
2148 w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] );
2151 this->sort(w, m, iperm);
2153 const bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
2156 for( i=0; i<recycledBlocks_; ++i ) {
2157 for( j=0; j<m; j++ ) {
2158 PP(j,i) = vr(j,iperm[i]);
2162 if(!scalarTypeIsComplex) {
2165 if (wi[iperm[recycledBlocks_-1]] != 0.0) {
2167 for ( i=0; i<recycledBlocks_; ++i ) {
2168 if (wi[iperm[i]] != 0.0)
2177 if (wi[iperm[recycledBlocks_-1]] > 0.0) {
2178 for( j=0; j<m; ++j ) {
2179 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
2183 for( j=0; j<m; ++j ) {
2184 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
2193 return recycledBlocks_+1;
2196 return recycledBlocks_;
2202 template<
class ScalarType,
class MV,
class OP>
2204 const Teuchos::SerialDenseMatrix<int,ScalarType>& HH,
2205 const Teuchos::RCP<const MV>& VV,
2206 Teuchos::SerialDenseMatrix<int,ScalarType>& PP) {
2208 int m2 = HH.numCols();
2209 bool xtraVec =
false;
2210 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
2211 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
2212 std::vector<int> index;
2215 std::vector<MagnitudeType> wr(m2), wi(m2);
2218 std::vector<MagnitudeType> w(m2);
2221 Teuchos::SerialDenseMatrix<int,ScalarType> vr(m2,m2,
false);
2224 std::vector<int> iperm(m2);
2227 builtRecycleSpace_ =
true;
2232 Teuchos::SerialDenseMatrix<int,ScalarType> B(m2,m2,
false);
2237 Teuchos::SerialDenseMatrix<int,ScalarType> A_tmp( keffloc+m+1, keffloc+m );
2240 index.resize(keffloc);
2241 for (i=0; i<keffloc; ++i) { index[i] = i; }
2242 Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
2243 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
2244 Teuchos::SerialDenseMatrix<int,ScalarType> A11( Teuchos::View, A_tmp, keffloc, keffloc );
2245 MVT::MvTransMv( one, *Ctmp, *Utmp, A11 );
2248 Teuchos::SerialDenseMatrix<int,ScalarType> A21( Teuchos::View, A_tmp, m+1, keffloc, keffloc );
2250 for (i=0; i < m+1; i++) { index[i] = i; }
2251 Teuchos::RCP<const MV> Vp = MVT::CloneView( *VV, index );
2252 MVT::MvTransMv( one, *Vp, *Utmp, A21 );
2255 for( i=keffloc; i<keffloc+m; i++ ) {
2260 Teuchos::SerialDenseMatrix<int,ScalarType> A( m2, A_tmp.numCols() );
2261 A.multiply(
Teuchos::TRANS, Teuchos::NO_TRANS, one, HH, A_tmp, zero );
2269 char balanc=
'P', jobvl=
'N', jobvr=
'V', sense=
'N';
2270 int ld = A.numRows();
2272 int ldvl = ld, ldvr = ld;
2273 int info = 0,ilo = 0,ihi = 0;
2274 MagnitudeType abnrm = 0.0, bbnrm = 0.0;
2276 std::vector<ScalarType> beta(ld);
2277 std::vector<ScalarType> work(lwork);
2278 std::vector<MagnitudeType> rwork(lwork);
2279 std::vector<MagnitudeType> lscale(ld), rscale(ld);
2280 std::vector<MagnitudeType> rconde(ld), rcondv(ld);
2281 std::vector<int> iwork(ld+6);
2286 lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0],
2287 &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
2288 &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &rwork[0],
2289 &iwork[0], bwork, &info);
2290 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK GGEVX failed to compute eigensolutions.");
2294 for( i=0; i<ld; i++ ) {
2295 w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot (wr[i]*wr[i] + wi[i]*wi[i]) /
2296 Teuchos::ScalarTraits<ScalarType>::magnitude (beta[i]);
2300 this->sort(w,ld,iperm);
2302 const bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
2305 for( i=0; i<recycledBlocks_; i++ ) {
2306 for( j=0; j<ld; j++ ) {
2307 PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
2311 if(!scalarTypeIsComplex) {
2314 if (wi[iperm[ld-recycledBlocks_]] != 0.0) {
2316 for ( i=ld-recycledBlocks_; i<ld; i++ ) {
2317 if (wi[iperm[i]] != 0.0)
2326 if (wi[iperm[ld-recycledBlocks_]] > 0.0) {
2327 for( j=0; j<ld; j++ ) {
2328 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
2332 for( j=0; j<ld; j++ ) {
2333 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
2342 return recycledBlocks_+1;
2345 return recycledBlocks_;
2352 template<
class ScalarType,
class MV,
class OP>
2354 int l, r, j, i, flag;
2356 MagnitudeType dRR, dK;
2383 if (dlist[j] > dlist[j - 1]) j = j + 1;
2385 if (dlist[j - 1] > dK) {
2386 dlist[i - 1] = dlist[j - 1];
2387 iperm[i - 1] = iperm[j - 1];
2401 dlist[r] = dlist[0];
2402 iperm[r] = iperm[0];
2417 template<
class ScalarType,
class MV,
class OP>
2419 std::ostringstream out;
2420 out <<
"Belos::GCRODRSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<
">";
2422 out <<
"Ortho Type: \"" << orthoType_ <<
"\"";
2423 out <<
", Num Blocks: " <<numBlocks_;
2424 out <<
", Num Recycle Blocks: " << recycledBlocks_;
2425 out <<
", Max Restarts: " << maxRestarts_;
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Collection of types and exceptions used within the Belos solvers.
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.
bool isLOADetected() const
Return whether a loss of accuracy was detected by this solver during the most current solve...
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.
Teuchos::RCP< const Teuchos::ParameterList > getDefaultParameters(const std::string &name) const
Default parameters for the given MatOrthoManager subclass.
void reset(const ResetType type)
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
An implementation of StatusTestResNorm using a family of residual norms.
GCRODRSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i.e.
GCRODRSolMgrRecyclingFailure(const std::string &what_arg)
This class implements the GCRODR iteration, where a single-std::vector Krylov subspace is constructed...
Belos::StatusTest class for specifying a maximum number of iterations.
int curDim
The current dimension of the reduction.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get current linear problem being solved for in this object.
A factory class for generating StatusTestOutput objects.
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 > > &outMan, const std::string &label, const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Return an instance of the specified MatOrthoManager subclass.
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
A Belos::StatusTest class for specifying a maximum number of iterations.
GCRODRSolMgrOrthoFailure(const std::string &what_arg)
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.
int getNumIters() const
Get the iteration count for the most recent call to solve().
bool isValidName(const std::string &name) const
Whether this factory recognizes the MatOrthoManager with the given name.
Teuchos::RCP< MV > U
The recycled subspace and its projection.
std::ostream & printValidNames(std::ostream &out) const
Print all recognized MatOrthoManager names to the given ostream.
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...
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.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
MagnitudeType achievedTol() const
Tolerance achieved by the last solve() invocation.
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.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem)
Set the linear problem that needs to be solved.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B
The projection of the Krylov subspace against the recycled subspace.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
A class for extending the status testing capabilities of Belos via logical combinations.
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...
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
Belos concrete class for performing the block, flexible GMRES iteration.