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) :
258 template<
class ScalarType,
class MV,
class OP>
264 typedef Teuchos::ScalarTraits<ScalarType> STS;
265 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
266 typedef Teuchos::ScalarTraits<MagnitudeType> STM;
309 const Teuchos::RCP<Teuchos::ParameterList>& pl);
326 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters()
const;
339 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers ()
const {
340 return Teuchos::tuple (timerSolve_);
402 void setParameters (
const Teuchos::RCP<Teuchos::ParameterList>& params);
414 problem_->setProblem ();
447 std::string description ()
const;
454 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
456 Teuchos::RCP<OutputManager<ScalarType> > printer_;
458 Teuchos::RCP<std::ostream> outputStream_;
461 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
462 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
463 Teuchos::RCP<LSQRStatusTest<ScalarType,MV,OP> > convTest_;
464 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
467 Teuchos::RCP<Teuchos::ParameterList> params_;
474 mutable Teuchos::RCP<const Teuchos::ParameterList> validParams_;
477 MagnitudeType lambda_;
478 MagnitudeType relRhsErr_;
479 MagnitudeType relMatErr_;
480 MagnitudeType condMax_;
481 int maxIters_, termIterMax_;
482 int verbosity_, outputStyle_, outputFreq_;
486 MagnitudeType matCondNum_;
487 MagnitudeType matNorm_;
488 MagnitudeType resNorm_;
489 MagnitudeType matResNorm_;
493 Teuchos::RCP<Teuchos::Time> timerSolve_;
500 template<
class ScalarType,
class MV,
class OP>
502 lambda_ (STM::zero ()),
503 relRhsErr_ (Teuchos::as<MagnitudeType> (10) * STM::squareroot (STM::eps ())),
504 relMatErr_ (Teuchos::as<MagnitudeType> (10) * STM::squareroot (STM::eps ())),
505 condMax_ (STM::one () / STM::eps ()),
512 matCondNum_ (STM::zero ()),
513 matNorm_ (STM::zero ()),
514 resNorm_ (STM::zero ()),
515 matResNorm_ (STM::zero ()),
520 template<
class ScalarType,
class MV,
class OP>
523 const Teuchos::RCP<Teuchos::ParameterList>& pl) :
525 lambda_ (STM::zero ()),
526 relRhsErr_ (Teuchos::as<MagnitudeType> (10) * STM::squareroot (STM::eps ())),
527 relMatErr_ (Teuchos::as<MagnitudeType> (10) * STM::squareroot (STM::eps ())),
528 condMax_ (STM::one () / STM::eps ()),
535 matCondNum_ (STM::zero ()),
536 matNorm_ (STM::zero ()),
537 resNorm_ (STM::zero ()),
538 matResNorm_ (STM::zero ()),
549 if (! pl.is_null ()) {
555 template<
class ScalarType,
class MV,
class OP>
556 Teuchos::RCP<const Teuchos::ParameterList>
559 using Teuchos::ParameterList;
560 using Teuchos::parameterList;
563 using Teuchos::rcpFromRef;
566 if (validParams_.is_null ()) {
570 const MagnitudeType ten = Teuchos::as<MagnitudeType> (10);
571 const MagnitudeType sqrtEps = STM::squareroot (STM::eps());
573 const MagnitudeType lambda = STM::zero();
574 RCP<std::ostream> outputStream = rcpFromRef (std::cout);
575 const MagnitudeType relRhsErr = ten * sqrtEps;
576 const MagnitudeType relMatErr = ten * sqrtEps;
577 const MagnitudeType condMax = STM::one() / STM::eps();
578 const int maxIters = 1000;
579 const int termIterMax = 1;
582 const int outputFreq = -1;
583 const std::string label (
"Belos");
585 RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
586 pl->set (
"Output Stream", outputStream,
"Teuchos::RCP<std::ostream> " 587 "(reference-counted pointer to the output stream) receiving " 588 "all solver output");
589 pl->set (
"Lambda", lambda,
"Damping parameter");
590 pl->set (
"Rel RHS Err", relRhsErr,
"Estimates the error in the data " 591 "defining the right-hand side");
592 pl->set (
"Rel Mat Err", relMatErr,
"Estimates the error in the data " 593 "defining the matrix.");
594 pl->set (
"Condition Limit", condMax,
"Bounds the estimated condition " 596 pl->set (
"Maximum Iterations", maxIters,
"Maximum number of iterations");
597 pl->set (
"Term Iter Max", termIterMax,
"The number of consecutive " 598 "iterations must that satisfy all convergence criteria in order " 599 "for LSQR to stop iterating");
600 pl->set (
"Verbosity", verbosity,
"Type(s) of solver information written to " 601 "the output stream");
602 pl->set (
"Output Style", outputStyle,
"Style of solver output");
603 pl->set (
"Output Frequency", outputFreq,
"Frequency at which information " 604 "is written to the output stream (-1 means \"not at all\")");
605 pl->set (
"Timer Label", label,
"String to use as a prefix for the timer " 608 pl->set (
"Block Size", 1,
"Block size parameter (currently, LSQR requires " 609 "this must always be 1)");
616 template<
class ScalarType,
class MV,
class OP>
621 using Teuchos::isParameterType;
622 using Teuchos::getParameter;
624 using Teuchos::ParameterList;
625 using Teuchos::parameterList;
628 using Teuchos::rcp_dynamic_cast;
629 using Teuchos::rcpFromRef;
631 using Teuchos::TimeMonitor;
632 using Teuchos::Exceptions::InvalidParameter;
633 using Teuchos::Exceptions::InvalidParameterName;
634 using Teuchos::Exceptions::InvalidParameterType;
636 TEUCHOS_TEST_FOR_EXCEPTION
637 (params.is_null (), std::invalid_argument,
638 "Belos::LSQRSolMgr::setParameters: The input ParameterList is null.");
639 RCP<const ParameterList> defaultParams = getValidParameters ();
660 if (params->isParameter (
"Lambda")) {
661 lambda_ = params->get<MagnitudeType> (
"Lambda");
662 }
else if (params->isParameter (
"lambda")) {
663 lambda_ = params->get<MagnitudeType> (
"lambda");
667 if (params->isParameter (
"Maximum Iterations")) {
668 maxIters_ = params->get<
int> (
"Maximum Iterations");
670 TEUCHOS_TEST_FOR_EXCEPTION
671 (maxIters_ < 0, std::invalid_argument,
"Belos::LSQRSolMgr::setParameters: " 672 "\"Maximum Iterations\" = " << maxIters_ <<
" < 0.");
676 const std::string newLabel =
677 params->isParameter (
"Timer Label") ?
678 params->get<std::string> (
"Timer Label") :
682 if (newLabel != label_) {
686 #ifdef BELOS_TEUCHOS_TIME_MONITOR 687 const std::string newSolveLabel = (newLabel !=
"") ?
688 (newLabel +
": Belos::LSQRSolMgr total solve time") :
689 std::string (
"Belos::LSQRSolMgr total solve time");
690 if (timerSolve_.is_null ()) {
692 timerSolve_ = TimeMonitor::getNewCounter (newSolveLabel);
701 const std::string oldSolveLabel = timerSolve_->name ();
703 if (oldSolveLabel != newSolveLabel) {
706 TimeMonitor::clearCounter (oldSolveLabel);
707 timerSolve_ = TimeMonitor::getNewCounter (newSolveLabel);
710 #endif // BELOS_TEUCHOS_TIME_MONITOR 714 if (params->isParameter (
"Verbosity")) {
715 int newVerbosity = 0;
723 }
catch (Teuchos::Exceptions::InvalidParameterType&) {
724 newVerbosity = params->get<
int> (
"Verbosity");
726 if (newVerbosity != verbosity_) {
727 verbosity_ = newVerbosity;
732 if (params->isParameter (
"Output Style")) {
733 outputStyle_ = params->get<
int> (
"Output Style");
740 if (params->isParameter (
"Output Stream")) {
741 outputStream_ = params->get<RCP<std::ostream> > (
"Output Stream");
748 if (outputStream_.is_null ()) {
749 outputStream_ = rcp (
new Teuchos::oblackholestream ());
754 if (params->isParameter (
"Output Frequency")) {
755 outputFreq_ = params->get<
int> (
"Output Frequency");
761 if (printer_.is_null ()) {
764 printer_->setVerbosity (verbosity_);
765 printer_->setOStream (outputStream_);
772 if (params->isParameter (
"Condition Limit")) {
773 condMax_ = params->get<MagnitudeType> (
"Condition Limit");
775 if (params->isParameter (
"Term Iter Max")) {
776 termIterMax_ = params->get<
int> (
"Term Iter Max");
778 if (params->isParameter (
"Rel RHS Err")) {
779 relRhsErr_ = params->get<MagnitudeType> (
"Rel RHS Err");
781 else if (params->isParameter (
"Convergence Tolerance")) {
784 relRhsErr_ = params->get<MagnitudeType> (
"Convergence Tolerance");
787 if (params->isParameter (
"Rel Mat Err")) {
788 relMatErr_ = params->get<MagnitudeType> (
"Rel Mat Err");
793 if (convTest_.is_null ()) {
796 relRhsErr_, relMatErr_));
798 convTest_->setCondLim (condMax_);
799 convTest_->setTermIterMax (termIterMax_);
800 convTest_->setRelRhsErr (relRhsErr_);
801 convTest_->setRelMatErr (relMatErr_);
808 if (maxIterTest_.is_null()) {
811 maxIterTest_->setMaxIters (maxIters_);
823 if (sTest_.is_null()) {
824 sTest_ = rcp (
new combo_type (combo_type::OR, maxIterTest_, convTest_));
827 if (outputTest_.is_null ()) {
831 outputTest_ = stoFactory.
create (printer_, sTest_, outputFreq_,
834 const std::string solverDesc =
" LSQR ";
835 outputTest_->setSolverDesc (solverDesc);
839 outputTest_->setOutputManager (printer_);
840 outputTest_->setChild (sTest_);
841 outputTest_->setOutputFrequency (outputFreq_);
857 template<
class ScalarType,
class MV,
class OP>
869 this->setParameters (Teuchos::parameterList (* (getValidParameters ())));
872 TEUCHOS_TEST_FOR_EXCEPTION
874 "Belos::LSQRSolMgr::solve: The linear problem to solve is null.");
875 TEUCHOS_TEST_FOR_EXCEPTION
877 "Belos::LSQRSolMgr::solve: The linear problem is not ready, " 878 "as its setProblem() method has not been called.");
879 TEUCHOS_TEST_FOR_EXCEPTION
880 (MVT::GetNumberVecs (*(problem_->getRHS ())) != 1,
882 "The current implementation of LSQR only knows how to solve problems " 883 "with one right-hand side, but the linear problem to solve has " 884 << MVT::GetNumberVecs (* (problem_->getRHS ()))
885 <<
" right-hand sides.");
901 std::vector<int> currRHSIdx (1, 0);
902 problem_->setLSIndex (currRHSIdx);
905 outputTest_->reset ();
909 bool isConverged =
false;
926 Teuchos::ParameterList plist;
933 plist.set (
"Lambda", lambda_);
936 RCP<iter_type> lsqr_iter =
937 rcp (
new iter_type (problem_, printer_, outputTest_, plist));
938 #ifdef BELOS_TEUCHOS_TIME_MONITOR 939 Teuchos::TimeMonitor slvtimer (*timerSolve_);
943 lsqr_iter->resetNumIters ();
945 outputTest_->resetNumCalls ();
948 lsqr_iter->initializeLSQR (newstate);
951 lsqr_iter->iterate ();
961 TEUCHOS_TEST_FOR_EXCEPTION
962 (
true, std::logic_error,
"Belos::LSQRSolMgr::solve: " 963 "LSQRIteration::iterate returned without either the convergence test " 964 "or the maximum iteration count test passing. " 965 "Please report this bug to the Belos developers.");
967 }
catch (
const std::exception& e) {
969 <<
"Error! Caught std::exception in LSQRIter::iterate at iteration " 970 << lsqr_iter->getNumIters () << std::endl << e.what () << std::endl;
975 problem_->setCurrLS();
981 #ifdef BELOS_TEUCHOS_TIME_MONITOR 987 #endif // BELOS_TEUCHOS_TIME_MONITOR 990 numIters_ = maxIterTest_->getNumIters();
991 matCondNum_ = convTest_->getMatCondNum();
992 matNorm_ = convTest_->getMatNorm();
993 resNorm_ = convTest_->getResidNorm();
994 matResNorm_ = convTest_->getLSResidNorm();
1004 template<
class ScalarType,
class MV,
class OP>
1007 std::ostringstream oss;
1008 oss <<
"LSQRSolMgr<...," << STS::name () <<
">";
1010 oss <<
"Lambda: " << lambda_;
1011 oss <<
", condition number limit: " << condMax_;
1012 oss <<
", relative RHS Error: " << relRhsErr_;
1013 oss <<
", relative Matrix Error: " << relMatErr_;
1014 oss <<
", maximum number of iterations: " << maxIters_;
1015 oss <<
", termIterMax: " << termIterMax_;
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
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...
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
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)
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.
int getNumIters() const
Iteration count from 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.
void reset(const ResetType type)
reset the solver manager as specified by the ResetType, informs the solver manager that the solver sh...
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem)
Set the linear problem that needs to be solved.
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.
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.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get current linear problem being solved for in this object.
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.
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...
bool isLOADetected() const
Whether a loss of accuracy was detected during the last solve.
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...
LSQRSolMgrLinearProblemFailure(const std::string &what_arg)