Belos  Version of the Day
 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 #include "Teuchos_BLAS.hpp"
63 #include "Teuchos_LAPACK.hpp"
64 
65 #ifdef BELOS_TEUCHOS_TIME_MONITOR
66 #include "Teuchos_TimeMonitor.hpp"
67 #endif
68 
69 namespace Belos {
70 
71 
73 
74 
82 public:
83  LSQRSolMgrLinearProblemFailure(const std::string& what_arg)
84  : BelosError(what_arg)
85  {}
86 };
87 
97 public:
98  LSQRSolMgrOrthoFailure(const std::string& what_arg)
99  : BelosError(what_arg)
100  {}
101 };
102 
110 public:
111  LSQRSolMgrBlockSizeFailure(const std::string& what_arg)
112  : BelosError(what_arg)
113  {}
114 };
115 
229 
230 
231 // Partial specialization for complex ScalarType.
232 // This contains a trivial implementation.
233 // See discussion in the class documentation above.
234 template<class ScalarType, class MV, class OP,
235  const bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex>
236 class LSQRSolMgr :
237  public Details::RealSolverManager<ScalarType, MV, OP,
238  Teuchos::ScalarTraits<ScalarType>::isComplex>
239 {
240  static const bool isComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
242 
243 public:
245  base_type ()
246  {}
249  base_type ()
250  {}
251  virtual ~LSQRSolMgr () {}
252 
256  }
257 };
258 
259 
260 // Partial specialization for real ScalarType.
261 // This contains the actual working implementation of LSQR.
262 // See discussion in the class documentation above.
263 template<class ScalarType, class MV, class OP>
264 class LSQRSolMgr<ScalarType, MV, OP, false> :
265  public Details::RealSolverManager<ScalarType, MV, OP, false> {
266 private:
270  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
272 
273 public:
274 
276 
277 
284  LSQRSolMgr ();
285 
315 
317  virtual ~LSQRSolMgr () {}
318 
322  }
324 
326 
329  const LinearProblem<ScalarType,MV,OP>& getProblem () const override {
330  return *problem_;
331  }
332 
335  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const override;
336 
340  return params_;
341  }
342 
349  return Teuchos::tuple (timerSolve_);
350  }
351 
353  int getNumIters () const override {
354  return numIters_;
355  }
356 
361  MagnitudeType getMatCondNum () const {
362  return matCondNum_;
363  }
364 
369  MagnitudeType getMatNorm () const {
370  return matNorm_;
371  }
372 
381  MagnitudeType getResNorm () const {
382  return resNorm_;
383  }
384 
386  MagnitudeType getMatResNorm () const {
387  return matResNorm_;
388  }
389 
398  bool isLOADetected () const override { return false; }
399 
401 
403 
404 
406  void setProblem (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >& problem) override {
407  problem_ = problem;
408  }
409 
411  void setParameters (const Teuchos::RCP<Teuchos::ParameterList>& params) override;
412 
414 
416 
417 
421  void reset (const ResetType type) override {
422  if ((type & Belos::Problem) && ! problem_.is_null ()) {
423  problem_->setProblem ();
424  }
425  }
426 
428 
430 
449  ReturnType solve() override;
450 
452 
454 
456  std::string description () const override;
457 
459 
460 private:
461 
467  Teuchos::RCP<std::ostream> outputStream_;
468 
474 
477 
483  mutable Teuchos::RCP<const Teuchos::ParameterList> validParams_;
484 
485  // Current solver input parameters
486  MagnitudeType lambda_;
487  MagnitudeType relRhsErr_;
488  MagnitudeType relMatErr_;
489  MagnitudeType condMax_;
490  int maxIters_, termIterMax_;
491  int verbosity_, outputStyle_, outputFreq_;
492 
493  // Terminal solver state values
494  int numIters_;
495  MagnitudeType matCondNum_;
496  MagnitudeType matNorm_;
497  MagnitudeType resNorm_;
498  MagnitudeType matResNorm_;
499 
500  // Timers.
501  std::string label_;
502  Teuchos::RCP<Teuchos::Time> timerSolve_;
503 
504  // Internal state variables.
505  bool isSet_;
506  bool loaDetected_;
507 };
508 
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 ()),
515  maxIters_ (1000),
516  termIterMax_ (1),
517  verbosity_ (Belos::Errors),
518  outputStyle_ (Belos::General),
519  outputFreq_ (-1),
520  numIters_ (0),
521  matCondNum_ (STM::zero ()),
522  matNorm_ (STM::zero ()),
523  resNorm_ (STM::zero ()),
524  matResNorm_ (STM::zero ()),
525  isSet_ (false),
526  loaDetected_ (false)
527 {}
528 
529 template<class ScalarType, class MV, class OP>
533  problem_ (problem),
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 ()),
538  maxIters_ (1000),
539  termIterMax_ (1),
540  verbosity_ (Belos::Errors),
541  outputStyle_ (Belos::General),
542  outputFreq_ (-1),
543  numIters_ (0),
544  matCondNum_ (STM::zero ()),
545  matNorm_ (STM::zero ()),
546  resNorm_ (STM::zero ()),
547  matResNorm_ (STM::zero ()),
548  isSet_ (false),
549  loaDetected_ (false)
550 {
551  // The linear problem to solve is allowed to be null here. The user
552  // must then set a nonnull linear problem (by calling setProblem())
553  // before calling solve().
554  //
555  // Similarly, users are allowed to set a null parameter list here,
556  // but they must first set a nonnull parameter list (by calling
557  // setParameters()) before calling solve().
558  if (! pl.is_null ()) {
559  setParameters (pl);
560  }
561 }
562 
563 
564 template<class ScalarType, class MV, class OP>
567 {
569  using Teuchos::parameterList;
570  using Teuchos::RCP;
571  using Teuchos::rcp;
572  using Teuchos::rcpFromRef;
573 
574  // Set all the valid parameters and their default values.
575  if (validParams_.is_null ()) {
576  // We use Teuchos::as just in case MagnitudeType doesn't have a
577  // constructor that takes an int. Otherwise, we could just write
578  // "MagnitudeType(10)".
579  const MagnitudeType ten = Teuchos::as<MagnitudeType> (10);
580  const MagnitudeType sqrtEps = STM::squareroot (STM::eps());
581 
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;
589  const int verbosity = Belos::Errors;
590  const int outputStyle = Belos::General;
591  const int outputFreq = -1;
592  const std::string label ("Belos");
593 
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 "
604  "number of Abar.");
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 "
615  "labels");
616  pl->set ("Block Size", 1, "Block size parameter (currently, LSQR requires "
617  "this must always be 1)");
618  validParams_ = pl;
619  }
620  return validParams_;
621 }
622 
623 
624 template<class ScalarType, class MV, class OP>
625 void
628 {
629  using Teuchos::isParameterType;
630  using Teuchos::getParameter;
631  using Teuchos::null;
633  using Teuchos::parameterList;
634  using Teuchos::RCP;
635  using Teuchos::rcp;
636  using Teuchos::rcp_dynamic_cast;
637  using Teuchos::rcpFromRef;
638  using Teuchos::Time;
639  using Teuchos::TimeMonitor;
643 
645  (params.is_null (), std::invalid_argument,
646  "Belos::LSQRSolMgr::setParameters: The input ParameterList is null.");
647  RCP<const ParameterList> defaultParams = getValidParameters ();
648 
649  // FIXME (mfh 29 Apr 2015) Our users would like to supply one
650  // ParameterList that works for both GMRES and LSQR. Thus, we want
651  // LSQR (the less-used solver) to ignore parameters it doesn't
652  // recognize). For now, therefore, it should not validate, since
653  // validation cannot distinguish between misspellings and
654  // unrecognized parameters. (Perhaps Belos should have a central
655  // facility for all parameters recognized by some solver in Belos,
656  // so we could use that for spell checking.)
657  //
658  //params->validateParameters (*defaultParams);
659 
660  // mfh 29 Apr 2015: The convention in Belos is that the input
661  // ParameterList is a "delta" from the current state. Thus, we
662  // don't fill in the input ParameterList with defaults, and we only
663  // change the current state if the corresponding parameter was
664  // explicitly set in the input ParameterList. We set up the solver
665  // with the default state on construction.
666 
667  // Get the damping (regularization) parameter lambda.
668  if (params->isParameter ("Lambda")) {
669  lambda_ = params->get<MagnitudeType> ("Lambda");
670  } else if (params->isParameter ("lambda")) {
671  lambda_ = params->get<MagnitudeType> ("lambda");
672  }
673 
674  // Get the maximum number of iterations.
675  if (params->isParameter ("Maximum Iterations")) {
676  maxIters_ = params->get<int> ("Maximum Iterations");
677  }
679  (maxIters_ < 0, std::invalid_argument, "Belos::LSQRSolMgr::setParameters: "
680  "\"Maximum Iterations\" = " << maxIters_ << " < 0.");
681 
682  // (Re)set the timer label.
683  {
684  const std::string newLabel =
685  params->isParameter ("Timer Label") ?
686  params->get<std::string> ("Timer Label") :
687  label_;
688 
689  // Update parameter in our list and solver timer
690  if (newLabel != label_) {
691  label_ = newLabel;
692  }
693 
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 ()) {
699  // Ask TimeMonitor for a new timer.
700  timerSolve_ = TimeMonitor::getNewCounter (newSolveLabel);
701  } else {
702  // We've already created a timer, but we may have changed its
703  // label. If we did change its name, then we have to forget
704  // about the old timer and create a new one with a different
705  // name. This is because Teuchos::Time doesn't give you a way
706  // to change a timer's name, once you've created it. We assume
707  // that if the user changed the timer's label, then the user
708  // wants to reset the timer's results.
709  const std::string oldSolveLabel = timerSolve_->name ();
710 
711  if (oldSolveLabel != newSolveLabel) {
712  // Tell TimeMonitor to forget about the old timer.
713  // TimeMonitor lets you clear timers by name.
714  TimeMonitor::clearCounter (oldSolveLabel);
715  timerSolve_ = TimeMonitor::getNewCounter (newSolveLabel);
716  }
717  }
718 #endif // BELOS_TEUCHOS_TIME_MONITOR
719  }
720 
721  // Check for a change in verbosity level
722  if (params->isParameter ("Verbosity")) {
723  int newVerbosity = 0;
724  // ParameterList gets confused sometimes about enums. This
725  // ensures that no matter how "Verbosity" was stored -- either an
726  // an int, or as a Belos::MsgType enum, we will be able to extract
727  // it. If it was stored as some other type, we let the exception
728  // through.
729  try {
730  newVerbosity = params->get<Belos::MsgType> ("Verbosity");
732  newVerbosity = params->get<int> ("Verbosity");
733  }
734  if (newVerbosity != verbosity_) {
735  verbosity_ = newVerbosity;
736  }
737  }
738 
739  // (Re)set the output style.
740  if (params->isParameter ("Output Style")) {
741  outputStyle_ = params->get<int> ("Output Style");
742  }
743 
744  // Get the output stream for the output manager.
745  //
746  // In case the output stream can't be read back in, we default to
747  // stdout (std::cout), just to ensure reasonable behavior.
748  if (params->isParameter ("Output Stream")) {
749  outputStream_ = params->get<RCP<std::ostream> > ("Output Stream");
750  }
751  // We assume that a null output stream indicates that the user
752  // doesn't want to print anything, so we replace it with a "black
753  // hole" stream that prints nothing sent to it. (We can't use a
754  // null output stream, since the output manager always sends
755  // things it wants to print to the output stream.)
756  if (outputStream_.is_null ()) {
757  outputStream_ = rcp (new Teuchos::oblackholestream ());
758  }
759 
760  // Get the frequency of solver output. (For example, -1 means
761  // "never," and 1 means "every iteration.")
762  if (params->isParameter ("Output Frequency")) {
763  outputFreq_ = params->get<int> ("Output Frequency");
764  }
765 
766  // Create output manager if we need to, using the verbosity level
767  // and output stream that we fetched above. Status tests (i.e.,
768  // stopping criteria) need this.
769  if (printer_.is_null ()) {
770  printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
771  } else {
772  printer_->setVerbosity (verbosity_);
773  printer_->setOStream (outputStream_);
774  }
775 
776  // Check for condition number limit, number of consecutive passed
777  // iterations, relative RHS error, and relative matrix error.
778  // Create the LSQR convergence test if necessary.
779  {
780  if (params->isParameter ("Condition Limit")) {
781  condMax_ = params->get<MagnitudeType> ("Condition Limit");
782  }
783  if (params->isParameter ("Term Iter Max")) {
784  termIterMax_ = params->get<int> ("Term Iter Max");
785  }
786  if (params->isParameter ("Rel RHS Err")) {
787  relRhsErr_ = params->get<MagnitudeType> ("Rel RHS Err");
788  }
789  else if (params->isParameter ("Convergence Tolerance")) {
790  // NOTE (mfh 29 Apr 2015) We accept this parameter as an alias
791  // for "Rel RHS Err".
792  relRhsErr_ = params->get<MagnitudeType> ("Convergence Tolerance");
793  }
794 
795  if (params->isParameter ("Rel Mat Err")) {
796  relMatErr_ = params->get<MagnitudeType> ("Rel Mat Err");
797  }
798 
799  // Create the LSQR convergence test if it doesn't exist yet.
800  // Otherwise, update its parameters.
801  if (convTest_.is_null ()) {
802  convTest_ =
803  rcp (new LSQRStatusTest<ScalarType,MV,OP> (condMax_, termIterMax_,
804  relRhsErr_, relMatErr_));
805  } else {
806  convTest_->setCondLim (condMax_);
807  convTest_->setTermIterMax (termIterMax_);
808  convTest_->setRelRhsErr (relRhsErr_);
809  convTest_->setRelMatErr (relMatErr_);
810  }
811  }
812 
813  // Create the status test for maximum number of iterations if
814  // necessary. Otherwise, update it with the new maximum iteration
815  // count.
816  if (maxIterTest_.is_null()) {
817  maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
818  } else {
819  maxIterTest_->setMaxIters (maxIters_);
820  }
821 
822  // The stopping criterion is an OR combination of the test for
823  // maximum number of iterations, and the LSQR convergence test.
824  // ("OR combination" means that both tests will always be evaluated,
825  // as opposed to a SEQ combination.)
826  typedef StatusTestCombo<ScalarType,MV,OP> combo_type;
827  // If sTest_ is not null, then maxIterTest_ and convTest_ were
828  // already constructed on entry to this routine, and sTest_ has
829  // their pointers. Thus, maxIterTest_ and convTest_ have gotten any
830  // parameter changes, so we don't need to do anything to sTest_.
831  if (sTest_.is_null()) {
832  sTest_ = rcp (new combo_type (combo_type::OR, maxIterTest_, convTest_));
833  }
834 
835  if (outputTest_.is_null ()) {
836  // Create the status test output class.
837  // This class manages and formats the output from the status test.
838  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
839  outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_,
840  Passed + Failed + Undefined);
841  // Set the solver string for the output test.
842  const std::string solverDesc = " LSQR ";
843  outputTest_->setSolverDesc (solverDesc);
844  } else {
845  // FIXME (mfh 18 Sep 2011) How do we change the output style of
846  // outputTest_, without destroying and recreating it?
847  outputTest_->setOutputManager (printer_);
848  outputTest_->setChild (sTest_);
849  outputTest_->setOutputFrequency (outputFreq_);
850  // Since outputTest_ can only be created here, I'm assuming that
851  // the fourth constructor argument ("printStates") was set
852  // correctly on constrution; I don't need to reset it (and I can't
853  // set it anyway, given StatusTestOutput's interface).
854  }
855 
856  // At this point, params is a valid ParameterList. Now we can
857  // "commit" it to our instance's ParameterList.
858  params_ = params;
859 
860  // Inform the solver manager that the current parameters were set.
861  isSet_ = true;
862 }
863 
864 
865 template<class ScalarType, class MV, class OP>
868 {
869  using Teuchos::RCP;
870  using Teuchos::rcp;
871 
872  // Set the current parameters if they were not set before. NOTE:
873  // This may occur if the user generated the solver manager with the
874  // default constructor, but did not set any parameters using
875  // setParameters().
876  if (! isSet_) {
877  this->setParameters (Teuchos::parameterList (* (getValidParameters ())));
878  }
879 
881  (problem_.is_null (), LSQRSolMgrLinearProblemFailure,
882  "Belos::LSQRSolMgr::solve: The linear problem to solve is null.");
884  (! problem_->isProblemSet (), LSQRSolMgrLinearProblemFailure,
885  "Belos::LSQRSolMgr::solve: The linear problem is not ready, "
886  "as its setProblem() method has not been called.");
888  (MVT::GetNumberVecs (*(problem_->getRHS ())) != 1,
889  LSQRSolMgrBlockSizeFailure, "Belos::LSQRSolMgr::solve: "
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.");
894 
895  // We've validated the LinearProblem instance above. If any of the
896  // StatusTests needed to be initialized using information from the
897  // LinearProblem, now would be the time to do so. (This is true of
898  // GMRES, where the residual convergence test(s) to instantiate
899  // depend on knowing whether there is a left preconditioner. This
900  // is why GMRES has an "isSTSet_" Boolean member datum, which tells
901  // you whether the status tests have been instantiated and are ready
902  // for use.
903 
904  // test isFlexible might go here.
905 
906  // Next the right-hand sides to solve are identified. Among other things,
907  // this enables getCurrLHSVec() to get the current initial guess vector,
908  // and getCurrRHSVec() to get the current right-hand side (in Iter).
909  std::vector<int> currRHSIdx (1, 0);
910  problem_->setLSIndex (currRHSIdx);
911 
912  // Reset the status test.
913  outputTest_->reset ();
914 
915  // Don't assume convergence unless we've verified that the
916  // convergence test passed.
917  bool isConverged = false;
918 
919  // FIXME: Currently we are setting the initial guess to zero, since
920  // the solver doesn't yet know how to handle a nonzero initial
921  // guess. This could be fixed by rewriting the solver to work with
922  // the residual and a delta.
923  //
924  // In a least squares problem with a nonzero initial guess, the
925  // minimzation problem involves the distance (in a norm depending on
926  // the preconditioner) between the solution and the the initial
927  // guess.
928 
930  // Solve the linear problem using LSQR
932 
933  // Parameter list for the LSQR iteration.
935 
936  // Use the solver manager's "Lambda" parameter to set the
937  // iteration's "Lambda" parameter. We know that the solver
938  // manager's parameter list (params_) does indeed contain the
939  // "Lambda" parameter, because solve() always ensures that
940  // setParameters() has been called.
941  plist.set ("Lambda", lambda_);
942 
943  typedef LSQRIter<ScalarType,MV,OP> iter_type;
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_);
948 #endif
949 
950  // Reset the number of iterations.
951  lsqr_iter->resetNumIters ();
952  // Reset the number of calls that the status test output knows about.
953  outputTest_->resetNumCalls ();
954  // Set the new state and initialize the solver.
956  lsqr_iter->initializeLSQR (newstate);
957  // tell lsqr_iter to iterate
958  try {
959  lsqr_iter->iterate ();
960 
961  // First check for convergence. If we didn't converge, then check
962  // whether we reached the maximum number of iterations. If
963  // neither of those happened, there must have been a bug.
964  if (convTest_->getStatus () == Belos::Passed) {
965  isConverged = true;
966  } else if (maxIterTest_->getStatus () == Belos::Passed) {
967  isConverged = false;
968  } else {
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.");
974  }
975  } catch (const std::exception& e) {
976  printer_->stream(Belos::Errors)
977  << "Error! Caught std::exception in LSQRIter::iterate at iteration "
978  << lsqr_iter->getNumIters () << std::endl << e.what () << std::endl;
979  throw;
980  }
981 
982  // identify current linear system as solved LinearProblem
983  problem_->setCurrLS();
984  // print final summary
985  sTest_->print (printer_->stream (Belos::FinalSummary));
986 
987  // Print timing information, if the corresponding compile-time and
988  // run-time options are enabled.
989 #ifdef BELOS_TEUCHOS_TIME_MONITOR
990  // Calling summarize() can be expensive, so don't call unless the
991  // user wants to print out timing details. summarize() will do all
992  // the work even if it's passed a "black hole" output stream.
993  if (verbosity_ & TimingDetails)
995 #endif // BELOS_TEUCHOS_TIME_MONITOR
996 
997  // A posteriori solve information
998  numIters_ = maxIterTest_->getNumIters();
999  matCondNum_ = convTest_->getMatCondNum();
1000  matNorm_ = convTest_->getMatNorm();
1001  resNorm_ = convTest_->getResidNorm();
1002  matResNorm_ = convTest_->getLSResidNorm();
1003 
1004  if (! isConverged) {
1005  return Belos::Unconverged;
1006  } else {
1007  return Belos::Converged;
1008  }
1009 }
1010 
1011 // LSQRSolMgr requires the solver manager to return an eponymous std::string.
1012 template<class ScalarType, class MV, class OP>
1014 {
1015  std::ostringstream oss;
1016  oss << "LSQRSolMgr<...," << STS::name () << ">";
1017  oss << "{";
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_;
1024  oss << "}";
1025  return oss.str ();
1026 }
1027 
1028 } // end Belos namespace
1029 
1030 #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...
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.
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.
T & get(const std::string &name, T def_value)
MsgType
Available message types recognized by the linear solvers.
Definition: BelosTypes.hpp:262
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.
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...
bool isParameter(const std::string &name) const
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
LSQRSolMgrOrthoFailure(const std::string &what_arg)
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. ...
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)
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.
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:155
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
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.
TypeTo as(const TypeFrom &t)
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...
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Get current linear problem being solved for in this object.
LSQRSolMgrLinearProblemFailure(const std::string &what_arg)
bool is_null() const

Generated on Fri Aug 14 2020 10:48:33 for Belos by doxygen 1.8.5