42 #ifndef BELOS_MINRES_SOLMGR_HPP
43 #define BELOS_MINRES_SOLMGR_HPP
60 #include "Teuchos_BLAS.hpp"
61 #include "Teuchos_LAPACK.hpp"
62 #ifdef BELOS_TEUCHOS_TIME_MONITOR
63 #include "Teuchos_TimeMonitor.hpp"
66 #include "Teuchos_StandardParameterEntryValidators.hpp"
114 template<
class ScalarType,
class MV,
class OP>
120 typedef Teuchos::ScalarTraits<ScalarType> SCT;
121 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
122 typedef Teuchos::ScalarTraits< MagnitudeType > MT;
184 const Teuchos::RCP<Teuchos::ParameterList> ¶ms);
190 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
205 if (defaultParams_.is_null()) {
208 return defaultParams_;
225 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
226 return Teuchos::tuple (timerSolve_);
262 setParameters (
const Teuchos::RCP<Teuchos::ParameterList>& params)
override;
273 problem_->setProblem ();
311 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
314 Teuchos::RCP<OutputManager<ScalarType> > printer_;
315 Teuchos::RCP<std::ostream> outputStream_;
323 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
328 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
333 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
338 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > impConvTest_;
343 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > expConvTest_;
349 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
354 mutable Teuchos::RCP<const Teuchos::ParameterList> defaultParams_;
357 Teuchos::RCP<Teuchos::ParameterList> params_;
360 MagnitudeType convtol_;
363 MagnitudeType achievedTol_;
387 Teuchos::RCP<Teuchos::Time> timerSolve_;
402 template<
class ScalarType,
class MV,
class OP>
403 Teuchos::RCP<const Teuchos::ParameterList>
406 using Teuchos::ParameterList;
407 using Teuchos::parameterList;
410 using Teuchos::rcpFromRef;
411 using Teuchos::EnhancedNumberValidator;
412 typedef MagnitudeType MT;
413 typedef Teuchos::ScalarTraits<MT> MST;
416 RCP<ParameterList> pl = parameterList (
"MINRES");
418 pl->set (
"Convergence Tolerance", MST::squareroot (MST::eps()),
419 "Relative residual tolerance that needs to be achieved by "
420 "the iterative solver, in order for the linear system to be "
421 "declared converged.",
422 rcp (
new EnhancedNumberValidator<MT> (MST::zero(), MST::rmax())));
423 pl->set (
"Maximum Iterations", static_cast<int>(1000),
424 "Maximum number of iterations allowed for each right-hand "
426 rcp (
new EnhancedNumberValidator<int> (0, INT_MAX)));
427 pl->set (
"Num Blocks", static_cast<int> (-1),
428 "Ignored, but permitted, for compatibility with other Belos "
430 pl->set (
"Block Size", static_cast<int> (1),
431 "Number of vectors in each block. WARNING: The current "
432 "implementation of MINRES only accepts a block size of 1, "
433 "since it can only solve for 1 right-hand side at a time.",
434 rcp (
new EnhancedNumberValidator<int> (1, 1)));
436 "The type(s) of solver information that should "
437 "be written to the output stream.");
439 "What style is used for the solver information written "
440 "to the output stream.");
441 pl->set (
"Output Frequency", static_cast<int>(-1),
442 "How often (in terms of number of iterations) intermediate "
443 "convergence information should be written to the output stream."
445 pl->set (
"Output Stream", rcpFromRef(std::cout),
446 "A reference-counted pointer to the output stream where all "
447 "solver output is sent. The output stream defaults to stdout.");
448 pl->set (
"Timer Label", std::string(
"Belos"),
449 "The string to use as a prefix for the timer labels.");
456 template<
class ScalarType,
class MV,
class OP>
466 parametersSet_ (false)
472 template<
class ScalarType,
class MV,
class OP>
475 const Teuchos::RCP<Teuchos::ParameterList>& params) :
478 parametersSet_ (false)
480 TEUCHOS_TEST_FOR_EXCEPTION(problem_.is_null(), std::invalid_argument,
481 "MinresSolMgr: The version of the constructor "
482 "that takes a LinearProblem to solve was given a "
483 "null LinearProblem.");
487 template<
class ScalarType,
class MV,
class OP>
492 TEUCHOS_TEST_FOR_EXCEPTION(problem.is_null(),
494 "MINRES requires that you have provided a nonnull LinearProblem to the "
495 "solver manager, before you call the solve() method.");
496 TEUCHOS_TEST_FOR_EXCEPTION(problem->getOperator().is_null(),
498 "MINRES requires a LinearProblem object with a non-null operator (the "
500 TEUCHOS_TEST_FOR_EXCEPTION(problem->getRHS().is_null(),
502 "MINRES requires a LinearProblem object with a non-null right-hand side.");
503 TEUCHOS_TEST_FOR_EXCEPTION( ! problem->isProblemSet(),
505 "MINRES requires that before you give it a LinearProblem to solve, you "
506 "must first call the linear problem's setProblem() method.");
509 template<
class ScalarType,
class MV,
class OP>
514 using Teuchos::ParameterList;
515 using Teuchos::parameterList;
518 using Teuchos::rcpFromRef;
520 using Teuchos::is_null;
525 if (params_.is_null()) {
526 params_ = parameterList (*getValidParameters());
528 RCP<ParameterList> pl = params;
529 pl->validateParametersAndSetDefaults (*params_);
535 blockSize_ = pl->get<
int> (
"Block Size");
536 verbosity_ = pl->get<
int> (
"Verbosity");
537 outputStyle_ = pl->get<
int> (
"Output Style");
538 outputFreq_ = pl->get<
int>(
"Output Frequency");
539 outputStream_ = pl->get<RCP<std::ostream> > (
"Output Stream");
540 convtol_ = pl->get<MagnitudeType> (
"Convergence Tolerance");
541 maxIters_ = pl->get<
int> (
"Maximum Iterations");
549 const string newLabel = pl->get<
string> (
"Timer Label");
551 if (newLabel != label_ || timerSolve_.is_null()) {
553 #ifdef BELOS_TEUCHOS_TIME_MONITOR
554 const string solveLabel = label_ +
": MinresSolMgr total solve time";
556 if (! timerSolve_.is_null()) {
557 Teuchos::TimeMonitor::clearCounter (label_);
558 timerSolve_ = Teuchos::null;
560 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
561 #endif // BELOS_TEUCHOS_TIME_MONITOR
566 bool recreatedPrinter =
false;
567 if (printer_.is_null()) {
569 recreatedPrinter =
true;
572 printer_->setVerbosity (verbosity_);
574 printer_->setOStream (outputStream_);
585 const bool allocatedConvergenceTests =
586 impConvTest_.is_null() || expConvTest_.is_null();
590 if (impConvTest_.is_null()) {
591 impConvTest_ = rcp (
new res_norm_type (convtol_));
592 impConvTest_->defineResForm (res_norm_type::Implicit,
TwoNorm);
597 impConvTest_->setTolerance (convtol_);
602 if (expConvTest_.is_null()) {
603 expConvTest_ = rcp (
new res_norm_type (convtol_));
604 expConvTest_->defineResForm (res_norm_type::Explicit,
TwoNorm);
609 expConvTest_->setTolerance (convtol_);
615 bool needToRecreateFullStatusTest = sTest_.is_null();
619 if (convTest_.is_null() || allocatedConvergenceTests) {
620 convTest_ = rcp (
new combo_type (combo_type::SEQ, impConvTest_, expConvTest_));
621 needToRecreateFullStatusTest =
true;
628 if (maxIterTest_.is_null()) {
630 needToRecreateFullStatusTest =
true;
632 maxIterTest_->setMaxIters (maxIters_);
643 if (needToRecreateFullStatusTest) {
644 sTest_ = rcp (
new combo_type (combo_type::OR, maxIterTest_, convTest_));
651 if (outputTest_.is_null() || needToRecreateFullStatusTest || recreatedPrinter) {
653 outputTest_ = stoFactory.
create (printer_, sTest_, outputFreq_,
656 outputTest_->setOutputFrequency (outputFreq_);
660 outputTest_->setSolverDesc (std::string (
" MINRES "));
663 parametersSet_ =
true;
665 if (verbosity_ &
Debug) {
668 std::ostream& dbg = printer_->stream (
Debug);
669 dbg <<
"MINRES parameters:" << endl << params_ << endl;
674 template<
class ScalarType,
class MV,
class OP>
679 using Teuchos::rcp_const_cast;
682 if (! parametersSet_) {
683 setParameters (params_);
685 std::ostream& dbg = printer_->stream (
Debug);
687 #ifdef BELOS_TEUCHOS_TIME_MONITOR
688 Teuchos::TimeMonitor solveTimerMonitor (*timerSolve_);
689 #endif // BELOS_TEUCHOS_TIME_MONITOR
692 validateProblem (problem_);
695 outputTest_->reset();
700 const int numRHS2Solve = MVT::GetNumberVecs (*(problem_->getRHS()));
705 RCP<iter_type> minres_iter =
706 rcp (
new iter_type (problem_, printer_, outputTest_, *params_));
712 std::vector<int> notConverged;
713 std::vector<int> currentIndices(1);
718 for (
int currentRHS = 0; currentRHS < numRHS2Solve; ++currentRHS) {
723 currentIndices[0] = currentRHS;
724 problem_->setLSIndex (currentIndices);
726 dbg <<
"-- Current right-hand side index being solved: "
727 << currentRHS << endl;
730 minres_iter->resetNumIters();
732 outputTest_->resetNumCalls();
738 newstate.
Y = MVT::CloneViewNonConst (*(rcp_const_cast<MV> (problem_->getInitResVec())), currentIndices);
739 minres_iter->initializeMinres (newstate);
745 minres_iter->iterate();
748 if (convTest_->getStatus() ==
Passed) {
749 dbg <<
"---- Converged after " << maxIterTest_->getNumIters()
750 <<
" iterations" << endl;
754 else if (maxIterTest_->getStatus() ==
Passed) {
755 dbg <<
"---- Did not converge after " << maxIterTest_->getNumIters()
756 <<
" iterations" << endl;
758 notConverged.push_back (currentRHS);
764 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
765 "Belos::MinresSolMgr::solve(): iterations neither converged, "
766 "nor reached the maximum number of iterations " << maxIters_
767 <<
". That means something went wrong.");
769 }
catch (
const std::exception &e) {
771 <<
"Error! Caught std::exception in MinresIter::iterate() at "
772 <<
"iteration " << minres_iter->getNumIters() << endl
781 problem_->setCurrLS();
785 numIters_ += maxIterTest_->getNumIters();
793 #ifdef BELOS_TEUCHOS_TIME_MONITOR
798 Teuchos::TimeMonitor::summarize (printer_->stream (
TimingDetails));
800 #endif // BELOS_TEUCHOS_TIME_MONITOR
812 const std::vector<MagnitudeType>* pTestValues = expConvTest_->getTestValue();
813 if (pTestValues == NULL || pTestValues->size() < 1) {
814 pTestValues = impConvTest_->getTestValue();
816 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
817 "Belos::MinresSolMgr::solve(): The implicit convergence test's getTestValue() "
818 "method returned NULL. Please report this bug to the Belos developers.");
819 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
820 "Belos::MinresSolMgr::solve(): The implicit convergence test's getTestValue() "
821 "method returned a vector of length zero. Please report this bug to the "
822 "Belos developers.");
827 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
830 if (notConverged.size() > 0) {
838 template<
class ScalarType,
class MV,
class OP>
841 std::ostringstream oss;
842 oss <<
"Belos::MinresSolMgr< "
843 << Teuchos::ScalarTraits<ScalarType>::name()