Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosPseudoBlockGmresSolMgr.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Belos: Block Linear Solvers Package
4 //
5 // Copyright 2004-2016 NTESS and the Belos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef BELOS_PSEUDO_BLOCK_GMRES_SOLMGR_HPP
11 #define BELOS_PSEUDO_BLOCK_GMRES_SOLMGR_HPP
12 
17 #include "BelosConfigDefs.hpp"
18 #include "BelosTypes.hpp"
19 
20 #include "BelosLinearProblem.hpp"
21 #include "BelosSolverManager.hpp"
22 
27 #include "BelosOutputManager.hpp"
28 #ifdef BELOS_TEUCHOS_TIME_MONITOR
29 #include "Teuchos_TimeMonitor.hpp"
30 #endif
31 
42 namespace Belos {
43 
45 
46 
54  PseudoBlockGmresSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
55  {}};
56 
64  PseudoBlockGmresSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
65  {}};
66 
90  template<class ScalarType, class MV, class OP>
91  class PseudoBlockGmresSolMgr : public SolverManager<ScalarType,MV,OP> {
92 
93  private:
97  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
99 
100  public:
101 
103 
104 
113 
226 
229 
233  }
235 
237 
238 
239  const LinearProblem<ScalarType,MV,OP>& getProblem() const override {
240  return *problem_;
241  }
242 
245 
248 
255  return Teuchos::tuple(timerSolve_);
256  }
257 
268  MagnitudeType achievedTol() const override {
269  return achievedTol_;
270  }
271 
273  int getNumIters() const override {
274  return numIters_;
275  }
276 
332  bool isLOADetected() const override { return loaDetected_; }
333 
336  getResidualStatusTest() const { return impConvTest_.getRawPtr(); }
338 
340 
341 
343  void setProblem (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem) override {
344  problem_ = problem;
345  }
346 
348  void setParameters (const Teuchos::RCP<Teuchos::ParameterList> &params) override;
349 
356  virtual void setUserConvStatusTest(
357  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &userConvStatusTest,
358  const typename StatusTestCombo<ScalarType,MV,OP>::ComboType &comboType =
360  ) override;
361 
363  void setDebugStatusTest( const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &debugStatusTest ) override;
364 
366 
368 
369 
373  void reset( const ResetType type ) override { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
375 
377 
378 
396  ReturnType solve() override;
397 
399 
402 
409  void
411  const Teuchos::EVerbosityLevel verbLevel =
413 
415  std::string description () const override;
416 
418 
419  private:
420 
435  bool checkStatusTest();
436 
439 
440  // Output manager.
442  Teuchos::RCP<std::ostream> outputStream_;
443 
444  // Status tests.
445  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > userConvStatusTest_;
450  Teuchos::RCP<StatusTestResNorm<ScalarType,MV,OP> > impConvTest_, expConvTest_;
454 
455  // Orthogonalization manager.
457 
458  // Current parameter list.
460 
461  // Default solver values.
462  static constexpr int maxRestarts_default_ = 20;
463  static constexpr int maxIters_default_ = 1000;
464  static constexpr bool showMaxResNormOnly_default_ = false;
465  static constexpr int blockSize_default_ = 1;
466  static constexpr int numBlocks_default_ = 300;
467  static constexpr int verbosity_default_ = Belos::Errors;
468  static constexpr int outputStyle_default_ = Belos::General;
469  static constexpr int outputFreq_default_ = -1;
470  static constexpr int defQuorum_default_ = 1;
471  static constexpr const char * impResScale_default_ = "Norm of Preconditioned Initial Residual";
472  static constexpr const char * expResScale_default_ = "Norm of Initial Residual";
473  static constexpr const char * label_default_ = "Belos";
474  static constexpr const char * orthoType_default_ = "ICGS";
475 
476  // Current solver values.
477  MagnitudeType convtol_, orthoKappa_, achievedTol_;
478  int maxRestarts_, maxIters_, numIters_;
479  int blockSize_, numBlocks_, verbosity_, outputStyle_, outputFreq_, defQuorum_;
480  bool showMaxResNormOnly_;
481  std::string orthoType_;
482  std::string impResScale_, expResScale_;
483  MagnitudeType resScaleFactor_;
484 
485  // Timers.
486  std::string label_;
487  Teuchos::RCP<Teuchos::Time> timerSolve_;
488 
489  // Internal state variables.
490  bool isSet_, isSTSet_, expResTest_;
491  bool loaDetected_;
492  };
493 
494 
495 // Empty Constructor
496 template<class ScalarType, class MV, class OP>
498  outputStream_(Teuchos::rcpFromRef(std::cout)),
499  taggedTests_(Teuchos::null),
500  convtol_(DefaultSolverParameters::convTol),
501  orthoKappa_(DefaultSolverParameters::orthoKappa),
502  achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
503  maxRestarts_(maxRestarts_default_),
504  maxIters_(maxIters_default_),
505  numIters_(0),
506  blockSize_(blockSize_default_),
507  numBlocks_(numBlocks_default_),
508  verbosity_(verbosity_default_),
509  outputStyle_(outputStyle_default_),
510  outputFreq_(outputFreq_default_),
511  defQuorum_(defQuorum_default_),
512  showMaxResNormOnly_(showMaxResNormOnly_default_),
513  orthoType_(orthoType_default_),
514  impResScale_(impResScale_default_),
515  expResScale_(expResScale_default_),
516  resScaleFactor_(DefaultSolverParameters::resScaleFactor),
517  label_(label_default_),
518  isSet_(false),
519  isSTSet_(false),
520  expResTest_(false),
521  loaDetected_(false)
522 {}
523 
524 // Basic Constructor
525 template<class ScalarType, class MV, class OP>
529  problem_(problem),
530  outputStream_(Teuchos::rcpFromRef(std::cout)),
531  taggedTests_(Teuchos::null),
532  convtol_(DefaultSolverParameters::convTol),
533  orthoKappa_(DefaultSolverParameters::orthoKappa),
534  achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
535  maxRestarts_(maxRestarts_default_),
536  maxIters_(maxIters_default_),
537  numIters_(0),
538  blockSize_(blockSize_default_),
539  numBlocks_(numBlocks_default_),
540  verbosity_(verbosity_default_),
541  outputStyle_(outputStyle_default_),
542  outputFreq_(outputFreq_default_),
543  defQuorum_(defQuorum_default_),
544  showMaxResNormOnly_(showMaxResNormOnly_default_),
545  orthoType_(orthoType_default_),
546  impResScale_(impResScale_default_),
547  expResScale_(expResScale_default_),
548  resScaleFactor_(DefaultSolverParameters::resScaleFactor),
549  label_(label_default_),
550  isSet_(false),
551  isSTSet_(false),
552  expResTest_(false),
553  loaDetected_(false)
554 {
555  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
556 
557  // If the parameter list pointer is null, then set the current parameters to the default parameter list.
558  if (!is_null(pl)) {
559  // Set the parameters using the list that was passed in.
560  setParameters( pl );
561  }
562 }
563 
564 template<class ScalarType, class MV, class OP>
565 void
568 {
570  using Teuchos::parameterList;
571  using Teuchos::rcp;
572  using Teuchos::rcp_dynamic_cast;
573 
574  // Create the internal parameter list if one doesn't already exist.
575  if (params_ == Teuchos::null) {
576  params_ = parameterList (*getValidParameters ());
577  } else {
578  // TAW: 3/8/2016: do not validate sub parameter lists as they
579  // might not have a pre-defined structure
580  // e.g. user-specified status tests
581  // The Belos Pseudo Block GMRES parameters on the first level are
582  // not affected and verified.
583  params->validateParameters (*getValidParameters (), 0);
584  }
585 
586  // Check for maximum number of restarts
587  if (params->isParameter ("Maximum Restarts")) {
588  maxRestarts_ = params->get ("Maximum Restarts", maxRestarts_default_);
589 
590  // Update parameter in our list.
591  params_->set ("Maximum Restarts", maxRestarts_);
592  }
593 
594  // Check for maximum number of iterations
595  if (params->isParameter ("Maximum Iterations")) {
596  maxIters_ = params->get ("Maximum Iterations", maxIters_default_);
597 
598  // Update parameter in our list and in status test.
599  params_->set ("Maximum Iterations", maxIters_);
600  if (! maxIterTest_.is_null ()) {
601  maxIterTest_->setMaxIters (maxIters_);
602  }
603  }
604 
605  // Check for blocksize
606  if (params->isParameter ("Block Size")) {
607  blockSize_ = params->get ("Block Size", blockSize_default_);
609  blockSize_ <= 0, std::invalid_argument,
610  "Belos::PseudoBlockGmresSolMgr::setParameters: "
611  "The \"Block Size\" parameter must be strictly positive, "
612  "but you specified a value of " << blockSize_ << ".");
613 
614  // Update parameter in our list.
615  params_->set ("Block Size", blockSize_);
616  }
617 
618  // Check for the maximum number of blocks.
619  if (params->isParameter ("Num Blocks")) {
620  numBlocks_ = params->get ("Num Blocks", numBlocks_default_);
622  numBlocks_ <= 0, std::invalid_argument,
623  "Belos::PseudoBlockGmresSolMgr::setParameters: "
624  "The \"Num Blocks\" parameter must be strictly positive, "
625  "but you specified a value of " << numBlocks_ << ".");
626 
627  // Update parameter in our list.
628  params_->set ("Num Blocks", numBlocks_);
629  }
630 
631  // Check to see if the timer label changed.
632  if (params->isParameter ("Timer Label")) {
633  const std::string tempLabel = params->get ("Timer Label", label_default_);
634 
635  // Update parameter in our list and solver timer
636  if (tempLabel != label_) {
637  label_ = tempLabel;
638  params_->set ("Timer Label", label_);
639  const std::string solveLabel =
640  label_ + ": PseudoBlockGmresSolMgr total solve time";
641 #ifdef BELOS_TEUCHOS_TIME_MONITOR
642  timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
643 #endif // BELOS_TEUCHOS_TIME_MONITOR
644  if (ortho_ != Teuchos::null) {
645  ortho_->setLabel( label_ );
646  }
647  }
648  }
649 
650 
651  // Check for a change in verbosity level
652  if (params->isParameter ("Verbosity")) {
653  if (Teuchos::isParameterType<int> (*params, "Verbosity")) {
654  verbosity_ = params->get ("Verbosity", verbosity_default_);
655  } else {
656  verbosity_ = (int) Teuchos::getParameter<Belos::MsgType> (*params, "Verbosity");
657  }
658 
659  // Update parameter in our list.
660  params_->set ("Verbosity", verbosity_);
661  if (! printer_.is_null ()) {
662  printer_->setVerbosity (verbosity_);
663  }
664  }
665 
666  // Check for a change in output style.
667  if (params->isParameter ("Output Style")) {
668  if (Teuchos::isParameterType<int> (*params, "Output Style")) {
669  outputStyle_ = params->get ("Output Style", outputStyle_default_);
670  } else {
671  outputStyle_ = (int) Teuchos::getParameter<Belos::OutputType> (*params, "Output Style");
672  }
673 
674  // Update parameter in our list.
675  params_->set ("Output Style", outputStyle_);
676  if (! outputTest_.is_null ()) {
677  isSTSet_ = false;
678  }
679 
680  }
681 
682  // Check if user has specified his own status tests
683  if (params->isSublist ("User Status Tests")) {
684  Teuchos::ParameterList userStatusTestsList = params->sublist("User Status Tests",true);
685  if ( userStatusTestsList.numParams() > 0 ) {
686  std::string userCombo_string = params->get<std::string>("User Status Tests Combo Type", "SEQ");
688  setUserConvStatusTest( testFactory->buildStatusTests(userStatusTestsList), testFactory->stringToComboType(userCombo_string) );
689  taggedTests_ = testFactory->getTaggedTests();
690  isSTSet_ = false;
691  }
692  }
693 
694  // output stream
695  if (params->isParameter ("Output Stream")) {
696  outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> > (*params, "Output Stream");
697 
698  // Update parameter in our list.
699  params_->set("Output Stream", outputStream_);
700  if (! printer_.is_null ()) {
701  printer_->setOStream (outputStream_);
702  }
703  }
704 
705  // frequency level
706  if (verbosity_ & Belos::StatusTestDetails) {
707  if (params->isParameter ("Output Frequency")) {
708  outputFreq_ = params->get ("Output Frequency", outputFreq_default_);
709  }
710 
711  // Update parameter in out list and output status test.
712  params_->set ("Output Frequency", outputFreq_);
713  if (! outputTest_.is_null ()) {
714  outputTest_->setOutputFrequency (outputFreq_);
715  }
716  }
717 
718  // Create output manager if we need to.
719  if (printer_.is_null ()) {
720  printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
721  }
722 
723  // Check if the orthogonalization changed.
724  bool changedOrthoType = false;
725  if (params->isParameter ("Orthogonalization")) {
726  std::string tempOrthoType = params->get ("Orthogonalization", orthoType_default_);
727  if (tempOrthoType != orthoType_) {
728  orthoType_ = tempOrthoType;
729  changedOrthoType = true;
730  }
731  }
732  params_->set("Orthogonalization", orthoType_);
733 
734  // Check which orthogonalization constant to use.
735  if (params->isParameter ("Orthogonalization Constant")) {
736  if (params->isType<MagnitudeType> ("Orthogonalization Constant")) {
737  orthoKappa_ = params->get ("Orthogonalization Constant",
738  static_cast<MagnitudeType> (DefaultSolverParameters::orthoKappa));
739  }
740  else {
741  orthoKappa_ = params->get ("Orthogonalization Constant",
743  }
744 
745  // Update parameter in our list.
746  params_->set ("Orthogonalization Constant", orthoKappa_);
747  if (orthoType_ == "DGKS") {
748  if (orthoKappa_ > 0 && ! ortho_.is_null() && !changedOrthoType) {
749  typedef DGKSOrthoManager<ScalarType, MV, OP> ortho_type;
750  rcp_dynamic_cast<ortho_type> (ortho_)->setDepTol (orthoKappa_);
751  }
752  }
753  }
754 
755  // Create orthogonalization manager if we need to.
756  if (ortho_.is_null() || changedOrthoType) {
758  Teuchos::RCP<Teuchos::ParameterList> paramsOrtho; // can be null
759  if (orthoType_=="DGKS" && orthoKappa_ > 0) {
760  paramsOrtho = Teuchos::rcp(new Teuchos::ParameterList());
761  paramsOrtho->set ("depTol", orthoKappa_ );
762  }
763 
764  ortho_ = factory.makeMatOrthoManager (orthoType_, Teuchos::null, printer_, label_, paramsOrtho);
765  }
766 
767  // Convergence
768 
769  // Check for convergence tolerance
770  if (params->isParameter ("Convergence Tolerance")) {
771  if (params->isType<MagnitudeType> ("Convergence Tolerance")) {
772  convtol_ = params->get ("Convergence Tolerance",
773  static_cast<MagnitudeType> (DefaultSolverParameters::convTol));
774  }
775  else {
776  convtol_ = params->get ("Convergence Tolerance", DefaultSolverParameters::convTol);
777  }
778 
779  // Update parameter in our list and residual tests.
780  params_->set ("Convergence Tolerance", convtol_);
781  if (! impConvTest_.is_null ()) {
782  impConvTest_->setTolerance (convtol_);
783  }
784  if (! expConvTest_.is_null ()) {
785  expConvTest_->setTolerance (convtol_);
786  }
787  }
788 
789  // Grab the user defined residual scaling
790  bool userDefinedResidualScalingUpdated = false;
791  if (params->isParameter ("User Defined Residual Scaling")) {
792  MagnitudeType tempResScaleFactor = DefaultSolverParameters::resScaleFactor;
793  if (params->isType<MagnitudeType> ("User Defined Residual Scaling")) {
794  tempResScaleFactor = params->get ("User Defined Residual Scaling",
795  static_cast<MagnitudeType> (DefaultSolverParameters::resScaleFactor));
796  }
797  else {
798  tempResScaleFactor = params->get ("User Defined Residual Scaling",
800  }
801 
802  // Only update the scaling if it's different.
803  if (resScaleFactor_ != tempResScaleFactor) {
804  resScaleFactor_ = tempResScaleFactor;
805  userDefinedResidualScalingUpdated = true;
806  }
807 
808  if(userDefinedResidualScalingUpdated)
809  {
810  if (! params->isParameter ("Implicit Residual Scaling") && ! impConvTest_.is_null ()) {
811  try {
812  if(impResScale_ == "User Provided")
813  impConvTest_->defineScaleForm (Belos::UserProvided, Belos::TwoNorm, resScaleFactor_);
814  }
815  catch (std::exception& e) {
816  // Make sure the convergence test gets constructed again.
817  isSTSet_ = false;
818  }
819  }
820  if (! params->isParameter ("Explicit Residual Scaling") && ! expConvTest_.is_null ()) {
821  try {
822  if(expResScale_ == "User Provided")
823  expConvTest_->defineScaleForm (Belos::UserProvided, Belos::TwoNorm, resScaleFactor_);
824  }
825  catch (std::exception& e) {
826  // Make sure the convergence test gets constructed again.
827  isSTSet_ = false;
828  }
829  }
830  }
831  }
832 
833  // Check for a change in scaling, if so we need to build new residual tests.
834  if (params->isParameter ("Implicit Residual Scaling")) {
835  const std::string tempImpResScale =
836  Teuchos::getParameter<std::string> (*params, "Implicit Residual Scaling");
837 
838  // Only update the scaling if it's different.
839  if (impResScale_ != tempImpResScale) {
840  Belos::ScaleType impResScaleType = convertStringToScaleType (tempImpResScale);
841  impResScale_ = tempImpResScale;
842 
843  // Update parameter in our list and residual tests
844  params_->set ("Implicit Residual Scaling", impResScale_);
845  if (! impConvTest_.is_null ()) {
846  try {
847  if(impResScale_ == "User Provided")
848  impConvTest_->defineScaleForm (impResScaleType, Belos::TwoNorm, resScaleFactor_);
849  else
850  impConvTest_->defineScaleForm (impResScaleType, Belos::TwoNorm);
851  }
852  catch (std::exception& e) {
853  // Make sure the convergence test gets constructed again.
854  isSTSet_ = false;
855  }
856  }
857  }
858  else if (userDefinedResidualScalingUpdated) {
859  Belos::ScaleType impResScaleType = convertStringToScaleType (impResScale_);
860 
861  if (! impConvTest_.is_null ()) {
862  try {
863  if(impResScale_ == "User Provided")
864  impConvTest_->defineScaleForm (impResScaleType, Belos::TwoNorm, resScaleFactor_);
865  }
866  catch (std::exception& e) {
867  // Make sure the convergence test gets constructed again.
868  isSTSet_ = false;
869  }
870  }
871  }
872  }
873 
874  if (params->isParameter ("Explicit Residual Scaling")) {
875  const std::string tempExpResScale =
876  Teuchos::getParameter<std::string> (*params, "Explicit Residual Scaling");
877 
878  // Only update the scaling if it's different.
879  if (expResScale_ != tempExpResScale) {
880  Belos::ScaleType expResScaleType = convertStringToScaleType (tempExpResScale);
881  expResScale_ = tempExpResScale;
882 
883  // Update parameter in our list and residual tests
884  params_->set ("Explicit Residual Scaling", expResScale_);
885  if (! expConvTest_.is_null ()) {
886  try {
887  if(expResScale_ == "User Provided")
888  expConvTest_->defineScaleForm (expResScaleType, Belos::TwoNorm, resScaleFactor_);
889  else
890  expConvTest_->defineScaleForm (expResScaleType, Belos::TwoNorm);
891  }
892  catch (std::exception& e) {
893  // Make sure the convergence test gets constructed again.
894  isSTSet_ = false;
895  }
896  }
897  }
898  else if (userDefinedResidualScalingUpdated) {
899  Belos::ScaleType expResScaleType = convertStringToScaleType (expResScale_);
900 
901  if (! expConvTest_.is_null ()) {
902  try {
903  if(expResScale_ == "User Provided")
904  expConvTest_->defineScaleForm (expResScaleType, Belos::TwoNorm, resScaleFactor_);
905  }
906  catch (std::exception& e) {
907  // Make sure the convergence test gets constructed again.
908  isSTSet_ = false;
909  }
910  }
911  }
912  }
913 
914 
915  if (params->isParameter ("Show Maximum Residual Norm Only")) {
916  showMaxResNormOnly_ =
917  Teuchos::getParameter<bool> (*params, "Show Maximum Residual Norm Only");
918 
919  // Update parameter in our list and residual tests
920  params_->set ("Show Maximum Residual Norm Only", showMaxResNormOnly_);
921  if (! impConvTest_.is_null ()) {
922  impConvTest_->setShowMaxResNormOnly (showMaxResNormOnly_);
923  }
924  if (! expConvTest_.is_null ()) {
925  expConvTest_->setShowMaxResNormOnly (showMaxResNormOnly_);
926  }
927  }
928 
929  // Create status tests if we need to.
930 
931  // Get the deflation quorum, or number of converged systems before deflation is allowed
932  if (params->isParameter("Deflation Quorum")) {
933  defQuorum_ = params->get("Deflation Quorum", defQuorum_);
935  defQuorum_ > blockSize_, std::invalid_argument,
936  "Belos::PseudoBlockGmresSolMgr::setParameters: "
937  "The \"Deflation Quorum\" parameter (= " << defQuorum_ << ") must not be "
938  "larger than \"Block Size\" (= " << blockSize_ << ").");
939  params_->set ("Deflation Quorum", defQuorum_);
940  if (! impConvTest_.is_null ()) {
941  impConvTest_->setQuorum (defQuorum_);
942  }
943  if (! expConvTest_.is_null ()) {
944  expConvTest_->setQuorum (defQuorum_);
945  }
946  }
947 
948  // Create the timer if we need to.
949  if (timerSolve_ == Teuchos::null) {
950  std::string solveLabel = label_ + ": PseudoBlockGmresSolMgr total solve time";
951 #ifdef BELOS_TEUCHOS_TIME_MONITOR
952  timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
953 #endif
954  }
955 
956  // Inform the solver manager that the current parameters were set.
957  isSet_ = true;
958 }
959 
960 
961 template<class ScalarType, class MV, class OP>
962 void
964  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &userConvStatusTest,
965  const typename StatusTestCombo<ScalarType,MV,OP>::ComboType &comboType
966  )
967 {
968  userConvStatusTest_ = userConvStatusTest;
969  comboType_ = comboType;
970 }
971 
972 template<class ScalarType, class MV, class OP>
973 void
975  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &debugStatusTest
976  )
977 {
978  debugStatusTest_ = debugStatusTest;
979 }
980 
981 
982 
983 template<class ScalarType, class MV, class OP>
986 {
988  if (is_null(validPL)) {
989  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
990  // Set all the valid parameters and their default values.
991 
992  // The static_cast is to resolve an issue with older clang versions which
993  // would cause the constexpr to link fail. With c++17 the problem is resolved.
995  pl->set("Convergence Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::convTol),
996  "The relative residual tolerance that needs to be achieved by the\n"
997  "iterative solver in order for the linear system to be declared converged.");
998  pl->set("Maximum Restarts", static_cast<int>(maxRestarts_default_),
999  "The maximum number of restarts allowed for each\n"
1000  "set of RHS solved.");
1001  pl->set("Maximum Iterations", static_cast<int>(maxIters_default_),
1002  "The maximum number of block iterations allowed for each\n"
1003  "set of RHS solved.");
1004  pl->set("Num Blocks", static_cast<int>(numBlocks_default_),
1005  "The maximum number of vectors allowed in the Krylov subspace\n"
1006  "for each set of RHS solved.");
1007  pl->set("Block Size", static_cast<int>(blockSize_default_),
1008  "The number of RHS solved simultaneously.");
1009  pl->set("Verbosity", static_cast<int>(verbosity_default_),
1010  "What type(s) of solver information should be outputted\n"
1011  "to the output stream.");
1012  pl->set("Output Style", static_cast<int>(outputStyle_default_),
1013  "What style is used for the solver information outputted\n"
1014  "to the output stream.");
1015  pl->set("Output Frequency", static_cast<int>(outputFreq_default_),
1016  "How often convergence information should be outputted\n"
1017  "to the output stream.");
1018  pl->set("Deflation Quorum", static_cast<int>(defQuorum_default_),
1019  "The number of linear systems that need to converge before\n"
1020  "they are deflated. This number should be <= block size.");
1021  pl->set("Output Stream", Teuchos::rcpFromRef(std::cout),
1022  "A reference-counted pointer to the output stream where all\n"
1023  "solver output is sent.");
1024  pl->set("Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
1025  "When convergence information is printed, only show the maximum\n"
1026  "relative residual norm when the block size is greater than one.");
1027  pl->set("Implicit Residual Scaling", static_cast<const char *>(impResScale_default_),
1028  "The type of scaling used in the implicit residual convergence test.");
1029  pl->set("Explicit Residual Scaling", static_cast<const char *>(expResScale_default_),
1030  "The type of scaling used in the explicit residual convergence test.");
1031  pl->set("Timer Label", static_cast<const char *>(label_default_),
1032  "The string to use as a prefix for the timer labels.");
1033  pl->set("Orthogonalization", static_cast<const char *>(orthoType_default_),
1034  "The type of orthogonalization to use.");
1035  pl->set("Orthogonalization Constant",static_cast<MagnitudeType>(DefaultSolverParameters::orthoKappa),
1036  "The constant used by DGKS orthogonalization to determine\n"
1037  "whether another step of classical Gram-Schmidt is necessary.");
1038  pl->sublist("User Status Tests");
1039  pl->set("User Status Tests Combo Type", "SEQ",
1040  "Type of logical combination operation of user-defined\n"
1041  "and/or solver-specific status tests.");
1042  validPL = pl;
1043  }
1044  return validPL;
1045 }
1046 
1047 // Check the status test versus the defined linear problem
1048 template<class ScalarType, class MV, class OP>
1050 
1051  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
1052  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestGenResNorm_t;
1053  typedef Belos::StatusTestImpResNorm<ScalarType,MV,OP> StatusTestImpResNorm_t;
1054 
1055  // Basic test checks maximum iterations and native residual.
1056  maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
1057 
1058  // If there is a left preconditioner, we create a combined status test that checks the implicit
1059  // and then explicit residual norm to see if we have convergence.
1060  if ( !Teuchos::is_null(problem_->getLeftPrec()) ) {
1061  expResTest_ = true;
1062  }
1063 
1064  if (expResTest_) {
1065 
1066  // Implicit residual test, using the native residual to determine if convergence was achieved.
1067  Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
1068  Teuchos::rcp( new StatusTestGenResNorm_t( convtol_, defQuorum_ ) );
1069  if(impResScale_ == "User Provided")
1070  tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm, resScaleFactor_ );
1071  else
1072  tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
1073 
1074  tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
1075  impConvTest_ = tmpImpConvTest;
1076 
1077  // Explicit residual test once the native residual is below the tolerance
1078  Teuchos::RCP<StatusTestGenResNorm_t> tmpExpConvTest =
1079  Teuchos::rcp( new StatusTestGenResNorm_t( convtol_, defQuorum_ ) );
1080  tmpExpConvTest->defineResForm( StatusTestGenResNorm_t::Explicit, Belos::TwoNorm );
1081  if(expResScale_ == "User Provided")
1082  tmpExpConvTest->defineScaleForm( convertStringToScaleType(expResScale_), Belos::TwoNorm, resScaleFactor_ );
1083  else
1084  tmpExpConvTest->defineScaleForm( convertStringToScaleType(expResScale_), Belos::TwoNorm );
1085  tmpExpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
1086  expConvTest_ = tmpExpConvTest;
1087 
1088  // The convergence test is a combination of the "cheap" implicit test and explicit test.
1089  convTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
1090  }
1091  else {
1092 
1093  // Implicit residual test, using the native residual to determine if convergence was achieved.
1094  // Use test that checks for loss of accuracy.
1095  Teuchos::RCP<StatusTestImpResNorm_t> tmpImpConvTest =
1096  Teuchos::rcp( new StatusTestImpResNorm_t( convtol_, defQuorum_ ) );
1097  if(impResScale_ == "User Provided")
1098  tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm, resScaleFactor_ );
1099  else
1100  tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
1101  tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
1102  impConvTest_ = tmpImpConvTest;
1103 
1104  // Set the explicit and total convergence test to this implicit test that checks for accuracy loss.
1105  expConvTest_ = impConvTest_;
1106  convTest_ = impConvTest_;
1107  }
1108 
1109  if (nonnull(userConvStatusTest_) ) {
1110  // Check if this is a combination of several StatusTestResNorm objects
1111  Teuchos::RCP<StatusTestCombo_t> tmpComboTest = Teuchos::rcp_dynamic_cast<StatusTestCombo_t>(userConvStatusTest_);
1112  if (tmpComboTest != Teuchos::null) {
1113  std::vector<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > tmpVec = tmpComboTest->getStatusTests();
1114  comboType_ = tmpComboTest->getComboType();
1115  const int numResTests = static_cast<int>(tmpVec.size());
1116  auto newConvTest = Teuchos::rcp(new StatusTestCombo_t(comboType_));
1117  newConvTest->addStatusTest(convTest_);
1118  for (int j = 0; j < numResTests; ++j) {
1119  newConvTest->addStatusTest(tmpVec[j]);
1120  }
1121  convTest_ = newConvTest;
1122  }
1123  else{
1124  // Override the overall convergence test with the users convergence test
1125  convTest_ = Teuchos::rcp(
1126  new StatusTestCombo_t( comboType_, convTest_, userConvStatusTest_ ) );
1127  // brief output style not compatible with more general combinations
1128  //outputStyle_ = Belos::General;
1129  // NOTE: Above, you have to run the other convergence tests also because
1130  // the logic in this class depends on it. This is very unfortunate.
1131  }
1132  }
1133 
1134  sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
1135 
1136  // Add debug status test, if one is provided by the user
1137  if (nonnull(debugStatusTest_) ) {
1138  // Add debug convergence test
1139  Teuchos::rcp_dynamic_cast<StatusTestCombo_t>(sTest_)->addStatusTest( debugStatusTest_ );
1140  }
1141 
1142  // Create the status test output class.
1143  // This class manages and formats the output from the status test.
1144  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_, taggedTests_ );
1145  outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
1146 
1147  // Set the solver string for the output test
1148  std::string solverDesc = " Pseudo Block Gmres ";
1149  outputTest_->setSolverDesc( solverDesc );
1150 
1151 
1152  // The status test is now set.
1153  isSTSet_ = true;
1154 
1155  return false;
1156 }
1157 
1158 
1159 // solve()
1160 template<class ScalarType, class MV, class OP>
1162 
1163  // Set the current parameters if they were not set before.
1164  // NOTE: This may occur if the user generated the solver manager with the default constructor and
1165  // then didn't set any parameters using setParameters().
1166  if (!isSet_) { setParameters( params_ ); }
1167 
1169  "Belos::PseudoBlockGmresSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
1170 
1171  // Check if we have to create the status tests.
1172  if (!isSTSet_ || (!expResTest_ && !Teuchos::is_null(problem_->getLeftPrec())) ) {
1174  "Belos::BlockGmresSolMgr::solve(): Linear problem and requested status tests are incompatible.");
1175  }
1176 
1177  // Create indices for the linear systems to be solved.
1178  int startPtr = 0;
1179  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1180  int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1181 
1182  std::vector<int> currIdx( numCurrRHS );
1183  blockSize_ = numCurrRHS;
1184  for (int i=0; i<numCurrRHS; ++i)
1185  { currIdx[i] = startPtr+i; }
1186 
1187  // Inform the linear problem of the current linear system to solve.
1188  problem_->setLSIndex( currIdx );
1189 
1191  // Parameter list
1192  Teuchos::ParameterList plist;
1193  plist.set("Num Blocks",numBlocks_);
1194 
1195  // Reset the status test.
1196  outputTest_->reset();
1197  loaDetected_ = false;
1198 
1199  // Assume convergence is achieved, then let any failed convergence set this to false.
1200  bool isConverged = true;
1201 
1203  // BlockGmres solver
1204 
1206  = Teuchos::rcp( new PseudoBlockGmresIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );
1207 
1208  // Enter solve() iterations
1209  {
1210 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1211  Teuchos::TimeMonitor slvtimer(*timerSolve_);
1212 #endif
1213 
1214  while ( numRHS2Solve > 0 ) {
1215 
1216  // Reset the active / converged vectors from this block
1217  std::vector<int> convRHSIdx;
1218  std::vector<int> currRHSIdx( currIdx );
1219  currRHSIdx.resize(numCurrRHS);
1220 
1221  // Set the current number of blocks with the pseudo Gmres iteration.
1222  block_gmres_iter->setNumBlocks( numBlocks_ );
1223 
1224  // Reset the number of iterations.
1225  block_gmres_iter->resetNumIters();
1226 
1227  // Reset the number of calls that the status test output knows about.
1228  outputTest_->resetNumCalls();
1229 
1230  // Get a new state struct and initialize the solver.
1232 
1233  // Create the first block in the current Krylov basis for each right-hand side.
1234  std::vector<int> index(1);
1235  Teuchos::RCP<MV> tmpV, R_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
1236  newState.V.resize( blockSize_ );
1237  newState.Z.resize( blockSize_ );
1238  for (int i=0; i<blockSize_; ++i) {
1239  index[0]=i;
1240  tmpV = MVT::CloneViewNonConst( *R_0, index );
1241 
1242  // Get a matrix to hold the orthonormalization coefficients.
1245 
1246  // Orthonormalize the new V_0
1247  int rank = ortho_->normalize( *tmpV, tmpZ );
1249  "Belos::PseudoBlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors.");
1250 
1251  newState.V[i] = tmpV;
1252  newState.Z[i] = tmpZ;
1253  }
1254 
1255  newState.curDim = 0;
1256  block_gmres_iter->initialize(newState);
1257  int numRestarts = 0;
1258 
1259  while(1) {
1260 
1261  // tell block_gmres_iter to iterate
1262  try {
1263  block_gmres_iter->iterate();
1264 
1266  //
1267  // check convergence first
1268  //
1270  if ( convTest_->getStatus() == Passed ) {
1271 
1272  if ( expConvTest_->getLOADetected() ) {
1273  // Oops! There was a loss of accuracy (LOA) for one or
1274  // more right-hand sides. That means the implicit
1275  // (a.k.a. "native") residuals claim convergence,
1276  // whereas the explicit (hence expConvTest_, i.e.,
1277  // "explicit convergence test") residuals have not
1278  // converged.
1279  //
1280  // We choose to deal with this situation by deflating
1281  // out the affected right-hand sides and moving on.
1282  loaDetected_ = true;
1283  printer_->stream(Warnings) <<
1284  "Belos::PseudoBlockGmresSolMgr::solve(): Warning! Solver has experienced a loss of accuracy!" << std::endl;
1285  isConverged = false;
1286  }
1287 
1288  // Figure out which linear systems converged.
1289  std::vector<int> convIdx = expConvTest_->convIndices();
1290 
1291  // If the number of converged linear systems is equal to the
1292  // number of current linear systems, then we are done with this block.
1293  if (convIdx.size() == currRHSIdx.size())
1294  break; // break from while(1){block_gmres_iter->iterate()}
1295 
1296  // Get a new state struct and initialize the solver.
1298 
1299  // Inform the linear problem that we are finished with this current linear system.
1300  problem_->setCurrLS();
1301 
1302  // Get the state.
1303  PseudoBlockGmresIterState<ScalarType,MV> oldState = block_gmres_iter->getState();
1304  int curDim = oldState.curDim;
1305 
1306  // Get a new state struct and reset currRHSIdx to have the right-hand sides that
1307  // are left to converge for this block.
1308  int have = 0;
1309  std::vector<int> oldRHSIdx( currRHSIdx );
1310  std::vector<int> defRHSIdx;
1311  for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
1312  bool found = false;
1313  for (unsigned int j=0; j<convIdx.size(); ++j) {
1314  if (currRHSIdx[i] == convIdx[j]) {
1315  found = true;
1316  break;
1317  }
1318  }
1319  if (found) {
1320  defRHSIdx.push_back( i );
1321  }
1322  else {
1323  defState.V.push_back( Teuchos::rcp_const_cast<MV>( oldState.V[i] ) );
1324  defState.Z.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,ScalarType> >( oldState.Z[i] ) );
1325  defState.H.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseMatrix<int,ScalarType> >( oldState.H[i] ) );
1326  defState.sn.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,ScalarType> >( oldState.sn[i] ) );
1327  defState.cs.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,MagnitudeType> >(oldState.cs[i] ) );
1328  currRHSIdx[have] = currRHSIdx[i];
1329  have++;
1330  }
1331  }
1332  defRHSIdx.resize(currRHSIdx.size()-have);
1333  currRHSIdx.resize(have);
1334 
1335  // Compute the current solution that needs to be deflated if this solver has taken any steps.
1336  if (curDim) {
1337  Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1338  Teuchos::RCP<MV> defUpdate = MVT::CloneViewNonConst( *update, defRHSIdx );
1339 
1340  // Set the deflated indices so we can update the solution.
1341  problem_->setLSIndex( convIdx );
1342 
1343  // Update the linear problem.
1344  problem_->updateSolution( defUpdate, true );
1345  }
1346 
1347  // Set the remaining indices after deflation.
1348  problem_->setLSIndex( currRHSIdx );
1349 
1350  // Set the dimension of the subspace, which is the same as the old subspace size.
1351  defState.curDim = curDim;
1352 
1353  // Initialize the solver with the deflated system.
1354  block_gmres_iter->initialize(defState);
1355  }
1357  //
1358  // check for maximum iterations
1359  //
1361  else if ( maxIterTest_->getStatus() == Passed ) {
1362  // we don't have convergence
1363  isConverged = false;
1364  break; // break from while(1){block_gmres_iter->iterate()}
1365  }
1367  //
1368  // check for restarting, i.e. the subspace is full
1369  //
1371  else if ( block_gmres_iter->getCurSubspaceDim() == block_gmres_iter->getMaxSubspaceDim() ) {
1372 
1373  if ( numRestarts >= maxRestarts_ ) {
1374  isConverged = false;
1375  break; // break from while(1){block_gmres_iter->iterate()}
1376  }
1377  numRestarts++;
1378 
1379  printer_->stream(Debug) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
1380 
1381  // Update the linear problem.
1382  Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1383  problem_->updateSolution( update, true );
1384 
1385  // Get the state.
1386  PseudoBlockGmresIterState<ScalarType,MV> oldState = block_gmres_iter->getState();
1387 
1388  // Set the new state.
1390  newstate.V.resize(currRHSIdx.size());
1391  newstate.Z.resize(currRHSIdx.size());
1392 
1393  // Compute the restart vectors
1394  // NOTE: Force the linear problem to update the current residual since the solution was updated.
1395  R_0 = MVT::Clone( *(problem_->getInitPrecResVec()), currRHSIdx.size() );
1396  problem_->computeCurrPrecResVec( &*R_0 );
1397  for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
1398  index[0] = i; // index(1) vector declared on line 891
1399 
1400  tmpV = MVT::CloneViewNonConst( *R_0, index );
1401 
1402  // Get a matrix to hold the orthonormalization coefficients.
1405 
1406  // Orthonormalize the new V_0
1407  int rank = ortho_->normalize( *tmpV, tmpZ );
1409  "Belos::PseudoBlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors after the restart.");
1410 
1411  newstate.V[i] = tmpV;
1412  newstate.Z[i] = tmpZ;
1413  }
1414 
1415  // Initialize the solver.
1416  newstate.curDim = 0;
1417  block_gmres_iter->initialize(newstate);
1418 
1419  } // end of restarting
1420 
1422  //
1423  // we returned from iterate(), but none of our status tests Passed.
1424  // something is wrong, and it is probably our fault.
1425  //
1427 
1428  else {
1429  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
1430  "Belos::PseudoBlockGmresSolMgr::solve(): Invalid return from PseudoBlockGmresIter::iterate().");
1431  }
1432  }
1433  catch (const PseudoBlockGmresIterOrthoFailure &e) {
1434 
1435  // Try to recover the most recent least-squares solution
1436  block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
1437 
1438  // Check to see if the most recent least-squares solution yielded convergence.
1439  sTest_->checkStatus( &*block_gmres_iter );
1440  if (convTest_->getStatus() != Passed)
1441  isConverged = false;
1442  break;
1443  }
1444  catch (const StatusTestNaNError& e) {
1445  // A NaN was detected in the solver. Set the solution to zero and return unconverged.
1446  achievedTol_ = MT::one();
1447  Teuchos::RCP<MV> X = problem_->getLHS();
1448  MVT::MvInit( *X, SCT::zero() );
1449  printer_->stream(Warnings) << "Belos::PseudoBlockGmresSolMgr::solve(): Warning! NaN has been detected!"
1450  << std::endl;
1451  return Unconverged;
1452  }
1453  catch (const std::exception &e) {
1454  printer_->stream(Errors) << "Error! Caught std::exception in PseudoBlockGmresIter::iterate() at iteration "
1455  << block_gmres_iter->getNumIters() << std::endl
1456  << e.what() << std::endl;
1457  throw;
1458  }
1459  }
1460 
1461  // Compute the current solution.
1462  // Update the linear problem.
1463  if (nonnull(userConvStatusTest_)) {
1464  //std::cout << "\nnonnull(userConvStatusTest_)\n";
1465  Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1466  problem_->updateSolution( update, true );
1467  }
1468  else if (nonnull(expConvTest_->getSolution())) {
1469  //std::cout << "\nexpConvTest_->getSolution()\n";
1470  Teuchos::RCP<MV> newX = expConvTest_->getSolution();
1471  Teuchos::RCP<MV> curX = problem_->getCurrLHSVec();
1472  MVT::MvAddMv( 0.0, *newX, 1.0, *newX, *curX );
1473  }
1474  else {
1475  //std::cout << "\nblock_gmres_iter->getCurrentUpdate()\n";
1476  Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1477  problem_->updateSolution( update, true );
1478  }
1479 
1480  // Inform the linear problem that we are finished with this block linear system.
1481  problem_->setCurrLS();
1482 
1483  // Update indices for the linear systems to be solved.
1484  startPtr += numCurrRHS;
1485  numRHS2Solve -= numCurrRHS;
1486  if ( numRHS2Solve > 0 ) {
1487  numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1488 
1489  blockSize_ = numCurrRHS;
1490  currIdx.resize( numCurrRHS );
1491  for (int i=0; i<numCurrRHS; ++i)
1492  { currIdx[i] = startPtr+i; }
1493 
1494  // Adapt the status test quorum if we need to.
1495  if (defQuorum_ > blockSize_) {
1496  if (impConvTest_ != Teuchos::null)
1497  impConvTest_->setQuorum( blockSize_ );
1498  if (expConvTest_ != Teuchos::null)
1499  expConvTest_->setQuorum( blockSize_ );
1500  }
1501 
1502  // Set the next indices.
1503  problem_->setLSIndex( currIdx );
1504  }
1505  else {
1506  currIdx.resize( numRHS2Solve );
1507  }
1508 
1509  }// while ( numRHS2Solve > 0 )
1510 
1511  }
1512 
1513  // print final summary
1514  sTest_->print( printer_->stream(FinalSummary) );
1515 
1516  // print timing information
1517 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1518  // Calling summarize() can be expensive, so don't call unless the
1519  // user wants to print out timing details. summarize() will do all
1520  // the work even if it's passed a "black hole" output stream.
1521  if (verbosity_ & TimingDetails)
1522  Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
1523 #endif
1524 
1525  // get iteration information for this solve
1526  numIters_ = maxIterTest_->getNumIters();
1527 
1528  // Save the convergence test value ("achieved tolerance") for this
1529  // solve. For this solver, convTest_ may either be a single
1530  // residual norm test, or a combination of two residual norm tests.
1531  // In the latter case, the master convergence test convTest_ is a
1532  // SEQ combo of the implicit resp. explicit tests. If the implicit
1533  // test never passes, then the explicit test won't ever be executed.
1534  // This manifests as expConvTest_->getTestValue()->size() < 1. We
1535  // deal with this case by using the values returned by
1536  // impConvTest_->getTestValue().
1537  {
1538  // We'll fetch the vector of residual norms one way or the other.
1539  const std::vector<MagnitudeType>* pTestValues = NULL;
1540  if (expResTest_) {
1541  pTestValues = expConvTest_->getTestValue();
1542  if (pTestValues == NULL || pTestValues->size() < 1) {
1543  pTestValues = impConvTest_->getTestValue();
1544  }
1545  }
1546  else {
1547  // Only the implicit residual norm test is being used.
1548  pTestValues = impConvTest_->getTestValue();
1549  }
1550  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1551  "Belos::PseudoBlockGmresSolMgr::solve(): The implicit convergence test's "
1552  "getTestValue() method returned NULL. Please report this bug to the "
1553  "Belos developers.");
1554  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1555  "Belos::PseudoBlockGmresSolMgr::solve(): The implicit convergence test's "
1556  "getTestValue() method returned a vector of length zero. Please report "
1557  "this bug to the Belos developers.");
1558 
1559  // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
1560  // achieved tolerances for all vectors in the current solve(), or
1561  // just for the vectors from the last deflation?
1562  achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1563  }
1564 
1565  if (!isConverged || loaDetected_) {
1566  return Unconverged; // return from PseudoBlockGmresSolMgr::solve()
1567  }
1568  return Converged; // return from PseudoBlockGmresSolMgr::solve()
1569 }
1570 
1571 
1572 template<class ScalarType, class MV, class OP>
1574 {
1575  std::ostringstream out;
1576 
1577  out << "\"Belos::PseudoBlockGmresSolMgr\": {";
1578  if (this->getObjectLabel () != "") {
1579  out << "Label: " << this->getObjectLabel () << ", ";
1580  }
1581  out << "Num Blocks: " << numBlocks_
1582  << ", Maximum Iterations: " << maxIters_
1583  << ", Maximum Restarts: " << maxRestarts_
1584  << ", Convergence Tolerance: " << convtol_
1585  << "}";
1586  return out.str ();
1587 }
1588 
1589 
1590 template<class ScalarType, class MV, class OP>
1591 void
1594  const Teuchos::EVerbosityLevel verbLevel) const
1595 {
1597  using Teuchos::VERB_DEFAULT;
1598  using Teuchos::VERB_NONE;
1599  using Teuchos::VERB_LOW;
1600  // using Teuchos::VERB_MEDIUM;
1601  // using Teuchos::VERB_HIGH;
1602  // using Teuchos::VERB_EXTREME;
1603  using std::endl;
1604 
1605  // Set default verbosity if applicable.
1606  const Teuchos::EVerbosityLevel vl =
1607  (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
1608 
1609  if (vl != VERB_NONE) {
1610  Teuchos::OSTab tab0 (out);
1611 
1612  out << "\"Belos::PseudoBlockGmresSolMgr\":" << endl;
1613  Teuchos::OSTab tab1 (out);
1614  out << "Template parameters:" << endl;
1615  {
1616  Teuchos::OSTab tab2 (out);
1617  out << "ScalarType: " << TypeNameTraits<ScalarType>::name () << endl
1618  << "MV: " << TypeNameTraits<MV>::name () << endl
1619  << "OP: " << TypeNameTraits<OP>::name () << endl;
1620  }
1621  if (this->getObjectLabel () != "") {
1622  out << "Label: " << this->getObjectLabel () << endl;
1623  }
1624  out << "Num Blocks: " << numBlocks_ << endl
1625  << "Maximum Iterations: " << maxIters_ << endl
1626  << "Maximum Restarts: " << maxRestarts_ << endl
1627  << "Convergence Tolerance: " << convtol_ << endl;
1628  }
1629 }
1630 
1631 } // end Belos namespace
1632 
1633 #endif /* BELOS_PSEUDO_BLOCK_GMRES_SOLMGR_HPP */
static const double orthoKappa
DGKS orthogonalization constant.
Definition: BelosTypes.hpp:267
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Definition: BelosTypes.cpp:74
std::vector< Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > > H
The current Hessenberg matrix.
Collection of types and exceptions used within the Belos solvers.
ComboType
The test can be either the AND of all the component tests, or the OR of all the component tests...
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > &params) override
Set the parameters the solver manager should use to solve the linear problem.
Class which manages the output and verbosity of the Belos solvers.
bool is_null(const boost::shared_ptr< T > &p)
ScaleType
The type of scaling to use on the residual norm value.
Definition: BelosTypes.hpp:88
Convergence test using the implicit residual norm(s), with an explicit residual norm(s) check for los...
T & get(const std::string &name, T def_value)
PseudoBlockGmresSolMgrOrthoFailure(const std::string &what_arg)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
static RCP< Time > getNewCounter(const std::string &name)
bool is_null(const std::shared_ptr< T > &p)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void setDebugStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &debugStatusTest) override
Set a debug status test that will be checked at the same time as the top-level status test...
An abstract class of StatusTest for stopping criteria using residual norms.
std::vector< Teuchos::RCP< const MV > > V
The current Krylov basis.
An implementation of StatusTestResNorm using a family of residual norms.
Ordinal numParams() const
static const double convTol
Default convergence tolerance.
Definition: BelosTypes.hpp:261
A pure virtual class for defining the status tests for the Belos iterative solvers.
Structure to contain pointers to PseudoBlockGmresIter state variables.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
A factory class for generating StatusTestOutput objects.
Traits class which defines basic operations on multivectors.
bool isParameter(const std::string &name) const
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
Definition: BelosTypes.hpp:174
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Pure virtual base class which describes the basic interface for a solver manager. ...
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const override
Print the object with the given verbosity level to a FancyOStream.
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
bool isSublist(const std::string &name) const
PseudoBlockGmresSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i...
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, MagnitudeType > > > cs
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
Interface to standard and &quot;pseudoblock&quot; GMRES.
int getNumIters() const override
Iteration count for the most recent call to solve().
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:123
The Belos::SolverManager is a templated virtual base class that defines the basic interface that any ...
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
A list of valid default parameters for this solver.
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
static const EVerbosityLevel verbLevel_default
bool isLOADetected() const override
Whether a &quot;loss of accuracy&quot; was detected during the last solve().
bool nonnull(const boost::shared_ptr< T > &p)
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, ScalarType > > > Z
The current right-hand side of the least squares system RY = Z.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem to solve.
Belos concrete class for performing the pseudo-block GMRES iteration.
bool isType(const std::string &name) const
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
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.
std::string description() const override
Return a one-line description of this object.
const StatusTestResNorm< ScalarType, MV, OP > * getResidualStatusTest() const
Return the residual status test.
Class which defines basic traits for the operator type.
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
int curDim
The current dimension of the reduction.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:28
ReturnType solve() override
This method performs possibly repeated calls to the underlying linear solver&#39;s iterate() routine unti...
Default parameters common to most Belos solvers.
Definition: BelosTypes.hpp:251
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, ScalarType > > > sn
The current Given&#39;s rotation coefficients.
static const double resScaleFactor
User-defined residual scaling factor.
Definition: BelosTypes.hpp:270
Belos header file which uses auto-configuration information to include necessary C++ headers...
This class implements the pseudo-block GMRES iteration, where a block Krylov subspace is constructed ...
Factory to build a set of status tests from a parameter list.
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
The current parameters for this solver.
PseudoBlockGmresIterOrthoFailure is thrown when the orthogonalization manager is unable to generate o...
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Return a reference to the linear problem being solved by this solver manager.
PseudoBlockGmresSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate...
virtual void setUserConvStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &userConvStatusTest, const typename StatusTestCombo< ScalarType, MV, OP >::ComboType &comboType=StatusTestCombo< ScalarType, MV, OP >::SEQ) override
Set a custom status test.
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.
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 > &params)
Return an instance of the specified MatOrthoManager subclass.

Generated on Mon Jul 15 2024 09:24:24 for Belos by doxygen 1.8.5