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"
217template<
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 {
246template<
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_;
453 Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
sTest_;
455 Teuchos::RCP<LSQRStatusTest<ScalarType,MV,OP> >
convTest_;
492template<
class ScalarType,
class MV,
class OP>
494 lambda_ (
STM::zero ()),
497 condMax_ (
STM::one () /
STM::eps ()),
504 matCondNum_ (
STM::zero ()),
505 matNorm_ (
STM::zero ()),
506 resNorm_ (
STM::zero ()),
507 matResNorm_ (
STM::zero ()),
512template<
class ScalarType,
class MV,
class OP>
515 const Teuchos::RCP<Teuchos::ParameterList>& pl) :
517 lambda_ (
STM::zero ()),
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 ()) {
547template<
class ScalarType,
class MV,
class OP>
548Teuchos::RCP<const Teuchos::ParameterList>
551 using Teuchos::ParameterList;
552 using Teuchos::parameterList;
555 using Teuchos::rcpFromRef;
558 if (validParams_.is_null ()) {
566 RCP<std::ostream> outputStream = rcpFromRef (std::cout);
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)");
607template<
class ScalarType,
class MV,
class OP>
610setParameters (
const Teuchos::RCP<Teuchos::ParameterList>& params)
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")) {
653 }
else if (params->isParameter (
"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);
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")) {
766 if (params->isParameter (
"Term Iter Max")) {
767 termIterMax_ = params->get<
int> (
"Term Iter Max");
769 if (params->isParameter (
"Rel RHS Err")) {
772 else if (params->isParameter (
"Convergence Tolerance")) {
775 relRhsErr_ = params->get<
MagnitudeType> (
"Convergence Tolerance");
778 if (params->isParameter (
"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_);
848template<
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
981 numIters_ = maxIterTest_->getNumIters();
982 matCondNum_ = convTest_->getMatCondNum();
983 matNorm_ = convTest_->getMatNorm();
984 resNorm_ = convTest_->getResidNorm();
985 matResNorm_ = convTest_->getLSResidNorm();
995template<
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_;
Belos header file which uses auto-configuration information to include necessary C++ headers.
Belos concrete class that iterates LSQR.
IterationState contains the data that defines the state of the LSQR solver at any given time.
Belos::StatusTest class defining LSQR convergence.
Class which describes the linear problem to be solved by the iterative solver.
Class which manages the output and verbosity of the Belos solvers.
Pure virtual base class which describes the basic interface for a solver manager.
Belos::StatusTest for logically combining several status tests.
Belos::StatusTest class for specifying a maximum number of iterations.
A factory class for generating StatusTestOutput objects.
Collection of types and exceptions used within the Belos solvers.
Parent class to all Belos exceptions.
Base class for Belos::SolverManager subclasses which normally can only compile for real ScalarType.
Implementation of the LSQR iteration.
LSQRSolMgrBlockSizeFailure is thrown when the linear problem has more than one RHS.
LSQRSolMgrBlockSizeFailure(const std::string &what_arg)
Belos::LSQRSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i....
LSQRSolMgrLinearProblemFailure(const std::string &what_arg)
MagnitudeType matResNorm_
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
The "master" status test (that includes all status tests).
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Get current linear problem being solved for in this object.
Teuchos::ScalarTraits< ScalarType > STS
MagnitudeType getMatCondNum() const
Estimated matrix condition number from the last solve.
MultiVecTraits< ScalarType, MV > MVT
Teuchos::RCP< Teuchos::Time > timerSolve_
Teuchos::ScalarTraits< MagnitudeType > STM
Teuchos::RCP< OutputManager< ScalarType > > printer_
The output manager.
int getNumIters() const override
Iteration count from the last solve.
virtual ~LSQRSolMgr()
Destructor (declared virtual for memory safety of base classes).
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...
Teuchos::RCP< std::ostream > outputStream_
Output stream to which to write status output.
bool isLOADetected() const override
Whether a loss of accuracy was detected during the last solve.
MagnitudeType getMatResNorm() const
Estimate of (residual vector ) from the last solve.
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::RCP< LSQRStatusTest< ScalarType, MV, OP > > convTest_
MagnitudeType getResNorm() const
Estimated residual norm from the last solve.
Teuchos::RCP< Teuchos::ParameterList > params_
Current parameter list.
MagnitudeType getMatNorm() const
Estimated matrix Frobenius norm from the last solve.
MagnitudeType matCondNum_
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
Teuchos::RCP< const Teuchos::ParameterList > validParams_
Default parameter list.
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
The linear problem to solve.
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
LSQR method (for linear systems and linear least-squares problems).
Details::RealSolverManager< ScalarType, MV, OP, isComplex > base_type
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
static const bool isComplex
LSQRSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
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.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
A class for extending the status testing capabilities of Belos via logical combinations.
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.
MsgType
Available message types recognized by the linear solvers.
ReturnType
Whether the Belos solve converged for all linear systems.
ResetType
How to reset the solver.
Structure to contain pointers to LSQRIteration state variables, ...