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

Generated for Belos by doxygen 1.8.14