Belos  Version of the Day
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 };
253 
254 
255 // Partial specialization for real ScalarType.
256 // This contains the actual working implementation of LSQR.
257 // See discussion in the class documentation above.
258 template<class ScalarType, class MV, class OP>
259 class LSQRSolMgr<ScalarType, MV, OP, false> :
260  public Details::RealSolverManager<ScalarType, MV, OP, false> {
261 private:
265  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
267 
268 public:
269 
271 
272 
279  LSQRSolMgr ();
280 
310 
312  virtual ~LSQRSolMgr () {}
313 
315 
317 
321  return *problem_;
322  }
323 
326  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
327 
331  return params_;
332  }
333 
340  return Teuchos::tuple (timerSolve_);
341  }
342 
344  int getNumIters () const {
345  return numIters_;
346  }
347 
352  MagnitudeType getMatCondNum () const {
353  return matCondNum_;
354  }
355 
360  MagnitudeType getMatNorm () const {
361  return matNorm_;
362  }
363 
372  MagnitudeType getResNorm () const {
373  return resNorm_;
374  }
375 
377  MagnitudeType getMatResNorm () const {
378  return matResNorm_;
379  }
380 
389  bool isLOADetected () const { return false; }
390 
392 
394 
395 
398  problem_ = problem;
399  }
400 
402  void setParameters (const Teuchos::RCP<Teuchos::ParameterList>& params);
403 
405 
407 
408 
412  void reset (const ResetType type) {
413  if ((type & Belos::Problem) && ! problem_.is_null ()) {
414  problem_->setProblem ();
415  }
416  }
417 
419 
421 
440  ReturnType solve();
441 
443 
445 
447  std::string description () const;
448 
450 
451 private:
452 
458  Teuchos::RCP<std::ostream> outputStream_;
459 
465 
468 
474  mutable Teuchos::RCP<const Teuchos::ParameterList> validParams_;
475 
476  // Current solver input parameters
477  MagnitudeType lambda_;
478  MagnitudeType relRhsErr_;
479  MagnitudeType relMatErr_;
480  MagnitudeType condMax_;
481  int maxIters_, termIterMax_;
482  int verbosity_, outputStyle_, outputFreq_;
483 
484  // Terminal solver state values
485  int numIters_;
486  MagnitudeType matCondNum_;
487  MagnitudeType matNorm_;
488  MagnitudeType resNorm_;
489  MagnitudeType matResNorm_;
490 
491  // Timers.
492  std::string label_;
493  Teuchos::RCP<Teuchos::Time> timerSolve_;
494 
495  // Internal state variables.
496  bool isSet_;
497  bool loaDetected_;
498 };
499 
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 ()),
506  maxIters_ (1000),
507  termIterMax_ (1),
508  verbosity_ (Belos::Errors),
509  outputStyle_ (Belos::General),
510  outputFreq_ (-1),
511  numIters_ (0),
512  matCondNum_ (STM::zero ()),
513  matNorm_ (STM::zero ()),
514  resNorm_ (STM::zero ()),
515  matResNorm_ (STM::zero ()),
516  isSet_ (false),
517  loaDetected_ (false)
518 {}
519 
520 template<class ScalarType, class MV, class OP>
524  problem_ (problem),
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 ()),
529  maxIters_ (1000),
530  termIterMax_ (1),
531  verbosity_ (Belos::Errors),
532  outputStyle_ (Belos::General),
533  outputFreq_ (-1),
534  numIters_ (0),
535  matCondNum_ (STM::zero ()),
536  matNorm_ (STM::zero ()),
537  resNorm_ (STM::zero ()),
538  matResNorm_ (STM::zero ()),
539  isSet_ (false),
540  loaDetected_ (false)
541 {
542  // The linear problem to solve is allowed to be null here. The user
543  // must then set a nonnull linear problem (by calling setProblem())
544  // before calling solve().
545  //
546  // Similarly, users are allowed to set a null parameter list here,
547  // but they must first set a nonnull parameter list (by calling
548  // setParameters()) before calling solve().
549  if (! pl.is_null ()) {
550  setParameters (pl);
551  }
552 }
553 
554 
555 template<class ScalarType, class MV, class OP>
558 {
560  using Teuchos::parameterList;
561  using Teuchos::RCP;
562  using Teuchos::rcp;
563  using Teuchos::rcpFromRef;
564 
565  // Set all the valid parameters and their default values.
566  if (validParams_.is_null ()) {
567  // We use Teuchos::as just in case MagnitudeType doesn't have a
568  // constructor that takes an int. Otherwise, we could just write
569  // "MagnitudeType(10)".
570  const MagnitudeType ten = Teuchos::as<MagnitudeType> (10);
571  const MagnitudeType sqrtEps = STM::squareroot (STM::eps());
572 
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;
580  const int verbosity = Belos::Errors;
581  const int outputStyle = Belos::General;
582  const int outputFreq = -1;
583  const std::string label ("Belos");
584 
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 "
595  "number of Abar.");
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 "
606  "labels");
607  // pl->set("Restart Timers", restartTimers_);
608  pl->set ("Block Size", 1, "Block size parameter (currently, LSQR requires "
609  "this must always be 1)");
610  validParams_ = pl;
611  }
612  return validParams_;
613 }
614 
615 
616 template<class ScalarType, class MV, class OP>
617 void
620 {
621  using Teuchos::isParameterType;
622  using Teuchos::getParameter;
623  using Teuchos::null;
625  using Teuchos::parameterList;
626  using Teuchos::RCP;
627  using Teuchos::rcp;
628  using Teuchos::rcp_dynamic_cast;
629  using Teuchos::rcpFromRef;
630  using Teuchos::Time;
631  using Teuchos::TimeMonitor;
635 
637  (params.is_null (), std::invalid_argument,
638  "Belos::LSQRSolMgr::setParameters: The input ParameterList is null.");
639  RCP<const ParameterList> defaultParams = getValidParameters ();
640 
641  // FIXME (mfh 29 Apr 2015) Our users would like to supply one
642  // ParameterList that works for both GMRES and LSQR. Thus, we want
643  // LSQR (the less-used solver) to ignore parameters it doesn't
644  // recognize). For now, therefore, it should not validate, since
645  // validation cannot distinguish between misspellings and
646  // unrecognized parameters. (Perhaps Belos should have a central
647  // facility for all parameters recognized by some solver in Belos,
648  // so we could use that for spell checking.)
649  //
650  //params->validateParameters (*defaultParams);
651 
652  // mfh 29 Apr 2015: The convention in Belos is that the input
653  // ParameterList is a "delta" from the current state. Thus, we
654  // don't fill in the input ParameterList with defaults, and we only
655  // change the current state if the corresponding parameter was
656  // explicitly set in the input ParameterList. We set up the solver
657  // with the default state on construction.
658 
659  // Get the damping (regularization) parameter lambda.
660  if (params->isParameter ("Lambda")) {
661  lambda_ = params->get<MagnitudeType> ("Lambda");
662  } else if (params->isParameter ("lambda")) {
663  lambda_ = params->get<MagnitudeType> ("lambda");
664  }
665 
666  // Get the maximum number of iterations.
667  if (params->isParameter ("Maximum Iterations")) {
668  maxIters_ = params->get<int> ("Maximum Iterations");
669  }
671  (maxIters_ < 0, std::invalid_argument, "Belos::LSQRSolMgr::setParameters: "
672  "\"Maximum Iterations\" = " << maxIters_ << " < 0.");
673 
674  // (Re)set the timer label.
675  {
676  const std::string newLabel =
677  params->isParameter ("Timer Label") ?
678  params->get<std::string> ("Timer Label") :
679  label_;
680 
681  // Update parameter in our list and solver timer
682  if (newLabel != label_) {
683  label_ = newLabel;
684  }
685 
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 ()) {
691  // Ask TimeMonitor for a new timer.
692  timerSolve_ = TimeMonitor::getNewCounter (newSolveLabel);
693  } else {
694  // We've already created a timer, but we may have changed its
695  // label. If we did change its name, then we have to forget
696  // about the old timer and create a new one with a different
697  // name. This is because Teuchos::Time doesn't give you a way
698  // to change a timer's name, once you've created it. We assume
699  // that if the user changed the timer's label, then the user
700  // wants to reset the timer's results.
701  const std::string oldSolveLabel = timerSolve_->name ();
702 
703  if (oldSolveLabel != newSolveLabel) {
704  // Tell TimeMonitor to forget about the old timer.
705  // TimeMonitor lets you clear timers by name.
706  TimeMonitor::clearCounter (oldSolveLabel);
707  timerSolve_ = TimeMonitor::getNewCounter (newSolveLabel);
708  }
709  }
710 #endif // BELOS_TEUCHOS_TIME_MONITOR
711  }
712 
713  // Check for a change in verbosity level
714  if (params->isParameter ("Verbosity")) {
715  int newVerbosity = 0;
716  // ParameterList gets confused sometimes about enums. This
717  // ensures that no matter how "Verbosity" was stored -- either an
718  // an int, or as a Belos::MsgType enum, we will be able to extract
719  // it. If it was stored as some other type, we let the exception
720  // through.
721  try {
722  newVerbosity = params->get<Belos::MsgType> ("Verbosity");
724  newVerbosity = params->get<int> ("Verbosity");
725  }
726  if (newVerbosity != verbosity_) {
727  verbosity_ = newVerbosity;
728  }
729  }
730 
731  // (Re)set the output style.
732  if (params->isParameter ("Output Style")) {
733  outputStyle_ = params->get<int> ("Output Style");
734  }
735 
736  // Get the output stream for the output manager.
737  //
738  // In case the output stream can't be read back in, we default to
739  // stdout (std::cout), just to ensure reasonable behavior.
740  if (params->isParameter ("Output Stream")) {
741  outputStream_ = params->get<RCP<std::ostream> > ("Output Stream");
742  }
743  // We assume that a null output stream indicates that the user
744  // doesn't want to print anything, so we replace it with a "black
745  // hole" stream that prints nothing sent to it. (We can't use a
746  // null output stream, since the output manager always sends
747  // things it wants to print to the output stream.)
748  if (outputStream_.is_null ()) {
749  outputStream_ = rcp (new Teuchos::oblackholestream ());
750  }
751 
752  // Get the frequency of solver output. (For example, -1 means
753  // "never," and 1 means "every iteration.")
754  if (params->isParameter ("Output Frequency")) {
755  outputFreq_ = params->get<int> ("Output Frequency");
756  }
757 
758  // Create output manager if we need to, using the verbosity level
759  // and output stream that we fetched above. Status tests (i.e.,
760  // stopping criteria) need this.
761  if (printer_.is_null ()) {
762  printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
763  } else {
764  printer_->setVerbosity (verbosity_);
765  printer_->setOStream (outputStream_);
766  }
767 
768  // Check for condition number limit, number of consecutive passed
769  // iterations, relative RHS error, and relative matrix error.
770  // Create the LSQR convergence test if necessary.
771  {
772  if (params->isParameter ("Condition Limit")) {
773  condMax_ = params->get<MagnitudeType> ("Condition Limit");
774  }
775  if (params->isParameter ("Term Iter Max")) {
776  termIterMax_ = params->get<int> ("Term Iter Max");
777  }
778  if (params->isParameter ("Rel RHS Err")) {
779  relRhsErr_ = params->get<MagnitudeType> ("Rel RHS Err");
780  }
781  else if (params->isParameter ("Convergence Tolerance")) {
782  // NOTE (mfh 29 Apr 2015) We accept this parameter as an alias
783  // for "Rel RHS Err".
784  relRhsErr_ = params->get<MagnitudeType> ("Convergence Tolerance");
785  }
786 
787  if (params->isParameter ("Rel Mat Err")) {
788  relMatErr_ = params->get<MagnitudeType> ("Rel Mat Err");
789  }
790 
791  // Create the LSQR convergence test if it doesn't exist yet.
792  // Otherwise, update its parameters.
793  if (convTest_.is_null ()) {
794  convTest_ =
795  rcp (new LSQRStatusTest<ScalarType,MV,OP> (condMax_, termIterMax_,
796  relRhsErr_, relMatErr_));
797  } else {
798  convTest_->setCondLim (condMax_);
799  convTest_->setTermIterMax (termIterMax_);
800  convTest_->setRelRhsErr (relRhsErr_);
801  convTest_->setRelMatErr (relMatErr_);
802  }
803  }
804 
805  // Create the status test for maximum number of iterations if
806  // necessary. Otherwise, update it with the new maximum iteration
807  // count.
808  if (maxIterTest_.is_null()) {
809  maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
810  } else {
811  maxIterTest_->setMaxIters (maxIters_);
812  }
813 
814  // The stopping criterion is an OR combination of the test for
815  // maximum number of iterations, and the LSQR convergence test.
816  // ("OR combination" means that both tests will always be evaluated,
817  // as opposed to a SEQ combination.)
818  typedef StatusTestCombo<ScalarType,MV,OP> combo_type;
819  // If sTest_ is not null, then maxIterTest_ and convTest_ were
820  // already constructed on entry to this routine, and sTest_ has
821  // their pointers. Thus, maxIterTest_ and convTest_ have gotten any
822  // parameter changes, so we don't need to do anything to sTest_.
823  if (sTest_.is_null()) {
824  sTest_ = rcp (new combo_type (combo_type::OR, maxIterTest_, convTest_));
825  }
826 
827  if (outputTest_.is_null ()) {
828  // Create the status test output class.
829  // This class manages and formats the output from the status test.
830  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
831  outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_,
832  Passed + Failed + Undefined);
833  // Set the solver string for the output test.
834  const std::string solverDesc = " LSQR ";
835  outputTest_->setSolverDesc (solverDesc);
836  } else {
837  // FIXME (mfh 18 Sep 2011) How do we change the output style of
838  // outputTest_, without destroying and recreating it?
839  outputTest_->setOutputManager (printer_);
840  outputTest_->setChild (sTest_);
841  outputTest_->setOutputFrequency (outputFreq_);
842  // Since outputTest_ can only be created here, I'm assuming that
843  // the fourth constructor argument ("printStates") was set
844  // correctly on constrution; I don't need to reset it (and I can't
845  // set it anyway, given StatusTestOutput's interface).
846  }
847 
848  // At this point, params is a valid ParameterList. Now we can
849  // "commit" it to our instance's ParameterList.
850  params_ = params;
851 
852  // Inform the solver manager that the current parameters were set.
853  isSet_ = true;
854 }
855 
856 
857 template<class ScalarType, class MV, class OP>
860 {
861  using Teuchos::RCP;
862  using Teuchos::rcp;
863 
864  // Set the current parameters if they were not set before. NOTE:
865  // This may occur if the user generated the solver manager with the
866  // default constructor, but did not set any parameters using
867  // setParameters().
868  if (! isSet_) {
869  this->setParameters (Teuchos::parameterList (* (getValidParameters ())));
870  }
871 
873  (problem_.is_null (), LSQRSolMgrLinearProblemFailure,
874  "Belos::LSQRSolMgr::solve: The linear problem to solve is null.");
876  (! problem_->isProblemSet (), LSQRSolMgrLinearProblemFailure,
877  "Belos::LSQRSolMgr::solve: The linear problem is not ready, "
878  "as its setProblem() method has not been called.");
880  (MVT::GetNumberVecs (*(problem_->getRHS ())) != 1,
881  LSQRSolMgrBlockSizeFailure, "Belos::LSQRSolMgr::solve: "
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.");
886 
887  // We've validated the LinearProblem instance above. If any of the
888  // StatusTests needed to be initialized using information from the
889  // LinearProblem, now would be the time to do so. (This is true of
890  // GMRES, where the residual convergence test(s) to instantiate
891  // depend on knowing whether there is a left preconditioner. This
892  // is why GMRES has an "isSTSet_" Boolean member datum, which tells
893  // you whether the status tests have been instantiated and are ready
894  // for use.
895 
896  // test isFlexible might go here.
897 
898  // Next the right-hand sides to solve are identified. Among other things,
899  // this enables getCurrLHSVec() to get the current initial guess vector,
900  // and getCurrRHSVec() to get the current right-hand side (in Iter).
901  std::vector<int> currRHSIdx (1, 0);
902  problem_->setLSIndex (currRHSIdx);
903 
904  // Reset the status test.
905  outputTest_->reset ();
906 
907  // Don't assume convergence unless we've verified that the
908  // convergence test passed.
909  bool isConverged = false;
910 
911  // FIXME: Currently we are setting the initial guess to zero, since
912  // the solver doesn't yet know how to handle a nonzero initial
913  // guess. This could be fixed by rewriting the solver to work with
914  // the residual and a delta.
915  //
916  // In a least squares problem with a nonzero initial guess, the
917  // minimzation problem involves the distance (in a norm depending on
918  // the preconditioner) between the solution and the the initial
919  // guess.
920 
922  // Solve the linear problem using LSQR
924 
925  // Parameter list for the LSQR iteration.
927 
928  // Use the solver manager's "Lambda" parameter to set the
929  // iteration's "Lambda" parameter. We know that the solver
930  // manager's parameter list (params_) does indeed contain the
931  // "Lambda" parameter, because solve() always ensures that
932  // setParameters() has been called.
933  plist.set ("Lambda", lambda_);
934 
935  typedef LSQRIter<ScalarType,MV,OP> iter_type;
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_);
940 #endif
941 
942  // Reset the number of iterations.
943  lsqr_iter->resetNumIters ();
944  // Reset the number of calls that the status test output knows about.
945  outputTest_->resetNumCalls ();
946  // Set the new state and initialize the solver.
948  lsqr_iter->initializeLSQR (newstate);
949  // tell lsqr_iter to iterate
950  try {
951  lsqr_iter->iterate ();
952 
953  // First check for convergence. If we didn't converge, then check
954  // whether we reached the maximum number of iterations. If
955  // neither of those happened, there must have been a bug.
956  if (convTest_->getStatus () == Belos::Passed) {
957  isConverged = true;
958  } else if (maxIterTest_->getStatus () == Belos::Passed) {
959  isConverged = false;
960  } else {
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.");
966  }
967  } catch (const std::exception& e) {
968  printer_->stream(Belos::Errors)
969  << "Error! Caught std::exception in LSQRIter::iterate at iteration "
970  << lsqr_iter->getNumIters () << std::endl << e.what () << std::endl;
971  throw;
972  }
973 
974  // identify current linear system as solved LinearProblem
975  problem_->setCurrLS();
976  // print final summary
977  sTest_->print (printer_->stream (Belos::FinalSummary));
978 
979  // Print timing information, if the corresponding compile-time and
980  // run-time options are enabled.
981 #ifdef BELOS_TEUCHOS_TIME_MONITOR
982  // Calling summarize() can be expensive, so don't call unless the
983  // user wants to print out timing details. summarize() will do all
984  // the work even if it's passed a "black hole" output stream.
985  if (verbosity_ & TimingDetails)
987 #endif // BELOS_TEUCHOS_TIME_MONITOR
988 
989  // A posteriori solve information
990  numIters_ = maxIterTest_->getNumIters();
991  matCondNum_ = convTest_->getMatCondNum();
992  matNorm_ = convTest_->getMatNorm();
993  resNorm_ = convTest_->getResidNorm();
994  matResNorm_ = convTest_->getLSResidNorm();
995 
996  if (! isConverged) {
997  return Belos::Unconverged;
998  } else {
999  return Belos::Converged;
1000  }
1001 }
1002 
1003 // LSQRSolMgr requires the solver manager to return an eponymous std::string.
1004 template<class ScalarType, class MV, class OP>
1006 {
1007  std::ostringstream oss;
1008  oss << "LSQRSolMgr<...," << STS::name () << ">";
1009  oss << "{";
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_;
1016  oss << "}";
1017  return oss.str ();
1018 }
1019 
1020 } // end Belos namespace
1021 
1022 #endif /* BELOS_LSQR_SOLMGR_HPP */
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&#39;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.
T & get(const std::string &name, T def_value)
MsgType
Available message types recognized by the linear solvers.
Definition: BelosTypes.hpp:257
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.
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.
bool is_null() const
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.
Definition: BelosTypes.hpp:205
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)
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.
Definition: BelosTypes.hpp:154
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.
TypeTo as(const TypeFrom &t)
A class for extending the status testing capabilities of Belos via logical combinations.
bool isParameter(const std::string &name) const
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...
LSQRSolMgrLinearProblemFailure(const std::string &what_arg)

Generated for Belos by doxygen 1.8.14