Belos Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosLSQRSolMgr.hpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
4 // Belos: Block Linear Solvers Package
5 // Copyright 2004 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
42 #ifndef BELOS_LSQR_SOLMGR_HPP
43 #define BELOS_LSQR_SOLMGR_HPP
44 
47 
48 #include "BelosConfigDefs.hpp"
49 #include "BelosTypes.hpp"
50 
51 #include "BelosLinearProblem.hpp"
52 #include "BelosSolverManager.hpp"
53 
54 #include "BelosLSQRIteration.hpp"
55 #include "BelosLSQRIter.hpp"
57 #include "BelosLSQRStatusTest.hpp"
58 #include "BelosStatusTestCombo.hpp"
60 #include "BelosOutputManager.hpp"
61 #include "Teuchos_as.hpp"
62 
63 #ifdef BELOS_TEUCHOS_TIME_MONITOR
64 #include "Teuchos_TimeMonitor.hpp"
65 #endif
66 
67 namespace Belos {
68 
69 
71 
72 
80 public:
81  LSQRSolMgrLinearProblemFailure(const std::string& what_arg)
82  : BelosError(what_arg)
83  {}
84 };
85 
93 public:
94  LSQRSolMgrBlockSizeFailure(const std::string& what_arg)
95  : BelosError(what_arg)
96  {}
97 };
98 
212 
213 
214 // Partial specialization for complex ScalarType.
215 // This contains a trivial implementation.
216 // See discussion in the class documentation above.
217 template<class ScalarType, class MV, class OP,
218  const bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex>
219 class LSQRSolMgr :
220  public Details::RealSolverManager<ScalarType, MV, OP,
221  Teuchos::ScalarTraits<ScalarType>::isComplex>
222 {
225 
226 public:
228  base_type ()
229  {}
232  base_type ()
233  {}
234  virtual ~LSQRSolMgr () {}
235 
239  }
240 };
241 
242 
243 // Partial specialization for real ScalarType.
244 // This contains the actual working implementation of LSQR.
245 // See discussion in the class documentation above.
246 template<class ScalarType, class MV, class OP>
247 class LSQRSolMgr<ScalarType, MV, OP, false> :
248  public Details::RealSolverManager<ScalarType, MV, OP, false> {
249 private:
255 
256 public:
257 
259 
260 
267  LSQRSolMgr ();
268 
298 
300  virtual ~LSQRSolMgr () {}
301 
305  }
307 
309 
312  const LinearProblem<ScalarType,MV,OP>& getProblem () const override {
313  return *problem_;
314  }
315 
318  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const override;
319 
323  return params_;
324  }
325 
332  return Teuchos::tuple (timerSolve_);
333  }
334 
336  int getNumIters () const override {
337  return numIters_;
338  }
339 
345  return matCondNum_;
346  }
347 
353  return matNorm_;
354  }
355 
365  return resNorm_;
366  }
367 
370  return matResNorm_;
371  }
372 
381  bool isLOADetected () const override { return false; }
382 
384 
386 
387 
389  void setProblem (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >& problem) override {
390  problem_ = problem;
391  }
392 
394  void setParameters (const Teuchos::RCP<Teuchos::ParameterList>& params) override;
395 
397 
399 
400 
404  void reset (const ResetType type) override {
405  if ((type & Belos::Problem) && ! problem_.is_null ()) {
406  problem_->setProblem ();
407  }
408  }
409 
411 
413 
432  ReturnType solve() override;
433 
435 
437 
439  std::string description () const override;
440 
442 
443 private:
444 
451 
457 
460 
467 
468  // Current solver input parameters
473  int maxIters_, termIterMax_;
474  int verbosity_, outputStyle_, outputFreq_;
475 
476  // Terminal solver state values
482 
483  // Timers.
484  std::string label_;
486 
487  // Internal state variables.
488  bool isSet_;
490 };
491 
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 ()),
498  maxIters_ (1000),
499  termIterMax_ (1),
500  verbosity_ (Belos::Errors),
501  outputStyle_ (Belos::General),
502  outputFreq_ (-1),
503  numIters_ (0),
504  matCondNum_ (STM::zero ()),
505  matNorm_ (STM::zero ()),
506  resNorm_ (STM::zero ()),
507  matResNorm_ (STM::zero ()),
508  isSet_ (false),
509  loaDetected_ (false)
510 {}
511 
512 template<class ScalarType, class MV, class OP>
516  problem_ (problem),
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 ()),
521  maxIters_ (1000),
522  termIterMax_ (1),
523  verbosity_ (Belos::Errors),
524  outputStyle_ (Belos::General),
525  outputFreq_ (-1),
526  numIters_ (0),
527  matCondNum_ (STM::zero ()),
528  matNorm_ (STM::zero ()),
529  resNorm_ (STM::zero ()),
530  matResNorm_ (STM::zero ()),
531  isSet_ (false),
532  loaDetected_ (false)
533 {
534  // The linear problem to solve is allowed to be null here. The user
535  // must then set a nonnull linear problem (by calling setProblem())
536  // before calling solve().
537  //
538  // Similarly, users are allowed to set a null parameter list here,
539  // but they must first set a nonnull parameter list (by calling
540  // setParameters()) before calling solve().
541  if (! pl.is_null ()) {
542  setParameters (pl);
543  }
544 }
545 
546 
547 template<class ScalarType, class MV, class OP>
550 {
552  using Teuchos::parameterList;
553  using Teuchos::RCP;
554  using Teuchos::rcp;
555  using Teuchos::rcpFromRef;
556 
557  // Set all the valid parameters and their default values.
558  if (validParams_.is_null ()) {
559  // We use Teuchos::as just in case MagnitudeType doesn't have a
560  // constructor that takes an int. Otherwise, we could just write
561  // "MagnitudeType(10)".
562  const MagnitudeType ten = Teuchos::as<MagnitudeType> (10);
563  const MagnitudeType sqrtEps = STM::squareroot (STM::eps());
564 
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;
572  const int verbosity = Belos::Errors;
573  const int outputStyle = Belos::General;
574  const int outputFreq = -1;
575  const std::string label ("Belos");
576 
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 "
587  "number of Abar.");
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 "
598  "labels");
599  pl->set ("Block Size", 1, "Block size parameter (currently, LSQR requires "
600  "this must always be 1)");
601  validParams_ = pl;
602  }
603  return validParams_;
604 }
605 
606 
607 template<class ScalarType, class MV, class OP>
608 void
611 {
612  using Teuchos::isParameterType;
613  using Teuchos::getParameter;
614  using Teuchos::null;
616  using Teuchos::parameterList;
617  using Teuchos::RCP;
618  using Teuchos::rcp;
619  using Teuchos::rcp_dynamic_cast;
620  using Teuchos::rcpFromRef;
621  using Teuchos::Time;
622  using Teuchos::TimeMonitor;
626 
628  (params.is_null (), std::invalid_argument,
629  "Belos::LSQRSolMgr::setParameters: The input ParameterList is null.");
630  RCP<const ParameterList> defaultParams = getValidParameters ();
631 
632  // FIXME (mfh 29 Apr 2015) Our users would like to supply one
633  // ParameterList that works for both GMRES and LSQR. Thus, we want
634  // LSQR (the less-used solver) to ignore parameters it doesn't
635  // recognize). For now, therefore, it should not validate, since
636  // validation cannot distinguish between misspellings and
637  // unrecognized parameters. (Perhaps Belos should have a central
638  // facility for all parameters recognized by some solver in Belos,
639  // so we could use that for spell checking.)
640  //
641  //params->validateParameters (*defaultParams);
642 
643  // mfh 29 Apr 2015: The convention in Belos is that the input
644  // ParameterList is a "delta" from the current state. Thus, we
645  // don't fill in the input ParameterList with defaults, and we only
646  // change the current state if the corresponding parameter was
647  // explicitly set in the input ParameterList. We set up the solver
648  // with the default state on construction.
649 
650  // Get the damping (regularization) parameter lambda.
651  if (params->isParameter ("Lambda")) {
652  lambda_ = params->get<MagnitudeType> ("Lambda");
653  } else if (params->isParameter ("lambda")) {
654  lambda_ = params->get<MagnitudeType> ("lambda");
655  }
656 
657  // Get the maximum number of iterations.
658  if (params->isParameter ("Maximum Iterations")) {
659  maxIters_ = params->get<int> ("Maximum Iterations");
660  }
662  (maxIters_ < 0, std::invalid_argument, "Belos::LSQRSolMgr::setParameters: "
663  "\"Maximum Iterations\" = " << maxIters_ << " < 0.");
664 
665  // (Re)set the timer label.
666  {
667  const std::string newLabel =
668  params->isParameter ("Timer Label") ?
669  params->get<std::string> ("Timer Label") :
670  label_;
671 
672  // Update parameter in our list and solver timer
673  if (newLabel != label_) {
674  label_ = newLabel;
675  }
676 
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 ()) {
682  // Ask TimeMonitor for a new timer.
683  timerSolve_ = TimeMonitor::getNewCounter (newSolveLabel);
684  } else {
685  // We've already created a timer, but we may have changed its
686  // label. If we did change its name, then we have to forget
687  // about the old timer and create a new one with a different
688  // name. This is because Teuchos::Time doesn't give you a way
689  // to change a timer's name, once you've created it. We assume
690  // that if the user changed the timer's label, then the user
691  // wants to reset the timer's results.
692  const std::string oldSolveLabel = timerSolve_->name ();
693 
694  if (oldSolveLabel != newSolveLabel) {
695  // Tell TimeMonitor to forget about the old timer.
696  // TimeMonitor lets you clear timers by name.
697  TimeMonitor::clearCounter (oldSolveLabel);
698  timerSolve_ = TimeMonitor::getNewCounter (newSolveLabel);
699  }
700  }
701 #endif // BELOS_TEUCHOS_TIME_MONITOR
702  }
703 
704  // Check for a change in verbosity level
705  if (params->isParameter ("Verbosity")) {
706  int newVerbosity = 0;
707  // ParameterList gets confused sometimes about enums. This
708  // ensures that no matter how "Verbosity" was stored -- either an
709  // an int, or as a Belos::MsgType enum, we will be able to extract
710  // it. If it was stored as some other type, we let the exception
711  // through.
712  try {
713  newVerbosity = params->get<Belos::MsgType> ("Verbosity");
715  newVerbosity = params->get<int> ("Verbosity");
716  }
717  if (newVerbosity != verbosity_) {
718  verbosity_ = newVerbosity;
719  }
720  }
721 
722  // (Re)set the output style.
723  if (params->isParameter ("Output Style")) {
724  outputStyle_ = params->get<int> ("Output Style");
725  }
726 
727  // Get the output stream for the output manager.
728  //
729  // In case the output stream can't be read back in, we default to
730  // stdout (std::cout), just to ensure reasonable behavior.
731  if (params->isParameter ("Output Stream")) {
732  outputStream_ = params->get<RCP<std::ostream> > ("Output Stream");
733  }
734  // We assume that a null output stream indicates that the user
735  // doesn't want to print anything, so we replace it with a "black
736  // hole" stream that prints nothing sent to it. (We can't use a
737  // null output stream, since the output manager always sends
738  // things it wants to print to the output stream.)
739  if (outputStream_.is_null ()) {
740  outputStream_ = rcp (new Teuchos::oblackholestream ());
741  }
742 
743  // Get the frequency of solver output. (For example, -1 means
744  // "never," and 1 means "every iteration.")
745  if (params->isParameter ("Output Frequency")) {
746  outputFreq_ = params->get<int> ("Output Frequency");
747  }
748 
749  // Create output manager if we need to, using the verbosity level
750  // and output stream that we fetched above. Status tests (i.e.,
751  // stopping criteria) need this.
752  if (printer_.is_null ()) {
753  printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
754  } else {
755  printer_->setVerbosity (verbosity_);
756  printer_->setOStream (outputStream_);
757  }
758 
759  // Check for condition number limit, number of consecutive passed
760  // iterations, relative RHS error, and relative matrix error.
761  // Create the LSQR convergence test if necessary.
762  {
763  if (params->isParameter ("Condition Limit")) {
764  condMax_ = params->get<MagnitudeType> ("Condition Limit");
765  }
766  if (params->isParameter ("Term Iter Max")) {
767  termIterMax_ = params->get<int> ("Term Iter Max");
768  }
769  if (params->isParameter ("Rel RHS Err")) {
770  relRhsErr_ = params->get<MagnitudeType> ("Rel RHS Err");
771  }
772  else if (params->isParameter ("Convergence Tolerance")) {
773  // NOTE (mfh 29 Apr 2015) We accept this parameter as an alias
774  // for "Rel RHS Err".
775  relRhsErr_ = params->get<MagnitudeType> ("Convergence Tolerance");
776  }
777 
778  if (params->isParameter ("Rel Mat Err")) {
779  relMatErr_ = params->get<MagnitudeType> ("Rel Mat Err");
780  }
781 
782  // Create the LSQR convergence test if it doesn't exist yet.
783  // Otherwise, update its parameters.
784  if (convTest_.is_null ()) {
785  convTest_ =
786  rcp (new LSQRStatusTest<ScalarType,MV,OP> (condMax_, termIterMax_,
787  relRhsErr_, relMatErr_));
788  } else {
789  convTest_->setCondLim (condMax_);
790  convTest_->setTermIterMax (termIterMax_);
791  convTest_->setRelRhsErr (relRhsErr_);
792  convTest_->setRelMatErr (relMatErr_);
793  }
794  }
795 
796  // Create the status test for maximum number of iterations if
797  // necessary. Otherwise, update it with the new maximum iteration
798  // count.
799  if (maxIterTest_.is_null()) {
800  maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
801  } else {
802  maxIterTest_->setMaxIters (maxIters_);
803  }
804 
805  // The stopping criterion is an OR combination of the test for
806  // maximum number of iterations, and the LSQR convergence test.
807  // ("OR combination" means that both tests will always be evaluated,
808  // as opposed to a SEQ combination.)
809  typedef StatusTestCombo<ScalarType,MV,OP> combo_type;
810  // If sTest_ is not null, then maxIterTest_ and convTest_ were
811  // already constructed on entry to this routine, and sTest_ has
812  // their pointers. Thus, maxIterTest_ and convTest_ have gotten any
813  // parameter changes, so we don't need to do anything to sTest_.
814  if (sTest_.is_null()) {
815  sTest_ = rcp (new combo_type (combo_type::OR, maxIterTest_, convTest_));
816  }
817 
818  if (outputTest_.is_null ()) {
819  // Create the status test output class.
820  // This class manages and formats the output from the status test.
821  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
822  outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_,
823  Passed + Failed + Undefined);
824  // Set the solver string for the output test.
825  const std::string solverDesc = " LSQR ";
826  outputTest_->setSolverDesc (solverDesc);
827  } else {
828  // FIXME (mfh 18 Sep 2011) How do we change the output style of
829  // outputTest_, without destroying and recreating it?
830  outputTest_->setOutputManager (printer_);
831  outputTest_->setChild (sTest_);
832  outputTest_->setOutputFrequency (outputFreq_);
833  // Since outputTest_ can only be created here, I'm assuming that
834  // the fourth constructor argument ("printStates") was set
835  // correctly on constrution; I don't need to reset it (and I can't
836  // set it anyway, given StatusTestOutput's interface).
837  }
838 
839  // At this point, params is a valid ParameterList. Now we can
840  // "commit" it to our instance's ParameterList.
841  params_ = params;
842 
843  // Inform the solver manager that the current parameters were set.
844  isSet_ = true;
845 }
846 
847 
848 template<class ScalarType, class MV, class OP>
851 {
852  using Teuchos::RCP;
853  using Teuchos::rcp;
854 
855  // Set the current parameters if they were not set before. NOTE:
856  // This may occur if the user generated the solver manager with the
857  // default constructor, but did not set any parameters using
858  // setParameters().
859  if (! isSet_) {
860  this->setParameters (Teuchos::parameterList (* (getValidParameters ())));
861  }
862 
864  (problem_.is_null (), LSQRSolMgrLinearProblemFailure,
865  "Belos::LSQRSolMgr::solve: The linear problem to solve is null.");
867  (! problem_->isProblemSet (), LSQRSolMgrLinearProblemFailure,
868  "Belos::LSQRSolMgr::solve: The linear problem is not ready, "
869  "as its setProblem() method has not been called.");
871  (MVT::GetNumberVecs (*(problem_->getRHS ())) != 1,
872  LSQRSolMgrBlockSizeFailure, "Belos::LSQRSolMgr::solve: "
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.");
877 
878  // We've validated the LinearProblem instance above. If any of the
879  // StatusTests needed to be initialized using information from the
880  // LinearProblem, now would be the time to do so. (This is true of
881  // GMRES, where the residual convergence test(s) to instantiate
882  // depend on knowing whether there is a left preconditioner. This
883  // is why GMRES has an "isSTSet_" Boolean member datum, which tells
884  // you whether the status tests have been instantiated and are ready
885  // for use.
886 
887  // test isFlexible might go here.
888 
889  // Next the right-hand sides to solve are identified. Among other things,
890  // this enables getCurrLHSVec() to get the current initial guess vector,
891  // and getCurrRHSVec() to get the current right-hand side (in Iter).
892  std::vector<int> currRHSIdx (1, 0);
893  problem_->setLSIndex (currRHSIdx);
894 
895  // Reset the status test.
896  outputTest_->reset ();
897 
898  // Don't assume convergence unless we've verified that the
899  // convergence test passed.
900  bool isConverged = false;
901 
902  // FIXME: Currently we are setting the initial guess to zero, since
903  // the solver doesn't yet know how to handle a nonzero initial
904  // guess. This could be fixed by rewriting the solver to work with
905  // the residual and a delta.
906  //
907  // In a least squares problem with a nonzero initial guess, the
908  // minimzation problem involves the distance (in a norm depending on
909  // the preconditioner) between the solution and the the initial
910  // guess.
911 
913  // Solve the linear problem using LSQR
915 
916  // Parameter list for the LSQR iteration.
918 
919  // Use the solver manager's "Lambda" parameter to set the
920  // iteration's "Lambda" parameter. We know that the solver
921  // manager's parameter list (params_) does indeed contain the
922  // "Lambda" parameter, because solve() always ensures that
923  // setParameters() has been called.
924  plist.set ("Lambda", lambda_);
925 
926  typedef LSQRIter<ScalarType,MV,OP> iter_type;
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_);
931 #endif
932 
933  // Reset the number of iterations.
934  lsqr_iter->resetNumIters ();
935  // Reset the number of calls that the status test output knows about.
936  outputTest_->resetNumCalls ();
937  // Set the new state and initialize the solver.
939  lsqr_iter->initializeLSQR (newstate);
940  // tell lsqr_iter to iterate
941  try {
942  lsqr_iter->iterate ();
943 
944  // First check for convergence. If we didn't converge, then check
945  // whether we reached the maximum number of iterations. If
946  // neither of those happened, there must have been a bug.
947  if (convTest_->getStatus () == Belos::Passed) {
948  isConverged = true;
949  } else if (maxIterTest_->getStatus () == Belos::Passed) {
950  isConverged = false;
951  } else {
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.");
957  }
958  } catch (const std::exception& e) {
959  printer_->stream(Belos::Errors)
960  << "Error! Caught std::exception in LSQRIter::iterate at iteration "
961  << lsqr_iter->getNumIters () << std::endl << e.what () << std::endl;
962  throw;
963  }
964 
965  // identify current linear system as solved LinearProblem
966  problem_->setCurrLS();
967  // print final summary
968  sTest_->print (printer_->stream (Belos::FinalSummary));
969 
970  // Print timing information, if the corresponding compile-time and
971  // run-time options are enabled.
972 #ifdef BELOS_TEUCHOS_TIME_MONITOR
973  // Calling summarize() can be expensive, so don't call unless the
974  // user wants to print out timing details. summarize() will do all
975  // the work even if it's passed a "black hole" output stream.
976  if (verbosity_ & TimingDetails)
978 #endif // BELOS_TEUCHOS_TIME_MONITOR
979 
980  // A posteriori solve information
981  numIters_ = maxIterTest_->getNumIters();
982  matCondNum_ = convTest_->getMatCondNum();
983  matNorm_ = convTest_->getMatNorm();
984  resNorm_ = convTest_->getResidNorm();
985  matResNorm_ = convTest_->getLSResidNorm();
986 
987  if (! isConverged) {
988  return Belos::Unconverged;
989  } else {
990  return Belos::Converged;
991  }
992 }
993 
994 // LSQRSolMgr requires the solver manager to return an eponymous std::string.
995 template<class ScalarType, class MV, class OP>
997 {
998  std::ostringstream oss;
999  oss << "LSQRSolMgr<...," << STS::name () << ">";
1000  oss << "{";
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_;
1007  oss << "}";
1008  return oss.str ();
1009 }
1010 
1011 } // end Belos namespace
1012 
1013 #endif /* BELOS_LSQR_SOLMGR_HPP */
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...
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
Belos concrete class that iterates LSQR.
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
Class which manages the output and verbosity of the Belos solvers.
LSQRSolMgrBlockSizeFailure(const std::string &what_arg)
MagnitudeType getMatCondNum() const
Estimated matrix condition number from the last solve.
OperatorTraits< ScalarType, MV, OP > OPT
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
Teuchos::ScalarTraits< MagnitudeType > STM
Details::RealSolverManager< ScalarType, MV, OP, isComplex > base_type
LSQRSolMgrBlockSizeFailure is thrown when the linear problem has more than one RHS.
T & get(ParameterList &l, const std::string &name)
MsgType
Available message types recognized by the linear solvers.
Definition: BelosTypes.hpp:254
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
A factory class for generating StatusTestOutput objects.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
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.
MagnitudeType getMatNorm() const
Estimated matrix Frobenius norm 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.
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
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...
Teuchos::RCP< Teuchos::ParameterList > params_
Current parameter list.
bool isParameter(const std::string &name) const
Teuchos::ScalarTraits< ScalarType > STS
static const bool isComplex
Teuchos::RCP< OutputManager< ScalarType > > printer_
The output manager.
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.
Definition: BelosTypes.hpp:206
Teuchos::RCP< LSQRStatusTest< ScalarType, MV, OP > > convTest_
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Pure virtual base class which describes the basic interface for a solver manager. ...
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
The &quot;master&quot; status test (that includes all status tests).
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
bool is_null(const RCP< T > &p)
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.
MagnitudeType getResNorm() const
Estimated residual norm from the last solve.
Teuchos::RCP< const Teuchos::ParameterList > validParams_
Default parameter list.
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:155
TypeTo as(const TypeFrom &t)
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
Teuchos::RCP< std::ostream > outputStream_
Output stream to which to write status output.
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.
MagnitudeType getMatResNorm() const
Estimate of (residual vector ) from the last solve.
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.
Definition: BelosTypes.hpp:60
Structure to contain pointers to LSQRIteration state variables, ...
Belos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
The linear problem to solve.
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Get current linear problem being solved for in this object.
LSQRSolMgrLinearProblemFailure(const std::string &what_arg)
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_