Belos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
53
55#include "BelosLSQRIter.hpp"
61#include "Teuchos_as.hpp"
62
63#ifdef BELOS_TEUCHOS_TIME_MONITOR
64#include "Teuchos_TimeMonitor.hpp"
65#endif
66
67namespace Belos {
68
69
71
72
80public:
81 LSQRSolMgrLinearProblemFailure(const std::string& what_arg)
82 : BelosError(what_arg)
83 {}
84};
85
93public:
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.
217template<class ScalarType, class MV, class OP,
218 const bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex>
220 public Details::RealSolverManager<ScalarType, MV, OP,
221 Teuchos::ScalarTraits<ScalarType>::isComplex>
222{
223 static const bool isComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
225
226public:
228 base_type ()
229 {}
230 LSQRSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
231 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
232 base_type ()
233 {}
234 virtual ~LSQRSolMgr () {}
235
237 Teuchos::RCP<SolverManager<ScalarType, MV, OP> > clone () const override {
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.
246template<class ScalarType, class MV, class OP>
247class LSQRSolMgr<ScalarType, MV, OP, false> :
248 public Details::RealSolverManager<ScalarType, MV, OP, false> {
249private:
252 typedef Teuchos::ScalarTraits<ScalarType> STS;
253 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
254 typedef Teuchos::ScalarTraits<MagnitudeType> STM;
255
256public:
257
259
260
267 LSQRSolMgr ();
268
296 LSQRSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >& problem,
297 const Teuchos::RCP<Teuchos::ParameterList>& pl);
298
300 virtual ~LSQRSolMgr () {}
301
303 Teuchos::RCP<SolverManager<ScalarType, MV, OP> > clone () const override {
304 return Teuchos::rcp(new LSQRSolMgr<ScalarType,MV,OP>);
305 }
307
309
313 return *problem_;
314 }
315
318 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const override;
319
322 Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const override {
323 return params_;
324 }
325
331 Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers () const {
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
443private:
444
446 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
448 Teuchos::RCP<OutputManager<ScalarType> > printer_;
450 Teuchos::RCP<std::ostream> outputStream_;
451
453 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
454 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
455 Teuchos::RCP<LSQRStatusTest<ScalarType,MV,OP> > convTest_;
456 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
457
459 Teuchos::RCP<Teuchos::ParameterList> params_;
460
466 mutable Teuchos::RCP<const Teuchos::ParameterList> validParams_;
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_;
485 Teuchos::RCP<Teuchos::Time> timerSolve_;
486
487 // Internal state variables.
488 bool isSet_;
490};
491
492template<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
512template<class ScalarType, class MV, class OP>
514LSQRSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >& problem,
515 const Teuchos::RCP<Teuchos::ParameterList>& pl) :
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
547template<class ScalarType, class MV, class OP>
548Teuchos::RCP<const Teuchos::ParameterList>
550{
551 using Teuchos::ParameterList;
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
607template<class ScalarType, class MV, class OP>
608void
610setParameters (const Teuchos::RCP<Teuchos::ParameterList>& params)
611{
612 using Teuchos::isParameterType;
613 using Teuchos::getParameter;
614 using Teuchos::null;
615 using Teuchos::ParameterList;
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;
623 using Teuchos::Exceptions::InvalidParameter;
624 using Teuchos::Exceptions::InvalidParameterName;
625 using Teuchos::Exceptions::InvalidParameterType;
626
627 TEUCHOS_TEST_FOR_EXCEPTION
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 }
661 TEUCHOS_TEST_FOR_EXCEPTION
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");
714 } catch (Teuchos::Exceptions::InvalidParameterType&) {
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_,
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
848template<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
863 TEUCHOS_TEST_FOR_EXCEPTION
864 (problem_.is_null (), LSQRSolMgrLinearProblemFailure,
865 "Belos::LSQRSolMgr::solve: The linear problem to solve is null.");
866 TEUCHOS_TEST_FOR_EXCEPTION
867 (! problem_->isProblemSet (), LSQRSolMgrLinearProblemFailure,
868 "Belos::LSQRSolMgr::solve: The linear problem is not ready, "
869 "as its setProblem() method has not been called.");
870 TEUCHOS_TEST_FOR_EXCEPTION
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.
917 Teuchos::ParameterList plist;
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 {
952 TEUCHOS_TEST_FOR_EXCEPTION
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)
977 Teuchos::TimeMonitor::summarize (printer_->stream (Belos::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.
995template<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 */
Belos header file which uses auto-configuration information to include necessary C++ headers.
Belos concrete class that iterates LSQR.
IterationState contains the data that defines the state of the LSQR solver at any given time.
Belos::StatusTest class defining LSQR convergence.
Class which describes the linear problem to be solved by the iterative solver.
Class which manages the output and verbosity of the Belos solvers.
Pure virtual base class which describes the basic interface for a solver manager.
Belos::StatusTest for logically combining several status tests.
Belos::StatusTest class for specifying a maximum number of iterations.
A factory class for generating StatusTestOutput objects.
Collection of types and exceptions used within the Belos solvers.
Parent class to all Belos exceptions.
Base class for Belos::SolverManager subclasses which normally can only compile for real ScalarType.
Implementation of the LSQR iteration.
LSQRSolMgrBlockSizeFailure is thrown when the linear problem has more than one RHS.
LSQRSolMgrBlockSizeFailure(const std::string &what_arg)
Belos::LSQRSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i....
LSQRSolMgrLinearProblemFailure(const std::string &what_arg)
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
The "master" status test (that includes all status tests).
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Get current linear problem being solved for in this object.
Teuchos::ScalarTraits< ScalarType > STS
MagnitudeType getMatCondNum() const
Estimated matrix condition number from the last solve.
Teuchos::ScalarTraits< MagnitudeType > STM
Teuchos::RCP< OutputManager< ScalarType > > printer_
The output manager.
int getNumIters() const override
Iteration count from the last solve.
virtual ~LSQRSolMgr()
Destructor (declared virtual for memory safety of base classes).
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
void reset(const ResetType type) override
reset the solver manager as specified by the ResetType, informs the solver manager that the solver sh...
Teuchos::RCP< std::ostream > outputStream_
Output stream to which to write status output.
bool isLOADetected() const override
Whether a loss of accuracy was detected during the last solve.
MagnitudeType getMatResNorm() const
Estimate of (residual vector ) from the last solve.
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::RCP< LSQRStatusTest< ScalarType, MV, OP > > convTest_
MagnitudeType getResNorm() const
Estimated residual norm from the last solve.
Teuchos::RCP< Teuchos::ParameterList > params_
Current parameter list.
MagnitudeType getMatNorm() const
Estimated matrix Frobenius norm from the last solve.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
Teuchos::RCP< const Teuchos::ParameterList > validParams_
Default parameter list.
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
The linear problem to solve.
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
LSQR method (for linear systems and linear least-squares problems).
Details::RealSolverManager< ScalarType, MV, OP, isComplex > base_type
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
static const bool isComplex
LSQRSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
A linear system to solve, and its associated information.
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
A class for extending the status testing capabilities of Belos via logical combinations.
A Belos::StatusTest class for specifying a maximum number of iterations.
A factory class for generating StatusTestOutput objects.
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.
MsgType
Available message types recognized by the linear solvers.
@ FinalSummary
@ TimingDetails
ReturnType
Whether the Belos solve converged for all linear systems.
@ Unconverged
ResetType
How to reset the solver.
Structure to contain pointers to LSQRIteration state variables, ...