42 #ifndef BELOS_LSQR_SOLMGR_HPP 43 #define BELOS_LSQR_SOLMGR_HPP 61 #include "Teuchos_as.hpp" 63 #ifdef BELOS_TEUCHOS_TIME_MONITOR 64 #include "Teuchos_TimeMonitor.hpp" 217 template<
class ScalarType,
class MV,
class OP,
218 const bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex>
221 Teuchos::ScalarTraits<ScalarType>::isComplex>
223 static const bool isComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
231 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
237 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
246 template<
class ScalarType,
class MV,
class OP>
252 typedef Teuchos::ScalarTraits<ScalarType> STS;
253 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
254 typedef Teuchos::ScalarTraits<MagnitudeType> STM;
297 const Teuchos::RCP<Teuchos::ParameterList>& pl);
303 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
318 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters()
const override;
331 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers ()
const {
332 return Teuchos::tuple (timerSolve_);
394 void setParameters (
const Teuchos::RCP<Teuchos::ParameterList>& params)
override;
406 problem_->setProblem ();
439 std::string description ()
const override;
446 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
448 Teuchos::RCP<OutputManager<ScalarType> > printer_;
450 Teuchos::RCP<std::ostream> outputStream_;
453 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
454 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
455 Teuchos::RCP<LSQRStatusTest<ScalarType,MV,OP> > convTest_;
456 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
459 Teuchos::RCP<Teuchos::ParameterList> params_;
466 mutable Teuchos::RCP<const Teuchos::ParameterList> validParams_;
469 MagnitudeType lambda_;
470 MagnitudeType relRhsErr_;
471 MagnitudeType relMatErr_;
472 MagnitudeType condMax_;
473 int maxIters_, termIterMax_;
474 int verbosity_, outputStyle_, outputFreq_;
478 MagnitudeType matCondNum_;
479 MagnitudeType matNorm_;
480 MagnitudeType resNorm_;
481 MagnitudeType matResNorm_;
485 Teuchos::RCP<Teuchos::Time> timerSolve_;
492 template<
class ScalarType,
class MV,
class OP>
494 lambda_ (STM::zero ()),
495 relRhsErr_ (Teuchos::as<MagnitudeType> (10) * STM::squareroot (STM::eps ())),
496 relMatErr_ (Teuchos::as<MagnitudeType> (10) * STM::squareroot (STM::eps ())),
497 condMax_ (STM::one () / STM::eps ()),
504 matCondNum_ (STM::zero ()),
505 matNorm_ (STM::zero ()),
506 resNorm_ (STM::zero ()),
507 matResNorm_ (STM::zero ()),
512 template<
class ScalarType,
class MV,
class OP>
515 const Teuchos::RCP<Teuchos::ParameterList>& pl) :
517 lambda_ (STM::zero ()),
518 relRhsErr_ (Teuchos::as<MagnitudeType> (10) * STM::squareroot (STM::eps ())),
519 relMatErr_ (Teuchos::as<MagnitudeType> (10) * STM::squareroot (STM::eps ())),
520 condMax_ (STM::one () / STM::eps ()),
527 matCondNum_ (STM::zero ()),
528 matNorm_ (STM::zero ()),
529 resNorm_ (STM::zero ()),
530 matResNorm_ (STM::zero ()),
541 if (! pl.is_null ()) {
547 template<
class ScalarType,
class MV,
class OP>
548 Teuchos::RCP<const Teuchos::ParameterList>
551 using Teuchos::ParameterList;
552 using Teuchos::parameterList;
555 using Teuchos::rcpFromRef;
558 if (validParams_.is_null ()) {
562 const MagnitudeType ten = Teuchos::as<MagnitudeType> (10);
563 const MagnitudeType sqrtEps = STM::squareroot (STM::eps());
565 const MagnitudeType lambda = STM::zero();
566 RCP<std::ostream> outputStream = rcpFromRef (std::cout);
567 const MagnitudeType relRhsErr = ten * sqrtEps;
568 const MagnitudeType relMatErr = ten * sqrtEps;
569 const MagnitudeType condMax = STM::one() / STM::eps();
570 const int maxIters = 1000;
571 const int termIterMax = 1;
574 const int outputFreq = -1;
575 const std::string label (
"Belos");
577 RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
578 pl->set (
"Output Stream", outputStream,
"Teuchos::RCP<std::ostream> " 579 "(reference-counted pointer to the output stream) receiving " 580 "all solver output");
581 pl->set (
"Lambda", lambda,
"Damping parameter");
582 pl->set (
"Rel RHS Err", relRhsErr,
"Estimates the error in the data " 583 "defining the right-hand side");
584 pl->set (
"Rel Mat Err", relMatErr,
"Estimates the error in the data " 585 "defining the matrix.");
586 pl->set (
"Condition Limit", condMax,
"Bounds the estimated condition " 588 pl->set (
"Maximum Iterations", maxIters,
"Maximum number of iterations");
589 pl->set (
"Term Iter Max", termIterMax,
"The number of consecutive " 590 "iterations must that satisfy all convergence criteria in order " 591 "for LSQR to stop iterating");
592 pl->set (
"Verbosity", verbosity,
"Type(s) of solver information written to " 593 "the output stream");
594 pl->set (
"Output Style", outputStyle,
"Style of solver output");
595 pl->set (
"Output Frequency", outputFreq,
"Frequency at which information " 596 "is written to the output stream (-1 means \"not at all\")");
597 pl->set (
"Timer Label", label,
"String to use as a prefix for the timer " 599 pl->set (
"Block Size", 1,
"Block size parameter (currently, LSQR requires " 600 "this must always be 1)");
607 template<
class ScalarType,
class MV,
class OP>
612 using Teuchos::isParameterType;
613 using Teuchos::getParameter;
615 using Teuchos::ParameterList;
616 using Teuchos::parameterList;
619 using Teuchos::rcp_dynamic_cast;
620 using Teuchos::rcpFromRef;
622 using Teuchos::TimeMonitor;
623 using Teuchos::Exceptions::InvalidParameter;
624 using Teuchos::Exceptions::InvalidParameterName;
625 using Teuchos::Exceptions::InvalidParameterType;
627 TEUCHOS_TEST_FOR_EXCEPTION
628 (params.is_null (), std::invalid_argument,
629 "Belos::LSQRSolMgr::setParameters: The input ParameterList is null.");
630 RCP<const ParameterList> defaultParams = getValidParameters ();
651 if (params->isParameter (
"Lambda")) {
652 lambda_ = params->get<MagnitudeType> (
"Lambda");
653 }
else if (params->isParameter (
"lambda")) {
654 lambda_ = params->get<MagnitudeType> (
"lambda");
658 if (params->isParameter (
"Maximum Iterations")) {
659 maxIters_ = params->get<
int> (
"Maximum Iterations");
661 TEUCHOS_TEST_FOR_EXCEPTION
662 (maxIters_ < 0, std::invalid_argument,
"Belos::LSQRSolMgr::setParameters: " 663 "\"Maximum Iterations\" = " << maxIters_ <<
" < 0.");
667 const std::string newLabel =
668 params->isParameter (
"Timer Label") ?
669 params->get<std::string> (
"Timer Label") :
673 if (newLabel != label_) {
677 #ifdef BELOS_TEUCHOS_TIME_MONITOR 678 const std::string newSolveLabel = (newLabel !=
"") ?
679 (newLabel +
": Belos::LSQRSolMgr total solve time") :
680 std::string (
"Belos::LSQRSolMgr total solve time");
681 if (timerSolve_.is_null ()) {
683 timerSolve_ = TimeMonitor::getNewCounter (newSolveLabel);
692 const std::string oldSolveLabel = timerSolve_->name ();
694 if (oldSolveLabel != newSolveLabel) {
697 TimeMonitor::clearCounter (oldSolveLabel);
698 timerSolve_ = TimeMonitor::getNewCounter (newSolveLabel);
701 #endif // BELOS_TEUCHOS_TIME_MONITOR 705 if (params->isParameter (
"Verbosity")) {
706 int newVerbosity = 0;
714 }
catch (Teuchos::Exceptions::InvalidParameterType&) {
715 newVerbosity = params->get<
int> (
"Verbosity");
717 if (newVerbosity != verbosity_) {
718 verbosity_ = newVerbosity;
723 if (params->isParameter (
"Output Style")) {
724 outputStyle_ = params->get<
int> (
"Output Style");
731 if (params->isParameter (
"Output Stream")) {
732 outputStream_ = params->get<RCP<std::ostream> > (
"Output Stream");
739 if (outputStream_.is_null ()) {
740 outputStream_ = rcp (
new Teuchos::oblackholestream ());
745 if (params->isParameter (
"Output Frequency")) {
746 outputFreq_ = params->get<
int> (
"Output Frequency");
752 if (printer_.is_null ()) {
755 printer_->setVerbosity (verbosity_);
756 printer_->setOStream (outputStream_);
763 if (params->isParameter (
"Condition Limit")) {
764 condMax_ = params->get<MagnitudeType> (
"Condition Limit");
766 if (params->isParameter (
"Term Iter Max")) {
767 termIterMax_ = params->get<
int> (
"Term Iter Max");
769 if (params->isParameter (
"Rel RHS Err")) {
770 relRhsErr_ = params->get<MagnitudeType> (
"Rel RHS Err");
772 else if (params->isParameter (
"Convergence Tolerance")) {
775 relRhsErr_ = params->get<MagnitudeType> (
"Convergence Tolerance");
778 if (params->isParameter (
"Rel Mat Err")) {
779 relMatErr_ = params->get<MagnitudeType> (
"Rel Mat Err");
784 if (convTest_.is_null ()) {
787 relRhsErr_, relMatErr_));
789 convTest_->setCondLim (condMax_);
790 convTest_->setTermIterMax (termIterMax_);
791 convTest_->setRelRhsErr (relRhsErr_);
792 convTest_->setRelMatErr (relMatErr_);
799 if (maxIterTest_.is_null()) {
802 maxIterTest_->setMaxIters (maxIters_);
814 if (sTest_.is_null()) {
815 sTest_ = rcp (
new combo_type (combo_type::OR, maxIterTest_, convTest_));
818 if (outputTest_.is_null ()) {
822 outputTest_ = stoFactory.
create (printer_, sTest_, outputFreq_,
825 const std::string solverDesc =
" LSQR ";
826 outputTest_->setSolverDesc (solverDesc);
830 outputTest_->setOutputManager (printer_);
831 outputTest_->setChild (sTest_);
832 outputTest_->setOutputFrequency (outputFreq_);
848 template<
class ScalarType,
class MV,
class OP>
860 this->setParameters (Teuchos::parameterList (* (getValidParameters ())));
863 TEUCHOS_TEST_FOR_EXCEPTION
865 "Belos::LSQRSolMgr::solve: The linear problem to solve is null.");
866 TEUCHOS_TEST_FOR_EXCEPTION
868 "Belos::LSQRSolMgr::solve: The linear problem is not ready, " 869 "as its setProblem() method has not been called.");
870 TEUCHOS_TEST_FOR_EXCEPTION
871 (MVT::GetNumberVecs (*(problem_->getRHS ())) != 1,
873 "The current implementation of LSQR only knows how to solve problems " 874 "with one right-hand side, but the linear problem to solve has " 875 << MVT::GetNumberVecs (* (problem_->getRHS ()))
876 <<
" right-hand sides.");
892 std::vector<int> currRHSIdx (1, 0);
893 problem_->setLSIndex (currRHSIdx);
896 outputTest_->reset ();
900 bool isConverged =
false;
917 Teuchos::ParameterList plist;
924 plist.set (
"Lambda", lambda_);
927 RCP<iter_type> lsqr_iter =
928 rcp (
new iter_type (problem_, printer_, outputTest_, plist));
929 #ifdef BELOS_TEUCHOS_TIME_MONITOR 930 Teuchos::TimeMonitor slvtimer (*timerSolve_);
934 lsqr_iter->resetNumIters ();
936 outputTest_->resetNumCalls ();
939 lsqr_iter->initializeLSQR (newstate);
942 lsqr_iter->iterate ();
952 TEUCHOS_TEST_FOR_EXCEPTION
953 (
true, std::logic_error,
"Belos::LSQRSolMgr::solve: " 954 "LSQRIteration::iterate returned without either the convergence test " 955 "or the maximum iteration count test passing. " 956 "Please report this bug to the Belos developers.");
958 }
catch (
const std::exception& e) {
960 <<
"Error! Caught std::exception in LSQRIter::iterate at iteration " 961 << lsqr_iter->getNumIters () << std::endl << e.what () << std::endl;
966 problem_->setCurrLS();
972 #ifdef BELOS_TEUCHOS_TIME_MONITOR 978 #endif // BELOS_TEUCHOS_TIME_MONITOR 981 numIters_ = maxIterTest_->getNumIters();
982 matCondNum_ = convTest_->getMatCondNum();
983 matNorm_ = convTest_->getMatNorm();
984 resNorm_ = convTest_->getResidNorm();
985 matResNorm_ = convTest_->getLSResidNorm();
995 template<
class ScalarType,
class MV,
class OP>
998 std::ostringstream oss;
999 oss <<
"LSQRSolMgr<...," << STS::name () <<
">";
1001 oss <<
"Lambda: " << lambda_;
1002 oss <<
", condition number limit: " << condMax_;
1003 oss <<
", relative RHS Error: " << relRhsErr_;
1004 oss <<
", relative Matrix Error: " << relMatErr_;
1005 oss <<
", maximum number of iterations: " << maxIters_;
1006 oss <<
", termIterMax: " << termIterMax_;
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
void reset(const ResetType type) override
reset the solver manager as specified by the ResetType, informs the solver manager that the solver sh...
Collection of types and exceptions used within the Belos solvers.
IterationState contains the data that defines the state of the LSQR solver at any given time...
Belos concrete class that iterates LSQR.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
MagnitudeType getResNorm() const
Estimated residual norm from the last solve.
Class which manages the output and verbosity of the Belos solvers.
LSQRSolMgrBlockSizeFailure(const std::string &what_arg)
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
LSQRSolMgrBlockSizeFailure is thrown when the linear problem has more than one RHS.
MsgType
Available message types recognized by the linear solvers.
A factory class for generating StatusTestOutput objects.
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
bool isLOADetected() const override
Whether a loss of accuracy was detected during the last solve.
Base class for Belos::SolverManager subclasses which normally can only compile for real ScalarType...
Belos::StatusTest class for specifying a maximum number of iterations.
Belos::StatusTest class defining LSQR convergence.
A factory class for generating StatusTestOutput objects.
LSQRSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
Belos::LSQRSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i...
A Belos::StatusTest class for specifying a maximum number of iterations.
int getNumIters() const override
Iteration count from the last solve.
ResetType
How to reset the solver.
Pure virtual base class which describes the basic interface for a solver manager. ...
MagnitudeType getMatResNorm() const
Estimate of (residual vector ) from the last solve.
A linear system to solve, and its associated information.
virtual ~LSQRSolMgr()
Destructor (declared virtual for memory safety of base classes).
Class which describes the linear problem to be solved by the iterative solver.
ReturnType
Whether the Belos solve converged for all linear systems.
MagnitudeType getMatNorm() const
Estimated matrix Frobenius norm from the last solve.
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
MagnitudeType getMatCondNum() const
Estimated matrix condition number from the last solve.
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.
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
A class for extending the status testing capabilities of Belos via logical combinations.
LSQR method (for linear systems and linear least-squares problems).
Class which defines basic traits for the operator type.
Implementation of the LSQR iteration.
Parent class to all Belos exceptions.
Structure to contain pointers to LSQRIteration state variables, ...
Belos header file which uses auto-configuration information to include necessary C++ headers...
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Get current linear problem being solved for in this object.
LSQRSolMgrLinearProblemFailure(const std::string &what_arg)