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

Generated on Thu Jul 21 2016 14:43:57 for Belos by doxygen 1.8.11