42 #ifndef BELOS_LSQR_SOLMGR_HPP 43 #define BELOS_LSQR_SOLMGR_HPP 61 #include "Teuchos_as.hpp" 62 #include "Teuchos_BLAS.hpp" 63 #include "Teuchos_LAPACK.hpp" 65 #ifdef BELOS_TEUCHOS_TIME_MONITOR 66 #include "Teuchos_TimeMonitor.hpp" 234 template<
class ScalarType,
class MV,
class OP,
235 const bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex>
238 Teuchos::ScalarTraits<ScalarType>::isComplex>
240 static const bool isComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
248 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
254 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
263 template<
class ScalarType,
class MV,
class OP>
269 typedef Teuchos::ScalarTraits<ScalarType> STS;
270 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
271 typedef Teuchos::ScalarTraits<MagnitudeType> STM;
314 const Teuchos::RCP<Teuchos::ParameterList>& pl);
320 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
335 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters()
const override;
348 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers ()
const {
349 return Teuchos::tuple (timerSolve_);
411 void setParameters (
const Teuchos::RCP<Teuchos::ParameterList>& params)
override;
423 problem_->setProblem ();
456 std::string description ()
const override;
463 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
465 Teuchos::RCP<OutputManager<ScalarType> > printer_;
467 Teuchos::RCP<std::ostream> outputStream_;
470 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
471 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
472 Teuchos::RCP<LSQRStatusTest<ScalarType,MV,OP> > convTest_;
473 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
476 Teuchos::RCP<Teuchos::ParameterList> params_;
483 mutable Teuchos::RCP<const Teuchos::ParameterList> validParams_;
486 MagnitudeType lambda_;
487 MagnitudeType relRhsErr_;
488 MagnitudeType relMatErr_;
489 MagnitudeType condMax_;
490 int maxIters_, termIterMax_;
491 int verbosity_, outputStyle_, outputFreq_;
495 MagnitudeType matCondNum_;
496 MagnitudeType matNorm_;
497 MagnitudeType resNorm_;
498 MagnitudeType matResNorm_;
502 Teuchos::RCP<Teuchos::Time> timerSolve_;
509 template<
class ScalarType,
class MV,
class OP>
511 lambda_ (STM::zero ()),
512 relRhsErr_ (Teuchos::as<MagnitudeType> (10) * STM::squareroot (STM::eps ())),
513 relMatErr_ (Teuchos::as<MagnitudeType> (10) * STM::squareroot (STM::eps ())),
514 condMax_ (STM::one () / STM::eps ()),
521 matCondNum_ (STM::zero ()),
522 matNorm_ (STM::zero ()),
523 resNorm_ (STM::zero ()),
524 matResNorm_ (STM::zero ()),
529 template<
class ScalarType,
class MV,
class OP>
532 const Teuchos::RCP<Teuchos::ParameterList>& pl) :
534 lambda_ (STM::zero ()),
535 relRhsErr_ (Teuchos::as<MagnitudeType> (10) * STM::squareroot (STM::eps ())),
536 relMatErr_ (Teuchos::as<MagnitudeType> (10) * STM::squareroot (STM::eps ())),
537 condMax_ (STM::one () / STM::eps ()),
544 matCondNum_ (STM::zero ()),
545 matNorm_ (STM::zero ()),
546 resNorm_ (STM::zero ()),
547 matResNorm_ (STM::zero ()),
558 if (! pl.is_null ()) {
564 template<
class ScalarType,
class MV,
class OP>
565 Teuchos::RCP<const Teuchos::ParameterList>
568 using Teuchos::ParameterList;
569 using Teuchos::parameterList;
572 using Teuchos::rcpFromRef;
575 if (validParams_.is_null ()) {
579 const MagnitudeType ten = Teuchos::as<MagnitudeType> (10);
580 const MagnitudeType sqrtEps = STM::squareroot (STM::eps());
582 const MagnitudeType lambda = STM::zero();
583 RCP<std::ostream> outputStream = rcpFromRef (std::cout);
584 const MagnitudeType relRhsErr = ten * sqrtEps;
585 const MagnitudeType relMatErr = ten * sqrtEps;
586 const MagnitudeType condMax = STM::one() / STM::eps();
587 const int maxIters = 1000;
588 const int termIterMax = 1;
591 const int outputFreq = -1;
592 const std::string label (
"Belos");
594 RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
595 pl->set (
"Output Stream", outputStream,
"Teuchos::RCP<std::ostream> " 596 "(reference-counted pointer to the output stream) receiving " 597 "all solver output");
598 pl->set (
"Lambda", lambda,
"Damping parameter");
599 pl->set (
"Rel RHS Err", relRhsErr,
"Estimates the error in the data " 600 "defining the right-hand side");
601 pl->set (
"Rel Mat Err", relMatErr,
"Estimates the error in the data " 602 "defining the matrix.");
603 pl->set (
"Condition Limit", condMax,
"Bounds the estimated condition " 605 pl->set (
"Maximum Iterations", maxIters,
"Maximum number of iterations");
606 pl->set (
"Term Iter Max", termIterMax,
"The number of consecutive " 607 "iterations must that satisfy all convergence criteria in order " 608 "for LSQR to stop iterating");
609 pl->set (
"Verbosity", verbosity,
"Type(s) of solver information written to " 610 "the output stream");
611 pl->set (
"Output Style", outputStyle,
"Style of solver output");
612 pl->set (
"Output Frequency", outputFreq,
"Frequency at which information " 613 "is written to the output stream (-1 means \"not at all\")");
614 pl->set (
"Timer Label", label,
"String to use as a prefix for the timer " 616 pl->set (
"Block Size", 1,
"Block size parameter (currently, LSQR requires " 617 "this must always be 1)");
624 template<
class ScalarType,
class MV,
class OP>
629 using Teuchos::isParameterType;
630 using Teuchos::getParameter;
632 using Teuchos::ParameterList;
633 using Teuchos::parameterList;
636 using Teuchos::rcp_dynamic_cast;
637 using Teuchos::rcpFromRef;
639 using Teuchos::TimeMonitor;
640 using Teuchos::Exceptions::InvalidParameter;
641 using Teuchos::Exceptions::InvalidParameterName;
642 using Teuchos::Exceptions::InvalidParameterType;
644 TEUCHOS_TEST_FOR_EXCEPTION
645 (params.is_null (), std::invalid_argument,
646 "Belos::LSQRSolMgr::setParameters: The input ParameterList is null.");
647 RCP<const ParameterList> defaultParams = getValidParameters ();
668 if (params->isParameter (
"Lambda")) {
669 lambda_ = params->get<MagnitudeType> (
"Lambda");
670 }
else if (params->isParameter (
"lambda")) {
671 lambda_ = params->get<MagnitudeType> (
"lambda");
675 if (params->isParameter (
"Maximum Iterations")) {
676 maxIters_ = params->get<
int> (
"Maximum Iterations");
678 TEUCHOS_TEST_FOR_EXCEPTION
679 (maxIters_ < 0, std::invalid_argument,
"Belos::LSQRSolMgr::setParameters: " 680 "\"Maximum Iterations\" = " << maxIters_ <<
" < 0.");
684 const std::string newLabel =
685 params->isParameter (
"Timer Label") ?
686 params->get<std::string> (
"Timer Label") :
690 if (newLabel != label_) {
694 #ifdef BELOS_TEUCHOS_TIME_MONITOR 695 const std::string newSolveLabel = (newLabel !=
"") ?
696 (newLabel +
": Belos::LSQRSolMgr total solve time") :
697 std::string (
"Belos::LSQRSolMgr total solve time");
698 if (timerSolve_.is_null ()) {
700 timerSolve_ = TimeMonitor::getNewCounter (newSolveLabel);
709 const std::string oldSolveLabel = timerSolve_->name ();
711 if (oldSolveLabel != newSolveLabel) {
714 TimeMonitor::clearCounter (oldSolveLabel);
715 timerSolve_ = TimeMonitor::getNewCounter (newSolveLabel);
718 #endif // BELOS_TEUCHOS_TIME_MONITOR 722 if (params->isParameter (
"Verbosity")) {
723 int newVerbosity = 0;
731 }
catch (Teuchos::Exceptions::InvalidParameterType&) {
732 newVerbosity = params->get<
int> (
"Verbosity");
734 if (newVerbosity != verbosity_) {
735 verbosity_ = newVerbosity;
740 if (params->isParameter (
"Output Style")) {
741 outputStyle_ = params->get<
int> (
"Output Style");
748 if (params->isParameter (
"Output Stream")) {
749 outputStream_ = params->get<RCP<std::ostream> > (
"Output Stream");
756 if (outputStream_.is_null ()) {
757 outputStream_ = rcp (
new Teuchos::oblackholestream ());
762 if (params->isParameter (
"Output Frequency")) {
763 outputFreq_ = params->get<
int> (
"Output Frequency");
769 if (printer_.is_null ()) {
772 printer_->setVerbosity (verbosity_);
773 printer_->setOStream (outputStream_);
780 if (params->isParameter (
"Condition Limit")) {
781 condMax_ = params->get<MagnitudeType> (
"Condition Limit");
783 if (params->isParameter (
"Term Iter Max")) {
784 termIterMax_ = params->get<
int> (
"Term Iter Max");
786 if (params->isParameter (
"Rel RHS Err")) {
787 relRhsErr_ = params->get<MagnitudeType> (
"Rel RHS Err");
789 else if (params->isParameter (
"Convergence Tolerance")) {
792 relRhsErr_ = params->get<MagnitudeType> (
"Convergence Tolerance");
795 if (params->isParameter (
"Rel Mat Err")) {
796 relMatErr_ = params->get<MagnitudeType> (
"Rel Mat Err");
801 if (convTest_.is_null ()) {
804 relRhsErr_, relMatErr_));
806 convTest_->setCondLim (condMax_);
807 convTest_->setTermIterMax (termIterMax_);
808 convTest_->setRelRhsErr (relRhsErr_);
809 convTest_->setRelMatErr (relMatErr_);
816 if (maxIterTest_.is_null()) {
819 maxIterTest_->setMaxIters (maxIters_);
831 if (sTest_.is_null()) {
832 sTest_ = rcp (
new combo_type (combo_type::OR, maxIterTest_, convTest_));
835 if (outputTest_.is_null ()) {
839 outputTest_ = stoFactory.
create (printer_, sTest_, outputFreq_,
842 const std::string solverDesc =
" LSQR ";
843 outputTest_->setSolverDesc (solverDesc);
847 outputTest_->setOutputManager (printer_);
848 outputTest_->setChild (sTest_);
849 outputTest_->setOutputFrequency (outputFreq_);
865 template<
class ScalarType,
class MV,
class OP>
877 this->setParameters (Teuchos::parameterList (* (getValidParameters ())));
880 TEUCHOS_TEST_FOR_EXCEPTION
882 "Belos::LSQRSolMgr::solve: The linear problem to solve is null.");
883 TEUCHOS_TEST_FOR_EXCEPTION
885 "Belos::LSQRSolMgr::solve: The linear problem is not ready, " 886 "as its setProblem() method has not been called.");
887 TEUCHOS_TEST_FOR_EXCEPTION
888 (MVT::GetNumberVecs (*(problem_->getRHS ())) != 1,
890 "The current implementation of LSQR only knows how to solve problems " 891 "with one right-hand side, but the linear problem to solve has " 892 << MVT::GetNumberVecs (* (problem_->getRHS ()))
893 <<
" right-hand sides.");
909 std::vector<int> currRHSIdx (1, 0);
910 problem_->setLSIndex (currRHSIdx);
913 outputTest_->reset ();
917 bool isConverged =
false;
934 Teuchos::ParameterList plist;
941 plist.set (
"Lambda", lambda_);
944 RCP<iter_type> lsqr_iter =
945 rcp (
new iter_type (problem_, printer_, outputTest_, plist));
946 #ifdef BELOS_TEUCHOS_TIME_MONITOR 947 Teuchos::TimeMonitor slvtimer (*timerSolve_);
951 lsqr_iter->resetNumIters ();
953 outputTest_->resetNumCalls ();
956 lsqr_iter->initializeLSQR (newstate);
959 lsqr_iter->iterate ();
969 TEUCHOS_TEST_FOR_EXCEPTION
970 (
true, std::logic_error,
"Belos::LSQRSolMgr::solve: " 971 "LSQRIteration::iterate returned without either the convergence test " 972 "or the maximum iteration count test passing. " 973 "Please report this bug to the Belos developers.");
975 }
catch (
const std::exception& e) {
977 <<
"Error! Caught std::exception in LSQRIter::iterate at iteration " 978 << lsqr_iter->getNumIters () << std::endl << e.what () << std::endl;
983 problem_->setCurrLS();
989 #ifdef BELOS_TEUCHOS_TIME_MONITOR 995 #endif // BELOS_TEUCHOS_TIME_MONITOR 998 numIters_ = maxIterTest_->getNumIters();
999 matCondNum_ = convTest_->getMatCondNum();
1000 matNorm_ = convTest_->getMatNorm();
1001 resNorm_ = convTest_->getResidNorm();
1002 matResNorm_ = convTest_->getLSResidNorm();
1004 if (! isConverged) {
1012 template<
class ScalarType,
class MV,
class OP>
1015 std::ostringstream oss;
1016 oss <<
"LSQRSolMgr<...," << STS::name () <<
">";
1018 oss <<
"Lambda: " << lambda_;
1019 oss <<
", condition number limit: " << condMax_;
1020 oss <<
", relative RHS Error: " << relRhsErr_;
1021 oss <<
", relative Matrix Error: " << relMatErr_;
1022 oss <<
", maximum number of iterations: " << maxIters_;
1023 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.
LSQRSolMgrOrthoFailure(const std::string &what_arg)
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.
LSQRSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate orthonormal...
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)