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

Generated for Belos by doxygen 1.8.14