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 #ifdef HAVE_BELOS_TSQR
60 # include "BelosTsqrOrthoManager.hpp"
61 #endif // HAVE_BELOS_TSQR
64 #include "BelosOutputManager.hpp"
65 #include "Teuchos_BLAS.hpp"
66 #ifdef BELOS_TEUCHOS_TIME_MONITOR
67 #include "Teuchos_TimeMonitor.hpp"
68 #endif
69 
77 namespace Belos {
78 
80 
81 
89  PseudoBlockGmresSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
90  {}};
91 
99  PseudoBlockGmresSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
100  {}};
101 
125  template<class ScalarType, class MV, class OP>
126  class PseudoBlockGmresSolMgr : public SolverManager<ScalarType,MV,OP> {
127 
128  private:
132  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
134 
135  public:
136 
138 
139 
148 
260 
263 
267  }
269 
271 
272 
273  const LinearProblem<ScalarType,MV,OP>& getProblem() const override {
274  return *problem_;
275  }
276 
279 
282 
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
442  const Teuchos::EVerbosityLevel verbLevel =
444 
446  std::string description () const override;
447 
449 
450  private:
451 
466  bool checkStatusTest();
467 
470 
471  // Output manager.
473  Teuchos::RCP<std::ostream> outputStream_;
474 
475  // Status tests.
476  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > userConvStatusTest_;
481  Teuchos::RCP<StatusTestResNorm<ScalarType,MV,OP> > impConvTest_, expConvTest_;
485 
486  // Orthogonalization manager.
488 
489  // Current parameter list.
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.
509  MagnitudeType convtol_, orthoKappa_, achievedTol_;
510  int maxRestarts_, maxIters_, numIters_;
511  int blockSize_, numBlocks_, verbosity_, outputStyle_, outputFreq_, defQuorum_;
512  bool showMaxResNormOnly_;
513  std::string orthoType_;
514  std::string impResScale_, expResScale_;
515  MagnitudeType resScaleFactor_;
516 
517  // Timers.
518  std::string label_;
519  Teuchos::RCP<Teuchos::Time> timerSolve_;
520 
521  // Internal state variables.
522  bool isSet_, isSTSet_, expResTest_;
523  bool loaDetected_;
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>
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
600 {
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_);
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_);
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
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
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");
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")) {
852  MagnitudeType tempResScaleFactor = DefaultSolverParameters::resScaleFactor;
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_);
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
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
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>
1092 {
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 
1260 
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())) ) {
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 
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.
1338 
1339  // Orthonormalize the new V_0
1340  int rank = ortho_->normalize( *tmpV, tmpZ );
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.
1498 
1499  // Orthonormalize the new V_0
1500  int rank = ortho_->normalize( *tmpV, tmpZ );
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
1678  const Teuchos::EVerbosityLevel verbLevel) const
1679 {
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.
bool is_null(const boost::shared_ptr< T > &p)
ScaleType
The type of scaling to use on the residual norm value.
Definition: BelosTypes.hpp:119
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...
T & get(const std::string &name, T def_value)
PseudoBlockGmresSolMgrOrthoFailure(const std::string &what_arg)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
static RCP< Time > getNewCounter(const std::string &name)
bool is_null(const std::shared_ptr< T > &p)
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
bool isSublist(const std::string &name) const
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.
Ordinal numParams() const
An implementation of StatusTestResNorm using a family of residual norms.
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.
Iterated Modified Gram-Schmidt (IMGS) implementation of the Belos::OrthoManager class.
Traits class which defines basic operations on multivectors.
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
Classical Gram-Schmidt (with DGKS correction) implementation of the Belos::OrthoManager class...
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
Definition: BelosTypes.hpp:205
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Pure virtual base class which describes the basic interface for a solver manager. ...
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.
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
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.
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...
bool isType(const std::string &name) const
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:154
The Belos::SolverManager is a templated virtual base class that defines the basic interface that any ...
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
A list of valid default parameters for this solver.
Iterated Classical Gram-Schmidt (ICGS) implementation of the Belos::OrthoManager class.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
static const EVerbosityLevel verbLevel_default
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
bool isLOADetected() const override
Whether a "loss of accuracy" was detected during the last solve().
bool nonnull(const boost::shared_ptr< T > &p)
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, ScalarType > > > Z
The current right-hand side of the least squares system RY = Z.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem to solve.
Belos concrete class for performing the pseudo-block GMRES iteration.
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.
bool isParameter(const std::string &name) const
Orthogonalization manager based on Tall Skinny QR (TSQR)
Class which defines basic traits for the operator type.
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
int curDim
The current dimension of the reduction.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:60
ReturnType solve() override
This method performs possibly repeated calls to the underlying linear solver&#39;s iterate() routine unti...
Default parameters common to most Belos solvers.
Definition: BelosTypes.hpp: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...
This class implements the pseudo-block GMRES iteration, where a block Krylov subspace is constructed ...
Factory to build a set of status tests from a parameter list.
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
The current parameters for this solver.
PseudoBlockGmresIterOrthoFailure is thrown when the orthogonalization manager is unable to generate o...
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Return a reference to the linear problem being solved by this solver manager.
PseudoBlockGmresSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate...
virtual void setUserConvStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &userConvStatusTest, const typename StatusTestCombo< ScalarType, MV, OP >::ComboType &comboType=StatusTestCombo< ScalarType, MV, OP >::SEQ) override
Set a custom status test.
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.

Generated for Belos by doxygen 1.8.14