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