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 //
4 // Belos: Block Linear Solvers Package
5 // Copyright 2004 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
42 #ifndef BELOS_PSEUDO_BLOCK_GMRES_SOLMGR_HPP
43 #define BELOS_PSEUDO_BLOCK_GMRES_SOLMGR_HPP
44 
49 #include "BelosConfigDefs.hpp"
50 #include "BelosTypes.hpp"
51 
52 #include "BelosLinearProblem.hpp"
53 #include "BelosSolverManager.hpp"
54 
59 #ifdef HAVE_BELOS_TSQR
60 # include "BelosTsqrOrthoManager.hpp"
61 #endif // HAVE_BELOS_TSQR
64 #include "BelosOutputManager.hpp"
65 #include "Teuchos_BLAS.hpp"
66 #ifdef BELOS_TEUCHOS_TIME_MONITOR
67 #include "Teuchos_TimeMonitor.hpp"
68 #endif
69 
77 namespace Belos {
78 
80 
81 
89  PseudoBlockGmresSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
90  {}};
91 
99  PseudoBlockGmresSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
100  {}};
101 
125  template<class ScalarType, class MV, class OP>
126  class PseudoBlockGmresSolMgr : public SolverManager<ScalarType,MV,OP> {
127 
128  private:
132  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
134 
135  public:
136 
138 
139 
148 
261 
264 
268  }
270 
272 
273 
274  const LinearProblem<ScalarType,MV,OP>& getProblem() const override {
275  return *problem_;
276  }
277 
280 
283 
290  return Teuchos::tuple(timerSolve_);
291  }
292 
303  MagnitudeType achievedTol() const override {
304  return achievedTol_;
305  }
306 
308  int getNumIters() const override {
309  return numIters_;
310  }
311 
367  bool isLOADetected() const override { return loaDetected_; }
368 
370 
372 
373 
375  void setProblem (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem) override {
376  problem_ = problem;
377  }
378 
380  void setParameters (const Teuchos::RCP<Teuchos::ParameterList> &params) override;
381 
388  virtual void setUserConvStatusTest(
389  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &userConvStatusTest,
390  const typename StatusTestCombo<ScalarType,MV,OP>::ComboType &comboType =
392  ) override;
393 
395  void setDebugStatusTest( const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &debugStatusTest ) override;
396 
398 
400 
401 
405  void reset( const ResetType type ) override { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
407 
409 
410 
428  ReturnType solve() override;
429 
431 
434 
441  void
443  const Teuchos::EVerbosityLevel verbLevel =
445 
447  std::string description () const override;
448 
450 
451  private:
452 
467  bool checkStatusTest();
468 
471 
472  // Output manager.
474  Teuchos::RCP<std::ostream> outputStream_;
475 
476  // Status tests.
477  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > userConvStatusTest_;
482  Teuchos::RCP<StatusTestResNorm<ScalarType,MV,OP> > impConvTest_, expConvTest_;
486 
487  // Orthogonalization manager.
489 
490  // Current parameter list.
492 
493  // Default solver values.
494  static constexpr int maxRestarts_default_ = 20;
495  static constexpr int maxIters_default_ = 1000;
496  static constexpr bool showMaxResNormOnly_default_ = false;
497  static constexpr int blockSize_default_ = 1;
498  static constexpr int numBlocks_default_ = 300;
499  static constexpr int verbosity_default_ = Belos::Errors;
500  static constexpr int outputStyle_default_ = Belos::General;
501  static constexpr int outputFreq_default_ = -1;
502  static constexpr int defQuorum_default_ = 1;
503  static constexpr const char * impResScale_default_ = "Norm of Preconditioned Initial Residual";
504  static constexpr const char * expResScale_default_ = "Norm of Initial Residual";
505  static constexpr const char * label_default_ = "Belos";
506  static constexpr const char * orthoType_default_ = "ICGS";
507  static constexpr std::ostream * outputStream_default_ = &std::cout;
508 
509  // Current solver values.
510  MagnitudeType convtol_, orthoKappa_, achievedTol_;
511  int maxRestarts_, maxIters_, numIters_;
512  int blockSize_, numBlocks_, verbosity_, outputStyle_, outputFreq_, defQuorum_;
513  bool showMaxResNormOnly_;
514  std::string orthoType_;
515  std::string impResScale_, expResScale_;
516  MagnitudeType resScaleFactor_;
517 
518  // Timers.
519  std::string label_;
520  Teuchos::RCP<Teuchos::Time> timerSolve_;
521 
522  // Internal state variables.
523  bool isSet_, isSTSet_, expResTest_;
524  bool loaDetected_;
525  };
526 
527 
528 // Empty Constructor
529 template<class ScalarType, class MV, class OP>
531  outputStream_(Teuchos::rcp(outputStream_default_,false)),
532  taggedTests_(Teuchos::null),
533  convtol_(DefaultSolverParameters::convTol),
534  orthoKappa_(DefaultSolverParameters::orthoKappa),
535  achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
536  maxRestarts_(maxRestarts_default_),
537  maxIters_(maxIters_default_),
538  numIters_(0),
539  blockSize_(blockSize_default_),
540  numBlocks_(numBlocks_default_),
541  verbosity_(verbosity_default_),
542  outputStyle_(outputStyle_default_),
543  outputFreq_(outputFreq_default_),
544  defQuorum_(defQuorum_default_),
545  showMaxResNormOnly_(showMaxResNormOnly_default_),
546  orthoType_(orthoType_default_),
547  impResScale_(impResScale_default_),
548  expResScale_(expResScale_default_),
549  resScaleFactor_(DefaultSolverParameters::resScaleFactor),
550  label_(label_default_),
551  isSet_(false),
552  isSTSet_(false),
553  expResTest_(false),
554  loaDetected_(false)
555 {}
556 
557 // Basic Constructor
558 template<class ScalarType, class MV, class OP>
562  problem_(problem),
563  outputStream_(Teuchos::rcp(outputStream_default_,false)),
564  taggedTests_(Teuchos::null),
565  convtol_(DefaultSolverParameters::convTol),
566  orthoKappa_(DefaultSolverParameters::orthoKappa),
567  achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
568  maxRestarts_(maxRestarts_default_),
569  maxIters_(maxIters_default_),
570  numIters_(0),
571  blockSize_(blockSize_default_),
572  numBlocks_(numBlocks_default_),
573  verbosity_(verbosity_default_),
574  outputStyle_(outputStyle_default_),
575  outputFreq_(outputFreq_default_),
576  defQuorum_(defQuorum_default_),
577  showMaxResNormOnly_(showMaxResNormOnly_default_),
578  orthoType_(orthoType_default_),
579  impResScale_(impResScale_default_),
580  expResScale_(expResScale_default_),
581  resScaleFactor_(DefaultSolverParameters::resScaleFactor),
582  label_(label_default_),
583  isSet_(false),
584  isSTSet_(false),
585  expResTest_(false),
586  loaDetected_(false)
587 {
588  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
589 
590  // If the parameter list pointer is null, then set the current parameters to the default parameter list.
591  if (!is_null(pl)) {
592  // Set the parameters using the list that was passed in.
593  setParameters( pl );
594  }
595 }
596 
597 template<class ScalarType, class MV, class OP>
598 void
601 {
603  using Teuchos::parameterList;
604  using Teuchos::rcp;
605  using Teuchos::rcp_dynamic_cast;
606 
607  // Create the internal parameter list if one doesn't already exist.
608  if (params_ == Teuchos::null) {
609  params_ = parameterList (*getValidParameters ());
610  } else {
611  // TAW: 3/8/2016: do not validate sub parameter lists as they
612  // might not have a pre-defined structure
613  // e.g. user-specified status tests
614  // The Belos Pseudo Block GMRES parameters on the first level are
615  // not affected and verified.
616  params->validateParameters (*getValidParameters (), 0);
617  }
618 
619  // Check for maximum number of restarts
620  if (params->isParameter ("Maximum Restarts")) {
621  maxRestarts_ = params->get ("Maximum Restarts", maxRestarts_default_);
622 
623  // Update parameter in our list.
624  params_->set ("Maximum Restarts", maxRestarts_);
625  }
626 
627  // Check for maximum number of iterations
628  if (params->isParameter ("Maximum Iterations")) {
629  maxIters_ = params->get ("Maximum Iterations", maxIters_default_);
630 
631  // Update parameter in our list and in status test.
632  params_->set ("Maximum Iterations", maxIters_);
633  if (! maxIterTest_.is_null ()) {
634  maxIterTest_->setMaxIters (maxIters_);
635  }
636  }
637 
638  // Check for blocksize
639  if (params->isParameter ("Block Size")) {
640  blockSize_ = params->get ("Block Size", blockSize_default_);
642  blockSize_ <= 0, std::invalid_argument,
643  "Belos::PseudoBlockGmresSolMgr::setParameters: "
644  "The \"Block Size\" parameter must be strictly positive, "
645  "but you specified a value of " << blockSize_ << ".");
646 
647  // Update parameter in our list.
648  params_->set ("Block Size", blockSize_);
649  }
650 
651  // Check for the maximum number of blocks.
652  if (params->isParameter ("Num Blocks")) {
653  numBlocks_ = params->get ("Num Blocks", numBlocks_default_);
655  numBlocks_ <= 0, std::invalid_argument,
656  "Belos::PseudoBlockGmresSolMgr::setParameters: "
657  "The \"Num Blocks\" parameter must be strictly positive, "
658  "but you specified a value of " << numBlocks_ << ".");
659 
660  // Update parameter in our list.
661  params_->set ("Num Blocks", numBlocks_);
662  }
663 
664  // Check to see if the timer label changed.
665  if (params->isParameter ("Timer Label")) {
666  const std::string tempLabel = params->get ("Timer Label", label_default_);
667 
668  // Update parameter in our list and solver timer
669  if (tempLabel != label_) {
670  label_ = tempLabel;
671  params_->set ("Timer Label", label_);
672  const std::string solveLabel =
673  label_ + ": PseudoBlockGmresSolMgr total solve time";
674 #ifdef BELOS_TEUCHOS_TIME_MONITOR
675  timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
676 #endif // BELOS_TEUCHOS_TIME_MONITOR
677  if (ortho_ != Teuchos::null) {
678  ortho_->setLabel( label_ );
679  }
680  }
681  }
682 
683  // Check if the orthogonalization changed.
684  if (params->isParameter ("Orthogonalization")) {
685  std::string tempOrthoType = params->get ("Orthogonalization", orthoType_default_);
686 #ifdef HAVE_BELOS_TSQR
688  tempOrthoType != "DGKS" && tempOrthoType != "ICGS" &&
689  tempOrthoType != "IMGS" && tempOrthoType != "TSQR",
690  std::invalid_argument,
691  "Belos::PseudoBlockGmresSolMgr::setParameters: "
692  "The \"Orthogonalization\" parameter must be one of \"DGKS\", \"ICGS\", "
693  "\"IMGS\", or \"TSQR\".");
694 #else
696  tempOrthoType != "DGKS" && tempOrthoType != "ICGS" &&
697  tempOrthoType != "IMGS",
698  std::invalid_argument,
699  "Belos::PseudoBlockGmresSolMgr::setParameters: "
700  "The \"Orthogonalization\" parameter must be one of \"DGKS\", \"ICGS\", "
701  "or \"IMGS\".");
702 #endif // HAVE_BELOS_TSQR
703 
704  if (tempOrthoType != orthoType_) {
705  orthoType_ = tempOrthoType;
706  params_->set("Orthogonalization", orthoType_);
707  // Create orthogonalization manager
708  if (orthoType_ == "DGKS") {
709  typedef DGKSOrthoManager<ScalarType, MV, OP> ortho_type;
710  if (orthoKappa_ <= 0) {
711  ortho_ = rcp (new ortho_type (label_));
712  }
713  else {
714  ortho_ = rcp (new ortho_type (label_));
715  rcp_dynamic_cast<ortho_type> (ortho_)->setDepTol (orthoKappa_);
716  }
717  }
718  else if (orthoType_ == "ICGS") {
719  typedef ICGSOrthoManager<ScalarType, MV, OP> ortho_type;
720  ortho_ = rcp (new ortho_type (label_));
721  }
722  else if (orthoType_ == "IMGS") {
723  typedef IMGSOrthoManager<ScalarType, MV, OP> ortho_type;
724  ortho_ = rcp (new ortho_type (label_));
725  }
726 #ifdef HAVE_BELOS_TSQR
727  else if (orthoType_ == "TSQR") {
728  typedef TsqrMatOrthoManager<ScalarType, MV, OP> ortho_type;
729  ortho_ = rcp (new ortho_type (label_));
730  }
731 #endif // HAVE_BELOS_TSQR
732  }
733  }
734 
735  // Check which orthogonalization constant to use.
736  if (params->isParameter ("Orthogonalization Constant")) {
737  if (params->isType<MagnitudeType> ("Orthogonalization Constant")) {
738  orthoKappa_ = params->get ("Orthogonalization Constant",
739  static_cast<MagnitudeType> (DefaultSolverParameters::orthoKappa));
740  }
741  else {
742  orthoKappa_ = params->get ("Orthogonalization Constant",
744  }
745 
746  // Update parameter in our list.
747  params_->set ("Orthogonalization Constant", orthoKappa_);
748  if (orthoType_ == "DGKS") {
749  if (orthoKappa_ > 0 && ! ortho_.is_null ()) {
750  typedef DGKSOrthoManager<ScalarType, MV, OP> ortho_type;
751  rcp_dynamic_cast<ortho_type> (ortho_)->setDepTol (orthoKappa_);
752  }
753  }
754  }
755 
756  // Check for a change in verbosity level
757  if (params->isParameter ("Verbosity")) {
758  if (Teuchos::isParameterType<int> (*params, "Verbosity")) {
759  verbosity_ = params->get ("Verbosity", verbosity_default_);
760  } else {
761  verbosity_ = (int) Teuchos::getParameter<Belos::MsgType> (*params, "Verbosity");
762  }
763 
764  // Update parameter in our list.
765  params_->set ("Verbosity", verbosity_);
766  if (! printer_.is_null ()) {
767  printer_->setVerbosity (verbosity_);
768  }
769  }
770 
771  // Check for a change in output style.
772  if (params->isParameter ("Output Style")) {
773  if (Teuchos::isParameterType<int> (*params, "Output Style")) {
774  outputStyle_ = params->get ("Output Style", outputStyle_default_);
775  } else {
776  outputStyle_ = (int) Teuchos::getParameter<Belos::OutputType> (*params, "Output Style");
777  }
778 
779  // Update parameter in our list.
780  params_->set ("Output Style", verbosity_);
781  if (! outputTest_.is_null ()) {
782  isSTSet_ = false;
783  }
784 
785  }
786 
787  // Check if user has specified his own status tests
788  if (params->isSublist ("User Status Tests")) {
789  Teuchos::ParameterList userStatusTestsList = params->sublist("User Status Tests",true);
790  if ( userStatusTestsList.numParams() > 0 ) {
791  std::string userCombo_string = params->get<std::string>("User Status Tests Combo Type", "SEQ");
793  setUserConvStatusTest( testFactory->buildStatusTests(userStatusTestsList), testFactory->stringToComboType(userCombo_string) );
794  taggedTests_ = testFactory->getTaggedTests();
795  isSTSet_ = false;
796  }
797  }
798 
799  // output stream
800  if (params->isParameter ("Output Stream")) {
801  outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> > (*params, "Output Stream");
802 
803  // Update parameter in our list.
804  params_->set("Output Stream", outputStream_);
805  if (! printer_.is_null ()) {
806  printer_->setOStream (outputStream_);
807  }
808  }
809 
810  // frequency level
811  if (verbosity_ & Belos::StatusTestDetails) {
812  if (params->isParameter ("Output Frequency")) {
813  outputFreq_ = params->get ("Output Frequency", outputFreq_default_);
814  }
815 
816  // Update parameter in out list and output status test.
817  params_->set ("Output Frequency", outputFreq_);
818  if (! outputTest_.is_null ()) {
819  outputTest_->setOutputFrequency (outputFreq_);
820  }
821  }
822 
823  // Create output manager if we need to.
824  if (printer_.is_null ()) {
825  printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
826  }
827 
828  // Convergence
829 
830  // Check for convergence tolerance
831  if (params->isParameter ("Convergence Tolerance")) {
832  if (params->isType<MagnitudeType> ("Convergence Tolerance")) {
833  convtol_ = params->get ("Convergence Tolerance",
834  static_cast<MagnitudeType> (DefaultSolverParameters::convTol));
835  }
836  else {
837  convtol_ = params->get ("Convergence Tolerance", DefaultSolverParameters::convTol);
838  }
839 
840  // Update parameter in our list and residual tests.
841  params_->set ("Convergence Tolerance", convtol_);
842  if (! impConvTest_.is_null ()) {
843  impConvTest_->setTolerance (convtol_);
844  }
845  if (! expConvTest_.is_null ()) {
846  expConvTest_->setTolerance (convtol_);
847  }
848  }
849 
850  // Grab the user defined residual scaling
851  bool userDefinedResidualScalingUpdated = false;
852  if (params->isParameter ("User Defined Residual Scaling")) {
853  MagnitudeType tempResScaleFactor = DefaultSolverParameters::resScaleFactor;
854  if (params->isType<MagnitudeType> ("User Defined Residual Scaling")) {
855  tempResScaleFactor = params->get ("User Defined Residual Scaling",
856  static_cast<MagnitudeType> (DefaultSolverParameters::resScaleFactor));
857  }
858  else {
859  tempResScaleFactor = params->get ("User Defined Residual Scaling",
861  }
862 
863  // Only update the scaling if it's different.
864  if (resScaleFactor_ != tempResScaleFactor) {
865  resScaleFactor_ = tempResScaleFactor;
866  userDefinedResidualScalingUpdated = true;
867  }
868 
869  if(userDefinedResidualScalingUpdated)
870  {
871  if (! params->isParameter ("Implicit Residual Scaling") && ! impConvTest_.is_null ()) {
872  try {
873  if(impResScale_ == "User Provided")
874  impConvTest_->defineScaleForm (Belos::UserProvided, Belos::TwoNorm, resScaleFactor_);
875  }
876  catch (std::exception& e) {
877  // Make sure the convergence test gets constructed again.
878  isSTSet_ = false;
879  }
880  }
881  if (! params->isParameter ("Explicit Residual Scaling") && ! expConvTest_.is_null ()) {
882  try {
883  if(expResScale_ == "User Provided")
884  expConvTest_->defineScaleForm (Belos::UserProvided, Belos::TwoNorm, resScaleFactor_);
885  }
886  catch (std::exception& e) {
887  // Make sure the convergence test gets constructed again.
888  isSTSet_ = false;
889  }
890  }
891  }
892  }
893 
894  // Check for a change in scaling, if so we need to build new residual tests.
895  if (params->isParameter ("Implicit Residual Scaling")) {
896  const std::string tempImpResScale =
897  Teuchos::getParameter<std::string> (*params, "Implicit Residual Scaling");
898 
899  // Only update the scaling if it's different.
900  if (impResScale_ != tempImpResScale) {
901  Belos::ScaleType impResScaleType = convertStringToScaleType (tempImpResScale);
902  impResScale_ = tempImpResScale;
903 
904  // Update parameter in our list and residual tests
905  params_->set ("Implicit Residual Scaling", impResScale_);
906  if (! impConvTest_.is_null ()) {
907  try {
908  if(impResScale_ == "User Provided")
909  impConvTest_->defineScaleForm (impResScaleType, Belos::TwoNorm, resScaleFactor_);
910  else
911  impConvTest_->defineScaleForm (impResScaleType, Belos::TwoNorm);
912  }
913  catch (std::exception& e) {
914  // Make sure the convergence test gets constructed again.
915  isSTSet_ = false;
916  }
917  }
918  }
919  else if (userDefinedResidualScalingUpdated) {
920  Belos::ScaleType impResScaleType = convertStringToScaleType (impResScale_);
921 
922  if (! impConvTest_.is_null ()) {
923  try {
924  if(impResScale_ == "User Provided")
925  impConvTest_->defineScaleForm (impResScaleType, Belos::TwoNorm, resScaleFactor_);
926  }
927  catch (std::exception& e) {
928  // Make sure the convergence test gets constructed again.
929  isSTSet_ = false;
930  }
931  }
932  }
933  }
934 
935  if (params->isParameter ("Explicit Residual Scaling")) {
936  const std::string tempExpResScale =
937  Teuchos::getParameter<std::string> (*params, "Explicit Residual Scaling");
938 
939  // Only update the scaling if it's different.
940  if (expResScale_ != tempExpResScale) {
941  Belos::ScaleType expResScaleType = convertStringToScaleType (tempExpResScale);
942  expResScale_ = tempExpResScale;
943 
944  // Update parameter in our list and residual tests
945  params_->set ("Explicit Residual Scaling", expResScale_);
946  if (! expConvTest_.is_null ()) {
947  try {
948  if(expResScale_ == "User Provided")
949  expConvTest_->defineScaleForm (expResScaleType, Belos::TwoNorm, resScaleFactor_);
950  else
951  expConvTest_->defineScaleForm (expResScaleType, Belos::TwoNorm);
952  }
953  catch (std::exception& e) {
954  // Make sure the convergence test gets constructed again.
955  isSTSet_ = false;
956  }
957  }
958  }
959  else if (userDefinedResidualScalingUpdated) {
960  Belos::ScaleType expResScaleType = convertStringToScaleType (expResScale_);
961 
962  if (! expConvTest_.is_null ()) {
963  try {
964  if(expResScale_ == "User Provided")
965  expConvTest_->defineScaleForm (expResScaleType, Belos::TwoNorm, resScaleFactor_);
966  }
967  catch (std::exception& e) {
968  // Make sure the convergence test gets constructed again.
969  isSTSet_ = false;
970  }
971  }
972  }
973  }
974 
975 
976  if (params->isParameter ("Show Maximum Residual Norm Only")) {
977  showMaxResNormOnly_ =
978  Teuchos::getParameter<bool> (*params, "Show Maximum Residual Norm Only");
979 
980  // Update parameter in our list and residual tests
981  params_->set ("Show Maximum Residual Norm Only", showMaxResNormOnly_);
982  if (! impConvTest_.is_null ()) {
983  impConvTest_->setShowMaxResNormOnly (showMaxResNormOnly_);
984  }
985  if (! expConvTest_.is_null ()) {
986  expConvTest_->setShowMaxResNormOnly (showMaxResNormOnly_);
987  }
988  }
989 
990  // Create status tests if we need to.
991 
992  // Get the deflation quorum, or number of converged systems before deflation is allowed
993  if (params->isParameter("Deflation Quorum")) {
994  defQuorum_ = params->get("Deflation Quorum", defQuorum_);
996  defQuorum_ > blockSize_, std::invalid_argument,
997  "Belos::PseudoBlockGmresSolMgr::setParameters: "
998  "The \"Deflation Quorum\" parameter (= " << defQuorum_ << ") must not be "
999  "larger than \"Block Size\" (= " << blockSize_ << ").");
1000  params_->set ("Deflation Quorum", defQuorum_);
1001  if (! impConvTest_.is_null ()) {
1002  impConvTest_->setQuorum (defQuorum_);
1003  }
1004  if (! expConvTest_.is_null ()) {
1005  expConvTest_->setQuorum (defQuorum_);
1006  }
1007  }
1008 
1009  // Create orthogonalization manager if we need to.
1010  if (ortho_.is_null ()) {
1011  params_->set("Orthogonalization", orthoType_);
1012  if (orthoType_ == "DGKS") {
1013  typedef DGKSOrthoManager<ScalarType, MV, OP> ortho_type;
1014  if (orthoKappa_ <= 0) {
1015  ortho_ = rcp (new ortho_type (label_));
1016  }
1017  else {
1018  ortho_ = rcp (new ortho_type (label_));
1019  rcp_dynamic_cast<ortho_type> (ortho_)->setDepTol (orthoKappa_);
1020  }
1021  }
1022  else if (orthoType_ == "ICGS") {
1023  typedef ICGSOrthoManager<ScalarType, MV, OP> ortho_type;
1024  ortho_ = rcp (new ortho_type (label_));
1025  }
1026  else if (orthoType_ == "IMGS") {
1027  typedef IMGSOrthoManager<ScalarType, MV, OP> ortho_type;
1028  ortho_ = rcp (new ortho_type (label_));
1029  }
1030 #ifdef HAVE_BELOS_TSQR
1031  else if (orthoType_ == "TSQR") {
1032  typedef TsqrMatOrthoManager<ScalarType, MV, OP> ortho_type;
1033  ortho_ = rcp (new ortho_type (label_));
1034  }
1035 #endif // HAVE_BELOS_TSQR
1036  else {
1037 #ifdef HAVE_BELOS_TSQR
1039  orthoType_ != "ICGS" && orthoType_ != "DGKS" &&
1040  orthoType_ != "IMGS" && orthoType_ != "TSQR",
1041  std::logic_error,
1042  "Belos::PseudoBlockGmresSolMgr::setParameters(): "
1043  "Invalid orthogonalization type \"" << orthoType_ << "\".");
1044 #else
1046  orthoType_ != "ICGS" && orthoType_ != "DGKS" &&
1047  orthoType_ != "IMGS",
1048  std::logic_error,
1049  "Belos::PseudoBlockGmresSolMgr::setParameters(): "
1050  "Invalid orthogonalization type \"" << orthoType_ << "\".");
1051 #endif // HAVE_BELOS_TSQR
1052  }
1053  }
1054 
1055  // Create the timer if we need to.
1056  if (timerSolve_ == Teuchos::null) {
1057  std::string solveLabel = label_ + ": PseudoBlockGmresSolMgr total solve time";
1058 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1059  timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
1060 #endif
1061  }
1062 
1063  // Inform the solver manager that the current parameters were set.
1064  isSet_ = true;
1065 }
1066 
1067 
1068 template<class ScalarType, class MV, class OP>
1069 void
1071  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &userConvStatusTest,
1072  const typename StatusTestCombo<ScalarType,MV,OP>::ComboType &comboType
1073  )
1074 {
1075  userConvStatusTest_ = userConvStatusTest;
1076  comboType_ = comboType;
1077 }
1078 
1079 template<class ScalarType, class MV, class OP>
1080 void
1082  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &debugStatusTest
1083  )
1084 {
1085  debugStatusTest_ = debugStatusTest;
1086 }
1087 
1088 
1089 
1090 template<class ScalarType, class MV, class OP>
1093 {
1095  if (is_null(validPL)) {
1096  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
1097  // Set all the valid parameters and their default values.
1098 
1099  // The static_cast is to resolve an issue with older clang versions which
1100  // would cause the constexpr to link fail. With c++17 the problem is resolved.
1101  pl= Teuchos::rcp( new Teuchos::ParameterList() );
1102  pl->set("Convergence Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::convTol),
1103  "The relative residual tolerance that needs to be achieved by the\n"
1104  "iterative solver in order for the linear system to be declared converged.");
1105  pl->set("Maximum Restarts", static_cast<int>(maxRestarts_default_),
1106  "The maximum number of restarts allowed for each\n"
1107  "set of RHS solved.");
1108  pl->set("Maximum Iterations", static_cast<int>(maxIters_default_),
1109  "The maximum number of block iterations allowed for each\n"
1110  "set of RHS solved.");
1111  pl->set("Num Blocks", static_cast<int>(numBlocks_default_),
1112  "The maximum number of vectors allowed in the Krylov subspace\n"
1113  "for each set of RHS solved.");
1114  pl->set("Block Size", static_cast<int>(blockSize_default_),
1115  "The number of RHS solved simultaneously.");
1116  pl->set("Verbosity", static_cast<int>(verbosity_default_),
1117  "What type(s) of solver information should be outputted\n"
1118  "to the output stream.");
1119  pl->set("Output Style", static_cast<int>(outputStyle_default_),
1120  "What style is used for the solver information outputted\n"
1121  "to the output stream.");
1122  pl->set("Output Frequency", static_cast<int>(outputFreq_default_),
1123  "How often convergence information should be outputted\n"
1124  "to the output stream.");
1125  pl->set("Deflation Quorum", static_cast<int>(defQuorum_default_),
1126  "The number of linear systems that need to converge before\n"
1127  "they are deflated. This number should be <= block size.");
1128  pl->set("Output Stream", Teuchos::rcp(outputStream_default_,false),
1129  "A reference-counted pointer to the output stream where all\n"
1130  "solver output is sent.");
1131  pl->set("Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
1132  "When convergence information is printed, only show the maximum\n"
1133  "relative residual norm when the block size is greater than one.");
1134  pl->set("Implicit Residual Scaling", static_cast<const char *>(impResScale_default_),
1135  "The type of scaling used in the implicit residual convergence test.");
1136  pl->set("Explicit Residual Scaling", static_cast<const char *>(expResScale_default_),
1137  "The type of scaling used in the explicit residual convergence test.");
1138  pl->set("Timer Label", static_cast<const char *>(label_default_),
1139  "The string to use as a prefix for the timer labels.");
1140  pl->set("Orthogonalization", static_cast<const char *>(orthoType_default_),
1141  "The type of orthogonalization to use.");
1142  pl->set("Orthogonalization Constant",static_cast<MagnitudeType>(DefaultSolverParameters::orthoKappa),
1143  "The constant used by DGKS orthogonalization to determine\n"
1144  "whether another step of classical Gram-Schmidt is necessary.");
1145  pl->sublist("User Status Tests");
1146  pl->set("User Status Tests Combo Type", "SEQ",
1147  "Type of logical combination operation of user-defined\n"
1148  "and/or solver-specific status tests.");
1149  validPL = pl;
1150  }
1151  return validPL;
1152 }
1153 
1154 // Check the status test versus the defined linear problem
1155 template<class ScalarType, class MV, class OP>
1157 
1158  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
1159  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestGenResNorm_t;
1160  typedef Belos::StatusTestImpResNorm<ScalarType,MV,OP> StatusTestImpResNorm_t;
1161 
1162  // Basic test checks maximum iterations and native residual.
1163  maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
1164 
1165  // If there is a left preconditioner, we create a combined status test that checks the implicit
1166  // and then explicit residual norm to see if we have convergence.
1167  if ( !Teuchos::is_null(problem_->getLeftPrec()) ) {
1168  expResTest_ = true;
1169  }
1170 
1171  if (expResTest_) {
1172 
1173  // Implicit residual test, using the native residual to determine if convergence was achieved.
1174  Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
1175  Teuchos::rcp( new StatusTestGenResNorm_t( convtol_, defQuorum_ ) );
1176  if(impResScale_ == "User Provided")
1177  tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm, resScaleFactor_ );
1178  else
1179  tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
1180 
1181  tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
1182  impConvTest_ = tmpImpConvTest;
1183 
1184  // Explicit residual test once the native residual is below the tolerance
1185  Teuchos::RCP<StatusTestGenResNorm_t> tmpExpConvTest =
1186  Teuchos::rcp( new StatusTestGenResNorm_t( convtol_, defQuorum_ ) );
1187  tmpExpConvTest->defineResForm( StatusTestGenResNorm_t::Explicit, Belos::TwoNorm );
1188  if(expResScale_ == "User Provided")
1189  tmpExpConvTest->defineScaleForm( convertStringToScaleType(expResScale_), Belos::TwoNorm, resScaleFactor_ );
1190  else
1191  tmpExpConvTest->defineScaleForm( convertStringToScaleType(expResScale_), Belos::TwoNorm );
1192  tmpExpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
1193  expConvTest_ = tmpExpConvTest;
1194 
1195  // The convergence test is a combination of the "cheap" implicit test and explicit test.
1196  convTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
1197  }
1198  else {
1199 
1200  // Implicit residual test, using the native residual to determine if convergence was achieved.
1201  // Use test that checks for loss of accuracy.
1202  Teuchos::RCP<StatusTestImpResNorm_t> tmpImpConvTest =
1203  Teuchos::rcp( new StatusTestImpResNorm_t( convtol_, defQuorum_ ) );
1204  if(impResScale_ == "User Provided")
1205  tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm, resScaleFactor_ );
1206  else
1207  tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
1208  tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
1209  impConvTest_ = tmpImpConvTest;
1210 
1211  // Set the explicit and total convergence test to this implicit test that checks for accuracy loss.
1212  expConvTest_ = impConvTest_;
1213  convTest_ = impConvTest_;
1214  }
1215 
1216  if (nonnull(userConvStatusTest_) ) {
1217  // Override the overall convergence test with the users convergence test
1218  convTest_ = Teuchos::rcp(
1219  new StatusTestCombo_t( comboType_, convTest_, userConvStatusTest_ ) );
1220  // brief output style not compatible with more general combinations
1221  //outputStyle_ = Belos::General;
1222  // NOTE: Above, you have to run the other convergence tests also because
1223  // the logic in this class depends on it. This is very unfortunate.
1224  }
1225 
1226  sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
1227 
1228  // Add debug status test, if one is provided by the user
1229  if (nonnull(debugStatusTest_) ) {
1230  // Add debug convergence test
1231  Teuchos::rcp_dynamic_cast<StatusTestCombo_t>(sTest_)->addStatusTest( debugStatusTest_ );
1232  }
1233 
1234  // Create the status test output class.
1235  // This class manages and formats the output from the status test.
1236  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_, taggedTests_ );
1237  outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
1238 
1239  // Set the solver string for the output test
1240  std::string solverDesc = " Pseudo Block Gmres ";
1241  outputTest_->setSolverDesc( solverDesc );
1242 
1243 
1244  // The status test is now set.
1245  isSTSet_ = true;
1246 
1247  return false;
1248 }
1249 
1250 
1251 // solve()
1252 template<class ScalarType, class MV, class OP>
1254 
1255  // Set the current parameters if they were not set before.
1256  // NOTE: This may occur if the user generated the solver manager with the default constructor and
1257  // then didn't set any parameters using setParameters().
1258  if (!isSet_) { setParameters( params_ ); }
1259 
1261 
1263  "Belos::PseudoBlockGmresSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
1264 
1265  // Check if we have to create the status tests.
1266  if (!isSTSet_ || (!expResTest_ && !Teuchos::is_null(problem_->getLeftPrec())) ) {
1268  "Belos::BlockGmresSolMgr::solve(): Linear problem and requested status tests are incompatible.");
1269  }
1270 
1271  // Create indices for the linear systems to be solved.
1272  int startPtr = 0;
1273  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1274  int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1275 
1276  std::vector<int> currIdx( numCurrRHS );
1277  blockSize_ = numCurrRHS;
1278  for (int i=0; i<numCurrRHS; ++i)
1279  { currIdx[i] = startPtr+i; }
1280 
1281  // Inform the linear problem of the current linear system to solve.
1282  problem_->setLSIndex( currIdx );
1283 
1285  // Parameter list
1286  Teuchos::ParameterList plist;
1287  plist.set("Num Blocks",numBlocks_);
1288 
1289  // Reset the status test.
1290  outputTest_->reset();
1291  loaDetected_ = false;
1292 
1293  // Assume convergence is achieved, then let any failed convergence set this to false.
1294  bool isConverged = true;
1295 
1297  // BlockGmres solver
1298 
1300  = Teuchos::rcp( new PseudoBlockGmresIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );
1301 
1302  // Enter solve() iterations
1303  {
1304 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1305  Teuchos::TimeMonitor slvtimer(*timerSolve_);
1306 #endif
1307 
1308  while ( numRHS2Solve > 0 ) {
1309 
1310  // Reset the active / converged vectors from this block
1311  std::vector<int> convRHSIdx;
1312  std::vector<int> currRHSIdx( currIdx );
1313  currRHSIdx.resize(numCurrRHS);
1314 
1315  // Set the current number of blocks with the pseudo Gmres iteration.
1316  block_gmres_iter->setNumBlocks( numBlocks_ );
1317 
1318  // Reset the number of iterations.
1319  block_gmres_iter->resetNumIters();
1320 
1321  // Reset the number of calls that the status test output knows about.
1322  outputTest_->resetNumCalls();
1323 
1324  // Get a new state struct and initialize the solver.
1326 
1327  // Create the first block in the current Krylov basis for each right-hand side.
1328  std::vector<int> index(1);
1329  Teuchos::RCP<MV> tmpV, R_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
1330  newState.V.resize( blockSize_ );
1331  newState.Z.resize( blockSize_ );
1332  for (int i=0; i<blockSize_; ++i) {
1333  index[0]=i;
1334  tmpV = MVT::CloneViewNonConst( *R_0, index );
1335 
1336  // Get a matrix to hold the orthonormalization coefficients.
1339 
1340  // Orthonormalize the new V_0
1341  int rank = ortho_->normalize( *tmpV, tmpZ );
1343  "Belos::PseudoBlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors.");
1344 
1345  newState.V[i] = tmpV;
1346  newState.Z[i] = tmpZ;
1347  }
1348 
1349  newState.curDim = 0;
1350  block_gmres_iter->initialize(newState);
1351  int numRestarts = 0;
1352 
1353  while(1) {
1354 
1355  // tell block_gmres_iter to iterate
1356  try {
1357  block_gmres_iter->iterate();
1358 
1360  //
1361  // check convergence first
1362  //
1364  if ( convTest_->getStatus() == Passed ) {
1365 
1366  if ( expConvTest_->getLOADetected() ) {
1367  // Oops! There was a loss of accuracy (LOA) for one or
1368  // more right-hand sides. That means the implicit
1369  // (a.k.a. "native") residuals claim convergence,
1370  // whereas the explicit (hence expConvTest_, i.e.,
1371  // "explicit convergence test") residuals have not
1372  // converged.
1373  //
1374  // We choose to deal with this situation by deflating
1375  // out the affected right-hand sides and moving on.
1376  loaDetected_ = true;
1377  printer_->stream(Warnings) <<
1378  "Belos::PseudoBlockGmresSolMgr::solve(): Warning! Solver has experienced a loss of accuracy!" << std::endl;
1379  isConverged = false;
1380  }
1381 
1382  // Figure out which linear systems converged.
1383  std::vector<int> convIdx = expConvTest_->convIndices();
1384 
1385  // If the number of converged linear systems is equal to the
1386  // number of current linear systems, then we are done with this block.
1387  if (convIdx.size() == currRHSIdx.size())
1388  break; // break from while(1){block_gmres_iter->iterate()}
1389 
1390  // Get a new state struct and initialize the solver.
1392 
1393  // Inform the linear problem that we are finished with this current linear system.
1394  problem_->setCurrLS();
1395 
1396  // Get the state.
1397  PseudoBlockGmresIterState<ScalarType,MV> oldState = block_gmres_iter->getState();
1398  int curDim = oldState.curDim;
1399 
1400  // Get a new state struct and reset currRHSIdx to have the right-hand sides that
1401  // are left to converge for this block.
1402  int have = 0;
1403  std::vector<int> oldRHSIdx( currRHSIdx );
1404  std::vector<int> defRHSIdx;
1405  for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
1406  bool found = false;
1407  for (unsigned int j=0; j<convIdx.size(); ++j) {
1408  if (currRHSIdx[i] == convIdx[j]) {
1409  found = true;
1410  break;
1411  }
1412  }
1413  if (found) {
1414  defRHSIdx.push_back( i );
1415  }
1416  else {
1417  defState.V.push_back( Teuchos::rcp_const_cast<MV>( oldState.V[i] ) );
1418  defState.Z.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,ScalarType> >( oldState.Z[i] ) );
1419  defState.H.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseMatrix<int,ScalarType> >( oldState.H[i] ) );
1420  defState.sn.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,ScalarType> >( oldState.sn[i] ) );
1421  defState.cs.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,MagnitudeType> >(oldState.cs[i] ) );
1422  currRHSIdx[have] = currRHSIdx[i];
1423  have++;
1424  }
1425  }
1426  defRHSIdx.resize(currRHSIdx.size()-have);
1427  currRHSIdx.resize(have);
1428 
1429  // Compute the current solution that needs to be deflated if this solver has taken any steps.
1430  if (curDim) {
1431  Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1432  Teuchos::RCP<MV> defUpdate = MVT::CloneViewNonConst( *update, defRHSIdx );
1433 
1434  // Set the deflated indices so we can update the solution.
1435  problem_->setLSIndex( convIdx );
1436 
1437  // Update the linear problem.
1438  problem_->updateSolution( defUpdate, true );
1439  }
1440 
1441  // Set the remaining indices after deflation.
1442  problem_->setLSIndex( currRHSIdx );
1443 
1444  // Set the dimension of the subspace, which is the same as the old subspace size.
1445  defState.curDim = curDim;
1446 
1447  // Initialize the solver with the deflated system.
1448  block_gmres_iter->initialize(defState);
1449  }
1451  //
1452  // check for maximum iterations
1453  //
1455  else if ( maxIterTest_->getStatus() == Passed ) {
1456  // we don't have convergence
1457  isConverged = false;
1458  break; // break from while(1){block_gmres_iter->iterate()}
1459  }
1461  //
1462  // check for restarting, i.e. the subspace is full
1463  //
1465  else if ( block_gmres_iter->getCurSubspaceDim() == block_gmres_iter->getMaxSubspaceDim() ) {
1466 
1467  if ( numRestarts >= maxRestarts_ ) {
1468  isConverged = false;
1469  break; // break from while(1){block_gmres_iter->iterate()}
1470  }
1471  numRestarts++;
1472 
1473  printer_->stream(Debug) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
1474 
1475  // Update the linear problem.
1476  Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1477  problem_->updateSolution( update, true );
1478 
1479  // Get the state.
1480  PseudoBlockGmresIterState<ScalarType,MV> oldState = block_gmres_iter->getState();
1481 
1482  // Set the new state.
1484  newstate.V.resize(currRHSIdx.size());
1485  newstate.Z.resize(currRHSIdx.size());
1486 
1487  // Compute the restart vectors
1488  // NOTE: Force the linear problem to update the current residual since the solution was updated.
1489  R_0 = MVT::Clone( *(problem_->getInitPrecResVec()), currRHSIdx.size() );
1490  problem_->computeCurrPrecResVec( &*R_0 );
1491  for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
1492  index[0] = i; // index(1) vector declared on line 891
1493 
1494  tmpV = MVT::CloneViewNonConst( *R_0, index );
1495 
1496  // Get a matrix to hold the orthonormalization coefficients.
1499 
1500  // Orthonormalize the new V_0
1501  int rank = ortho_->normalize( *tmpV, tmpZ );
1503  "Belos::PseudoBlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors after the restart.");
1504 
1505  newstate.V[i] = tmpV;
1506  newstate.Z[i] = tmpZ;
1507  }
1508 
1509  // Initialize the solver.
1510  newstate.curDim = 0;
1511  block_gmres_iter->initialize(newstate);
1512 
1513  } // end of restarting
1514 
1516  //
1517  // we returned from iterate(), but none of our status tests Passed.
1518  // something is wrong, and it is probably our fault.
1519  //
1521 
1522  else {
1523  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
1524  "Belos::PseudoBlockGmresSolMgr::solve(): Invalid return from PseudoBlockGmresIter::iterate().");
1525  }
1526  }
1527  catch (const PseudoBlockGmresIterOrthoFailure &e) {
1528 
1529  // Try to recover the most recent least-squares solution
1530  block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
1531 
1532  // Check to see if the most recent least-squares solution yielded convergence.
1533  sTest_->checkStatus( &*block_gmres_iter );
1534  if (convTest_->getStatus() != Passed)
1535  isConverged = false;
1536  break;
1537  }
1538  catch (const std::exception &e) {
1539  printer_->stream(Errors) << "Error! Caught std::exception in PseudoBlockGmresIter::iterate() at iteration "
1540  << block_gmres_iter->getNumIters() << std::endl
1541  << e.what() << std::endl;
1542  throw;
1543  }
1544  }
1545 
1546  // Compute the current solution.
1547  // Update the linear problem.
1548  if (nonnull(userConvStatusTest_)) {
1549  //std::cout << "\nnonnull(userConvStatusTest_)\n";
1550  Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1551  problem_->updateSolution( update, true );
1552  }
1553  else if (nonnull(expConvTest_->getSolution())) {
1554  //std::cout << "\nexpConvTest_->getSolution()\n";
1555  Teuchos::RCP<MV> newX = expConvTest_->getSolution();
1556  Teuchos::RCP<MV> curX = problem_->getCurrLHSVec();
1557  MVT::MvAddMv( 0.0, *newX, 1.0, *newX, *curX );
1558  }
1559  else {
1560  //std::cout << "\nblock_gmres_iter->getCurrentUpdate()\n";
1561  Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1562  problem_->updateSolution( update, true );
1563  }
1564 
1565  // Inform the linear problem that we are finished with this block linear system.
1566  problem_->setCurrLS();
1567 
1568  // Update indices for the linear systems to be solved.
1569  startPtr += numCurrRHS;
1570  numRHS2Solve -= numCurrRHS;
1571  if ( numRHS2Solve > 0 ) {
1572  numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1573 
1574  blockSize_ = numCurrRHS;
1575  currIdx.resize( numCurrRHS );
1576  for (int i=0; i<numCurrRHS; ++i)
1577  { currIdx[i] = startPtr+i; }
1578 
1579  // Adapt the status test quorum if we need to.
1580  if (defQuorum_ > blockSize_) {
1581  if (impConvTest_ != Teuchos::null)
1582  impConvTest_->setQuorum( blockSize_ );
1583  if (expConvTest_ != Teuchos::null)
1584  expConvTest_->setQuorum( blockSize_ );
1585  }
1586 
1587  // Set the next indices.
1588  problem_->setLSIndex( currIdx );
1589  }
1590  else {
1591  currIdx.resize( numRHS2Solve );
1592  }
1593 
1594  }// while ( numRHS2Solve > 0 )
1595 
1596  }
1597 
1598  // print final summary
1599  sTest_->print( printer_->stream(FinalSummary) );
1600 
1601  // print timing information
1602 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1603  // Calling summarize() can be expensive, so don't call unless the
1604  // user wants to print out timing details. summarize() will do all
1605  // the work even if it's passed a "black hole" output stream.
1606  if (verbosity_ & TimingDetails)
1607  Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
1608 #endif
1609 
1610  // get iteration information for this solve
1611  numIters_ = maxIterTest_->getNumIters();
1612 
1613  // Save the convergence test value ("achieved tolerance") for this
1614  // solve. For this solver, convTest_ may either be a single
1615  // residual norm test, or a combination of two residual norm tests.
1616  // In the latter case, the master convergence test convTest_ is a
1617  // SEQ combo of the implicit resp. explicit tests. If the implicit
1618  // test never passes, then the explicit test won't ever be executed.
1619  // This manifests as expConvTest_->getTestValue()->size() < 1. We
1620  // deal with this case by using the values returned by
1621  // impConvTest_->getTestValue().
1622  {
1623  // We'll fetch the vector of residual norms one way or the other.
1624  const std::vector<MagnitudeType>* pTestValues = NULL;
1625  if (expResTest_) {
1626  pTestValues = expConvTest_->getTestValue();
1627  if (pTestValues == NULL || pTestValues->size() < 1) {
1628  pTestValues = impConvTest_->getTestValue();
1629  }
1630  }
1631  else {
1632  // Only the implicit residual norm test is being used.
1633  pTestValues = impConvTest_->getTestValue();
1634  }
1635  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1636  "Belos::PseudoBlockGmresSolMgr::solve(): The implicit convergence test's "
1637  "getTestValue() method returned NULL. Please report this bug to the "
1638  "Belos developers.");
1639  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1640  "Belos::PseudoBlockGmresSolMgr::solve(): The implicit convergence test's "
1641  "getTestValue() method returned a vector of length zero. Please report "
1642  "this bug to the Belos developers.");
1643 
1644  // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
1645  // achieved tolerances for all vectors in the current solve(), or
1646  // just for the vectors from the last deflation?
1647  achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1648  }
1649 
1650  if (!isConverged || loaDetected_) {
1651  return Unconverged; // return from PseudoBlockGmresSolMgr::solve()
1652  }
1653  return Converged; // return from PseudoBlockGmresSolMgr::solve()
1654 }
1655 
1656 
1657 template<class ScalarType, class MV, class OP>
1659 {
1660  std::ostringstream out;
1661 
1662  out << "\"Belos::PseudoBlockGmresSolMgr\": {";
1663  if (this->getObjectLabel () != "") {
1664  out << "Label: " << this->getObjectLabel () << ", ";
1665  }
1666  out << "Num Blocks: " << numBlocks_
1667  << ", Maximum Iterations: " << maxIters_
1668  << ", Maximum Restarts: " << maxRestarts_
1669  << ", Convergence Tolerance: " << convtol_
1670  << "}";
1671  return out.str ();
1672 }
1673 
1674 
1675 template<class ScalarType, class MV, class OP>
1676 void
1679  const Teuchos::EVerbosityLevel verbLevel) const
1680 {
1682  using Teuchos::VERB_DEFAULT;
1683  using Teuchos::VERB_NONE;
1684  using Teuchos::VERB_LOW;
1685  // using Teuchos::VERB_MEDIUM;
1686  // using Teuchos::VERB_HIGH;
1687  // using Teuchos::VERB_EXTREME;
1688  using std::endl;
1689 
1690  // Set default verbosity if applicable.
1691  const Teuchos::EVerbosityLevel vl =
1692  (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
1693 
1694  if (vl != VERB_NONE) {
1695  Teuchos::OSTab tab0 (out);
1696 
1697  out << "\"Belos::PseudoBlockGmresSolMgr\":" << endl;
1698  Teuchos::OSTab tab1 (out);
1699  out << "Template parameters:" << endl;
1700  {
1701  Teuchos::OSTab tab2 (out);
1702  out << "ScalarType: " << TypeNameTraits<ScalarType>::name () << endl
1703  << "MV: " << TypeNameTraits<MV>::name () << endl
1704  << "OP: " << TypeNameTraits<OP>::name () << endl;
1705  }
1706  if (this->getObjectLabel () != "") {
1707  out << "Label: " << this->getObjectLabel () << endl;
1708  }
1709  out << "Num Blocks: " << numBlocks_ << endl
1710  << "Maximum Iterations: " << maxIters_ << endl
1711  << "Maximum Restarts: " << maxRestarts_ << endl
1712  << "Convergence Tolerance: " << convtol_ << endl;
1713  }
1714 }
1715 
1716 } // end Belos namespace
1717 
1718 #endif /* BELOS_PSEUDO_BLOCK_GMRES_SOLMGR_HPP */
static const double orthoKappa
DGKS orthogonalization constant.
Definition: BelosTypes.hpp:307
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Definition: BelosTypes.cpp:106
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:120
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...
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:301
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.
Iterated Modified Gram-Schmidt (IMGS) implementation of the Belos::OrthoManager class.
Traits class which defines basic operations on multivectors.
bool isParameter(const std::string &name) const
Classical Gram-Schmidt (with DGKS correction) implementation of the Belos::OrthoManager class...
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
Definition: BelosTypes.hpp:206
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:155
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.
Iterated Classical Gram-Schmidt (ICGS) implementation of the Belos::OrthoManager class.
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
static const EVerbosityLevel verbLevel_default
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
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.
Orthogonalization manager based on Tall Skinny QR (TSQR)
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:60
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:291
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:310
MatOrthoManager subclass using TSQR or DGKS.
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.

Generated on Fri Jun 5 2020 10:20:45 for Belos by doxygen 1.8.5