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