Belos Version of the Day
BelosLSQRSolMgr.hpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 //
00004 //                 Belos: Block Linear Solvers Package
00005 //                  Copyright 2004 Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 #ifndef BELOS_LSQR_SOLMGR_HPP
00043 #define BELOS_LSQR_SOLMGR_HPP
00044 
00047 
00048 #include "BelosConfigDefs.hpp"
00049 #include "BelosTypes.hpp"
00050 
00051 #include "BelosLinearProblem.hpp"
00052 #include "BelosSolverManager.hpp"
00053 
00054 #include "BelosLSQRIteration.hpp"
00055 #include "BelosLSQRIter.hpp"
00056 #include "BelosOrthoManagerFactory.hpp"
00057 #include "BelosStatusTestMaxIters.hpp"
00058 #include "BelosLSQRStatusTest.hpp"
00059 #include "BelosStatusTestCombo.hpp"
00060 #include "BelosStatusTestOutputFactory.hpp"
00061 #include "BelosOutputManager.hpp"
00062 #include "Teuchos_as.hpp"
00063 #include "Teuchos_BLAS.hpp"
00064 #include "Teuchos_LAPACK.hpp"
00065 
00066 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00067 #include "Teuchos_TimeMonitor.hpp"
00068 #endif
00069 
00070 namespace Belos {
00071 
00072 
00074 
00075 
00082 class LSQRSolMgrLinearProblemFailure : public BelosError {
00083 public:
00084   LSQRSolMgrLinearProblemFailure(const std::string& what_arg)
00085     : BelosError(what_arg)
00086   {}
00087 };
00088 
00095 class LSQRSolMgrOrthoFailure : public BelosError {
00096 public:
00097   LSQRSolMgrOrthoFailure(const std::string& what_arg)
00098     : BelosError(what_arg)
00099   {}
00100 };
00101 
00108 class LSQRSolMgrBlockSizeFailure : public BelosError {
00109 public:
00110   LSQRSolMgrBlockSizeFailure(const std::string& what_arg)
00111   : BelosError(what_arg)
00112   {}
00113 };
00114 
00228 
00229 
00230 // Partial specialization for complex ScalarType.
00231 // This contains a trivial implementation.
00232 // See discussion in the class documentation above.
00233 template<class ScalarType, class MV, class OP,
00234          const bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex>
00235 class LSQRSolMgr :
00236     public Details::RealSolverManager<ScalarType, MV, OP,
00237                                       Teuchos::ScalarTraits<ScalarType>::isComplex>
00238 {
00239   static const bool isComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
00240   typedef Details::RealSolverManager<ScalarType, MV, OP, isComplex> base_type;
00241 
00242 public:
00243   LSQRSolMgr () :
00244     base_type ()
00245   {}
00246   LSQRSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00247               const Teuchos::RCP<Teuchos::ParameterList> &pl) :
00248     base_type ()
00249   {}
00250   virtual ~LSQRSolMgr () {}
00251 };
00252 
00253 
00254 // Partial specialization for real ScalarType.
00255 // This contains the actual working implementation of LSQR.
00256 // See discussion in the class documentation above.
00257 template<class ScalarType, class MV, class OP>
00258 class LSQRSolMgr<ScalarType, MV, OP, false> :
00259     public Details::RealSolverManager<ScalarType, MV, OP, false> {
00260 private:
00261   typedef MultiVecTraits<ScalarType,MV> MVT;
00262   typedef OperatorTraits<ScalarType,MV,OP> OPT;
00263   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00264   typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00265   typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00266 
00267 public:
00268 
00270 
00271 
00278   LSQRSolMgr ();
00279 
00302   // This LSQR implementation only supports block size 1.  Like CG, LSQR is a short
00303   // recurrence method that, in finite precision arithmetic and without reorthogonalization,
00304   // does not have the "n" step convergence property.  Without either blocks or
00305   // reorthogonalization, there is nothing to "Orthogonalize."
00306 
00307   LSQRSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >& problem,
00308               const Teuchos::RCP<Teuchos::ParameterList>& pl);
00309 
00311   virtual ~LSQRSolMgr () {};
00312 
00314 
00315 
00316 
00319   const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00320     return *problem_;
00321   }
00322 
00325   Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00326 
00329   Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
00330 
00336   Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00337     return Teuchos::tuple(timerSolve_);
00338   }
00339 
00341   int getNumIters() const {
00342     return numIters_;
00343   }
00344 
00349   MagnitudeType getMatCondNum () const {
00350     return matCondNum_;
00351   }
00352 
00357   MagnitudeType getMatNorm () const {
00358     return matNorm_;
00359   }
00360 
00369   MagnitudeType getResNorm () const {
00370     return resNorm_;
00371   }
00372 
00374   MagnitudeType getMatResNorm () const {
00375     return matResNorm_;
00376   }
00377 
00386   bool isLOADetected () const { return false; }
00387 
00389 
00391 
00392 
00394   void setProblem (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >& problem) {
00395     problem_ = problem;
00396   }
00397 
00399   void setParameters (const Teuchos::RCP<Teuchos::ParameterList>& params);
00400 
00402 
00404 
00405 
00409   void reset (const ResetType type) {
00410     if ((type & Belos::Problem) && ! problem_.is_null ()) {
00411       problem_->setProblem ();
00412     }
00413   }
00414 
00416 
00417 
00418 
00437   ReturnType solve();
00438 
00440 
00441 
00442 
00444   std::string description() const;
00445 
00447 
00448 private:
00449 
00451   Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00453   Teuchos::RCP<OutputManager<ScalarType> > printer_;
00455   Teuchos::RCP<std::ostream> outputStream_;
00456 
00458   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00459   Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00460   Teuchos::RCP<LSQRStatusTest<ScalarType,MV,OP> > convTest_;
00461   Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00462 
00464   Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
00465 
00467   Teuchos::RCP<Teuchos::ParameterList> params_;
00468 
00474   mutable Teuchos::RCP<const Teuchos::ParameterList> validParams_;
00475 
00476   // Current solver input parameters
00477   MagnitudeType lambda_;
00478   MagnitudeType relRhsErr_;
00479   MagnitudeType relMatErr_;
00480   MagnitudeType condMax_;
00481   int maxIters_, termIterMax_;
00482   std::string orthoType_;
00483   MagnitudeType orthoKappa_;
00484   int verbosity_, outputStyle_, outputFreq_;
00485 
00486   // Terminal solver state values
00487   int numIters_;
00488   MagnitudeType matCondNum_;
00489   MagnitudeType matNorm_;
00490   MagnitudeType resNorm_;
00491   MagnitudeType matResNorm_;
00492 
00493   // Timers.
00494   std::string label_;
00495   Teuchos::RCP<Teuchos::Time> timerSolve_;
00496 
00497   // Internal state variables.
00498   bool isSet_;
00499   bool loaDetected_;
00500 };
00501 
00502 // Empty Constructor
00503 template<class ScalarType, class MV, class OP>
00504 LSQRSolMgr<ScalarType,MV,OP,false>::LSQRSolMgr() :
00505   lambda_(Teuchos::ScalarTraits<MagnitudeType>::zero ()),
00506   relRhsErr_(Teuchos::ScalarTraits<MagnitudeType>::zero ()),
00507   relMatErr_(Teuchos::ScalarTraits<MagnitudeType>::zero ()),
00508   condMax_(Teuchos::ScalarTraits<MagnitudeType>::zero ()),
00509   maxIters_(0),
00510   termIterMax_(0),
00511   orthoKappa_(Teuchos::ScalarTraits<MagnitudeType>::zero ()),
00512   verbosity_(0),
00513   outputStyle_(0),
00514   outputFreq_(0),
00515   numIters_(0),
00516   matCondNum_(Teuchos::ScalarTraits<MagnitudeType>::zero ()),
00517   matNorm_(Teuchos::ScalarTraits<MagnitudeType>::zero ()),
00518   resNorm_(Teuchos::ScalarTraits<MagnitudeType>::zero ()),
00519   matResNorm_(Teuchos::ScalarTraits<MagnitudeType>::zero ()),
00520   isSet_(false),
00521   loaDetected_(false)
00522 {}
00523 
00524 
00525 // Basic Constructor
00526 template<class ScalarType, class MV, class OP>
00527 LSQRSolMgr<ScalarType,MV,OP,false>::
00528 LSQRSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00529             const Teuchos::RCP<Teuchos::ParameterList> &pl) :
00530   problem_(problem),
00531   lambda_(Teuchos::ScalarTraits<MagnitudeType>::zero ()),
00532   relRhsErr_(Teuchos::ScalarTraits<MagnitudeType>::zero ()),
00533   relMatErr_(Teuchos::ScalarTraits<MagnitudeType>::zero ()),
00534   condMax_(Teuchos::ScalarTraits<MagnitudeType>::zero ()),
00535   maxIters_(0),
00536   termIterMax_(0),
00537   orthoKappa_(Teuchos::ScalarTraits<MagnitudeType>::zero ()),
00538   verbosity_(0),
00539   outputStyle_(0),
00540   outputFreq_(0),
00541   numIters_(0),
00542   matCondNum_(Teuchos::ScalarTraits<MagnitudeType>::zero ()),
00543   matNorm_(Teuchos::ScalarTraits<MagnitudeType>::zero ()),
00544   resNorm_(Teuchos::ScalarTraits<MagnitudeType>::zero ()),
00545   matResNorm_(Teuchos::ScalarTraits<MagnitudeType>::zero ()),
00546   isSet_(false),
00547   loaDetected_(false)
00548 {
00549   // The linear problem to solve is allowed to be null here.  The user
00550   // must then set a nonnull linear problem (by calling setProblem())
00551   // before calling solve().
00552   //
00553   // Similarly, users are allowed to set a null parameter list here,
00554   // but they must first set a nonnull parameter list (by calling
00555   // setParameters()) before calling solve().
00556   if (! is_null (pl)) {
00557     setParameters (pl);
00558   }
00559 }
00560 
00561 
00562 template<class ScalarType, class MV, class OP>
00563 Teuchos::RCP<const Teuchos::ParameterList>
00564 LSQRSolMgr<ScalarType,MV,OP,false>::getValidParameters() const
00565 {
00566   using Teuchos::ParameterList;
00567   using Teuchos::parameterList;
00568   using Teuchos::RCP;
00569   using Teuchos::rcp;
00570   using Teuchos::rcpFromRef;
00571   typedef Teuchos::ScalarTraits<MagnitudeType> STM;
00572 
00573   // Set all the valid parameters and their default values.
00574   if (is_null (validParams_)) {
00575     // We use Teuchos::as just in case MagnitudeType doesn't have a
00576     // constructor that takes an int.  Otherwise, we could just write
00577     // "MagnitudeType(10)".
00578     const MagnitudeType ten = Teuchos::as<MagnitudeType> (10);
00579     const MagnitudeType sqrtEps = STM::squareroot (STM::eps());
00580 
00581     const MagnitudeType lambda = STM::zero();
00582     RCP<std::ostream> outputStream = rcpFromRef (std::cout);
00583     const MagnitudeType relRhsErr = ten * sqrtEps;
00584     const MagnitudeType relMatErr = ten * sqrtEps;
00585     const MagnitudeType condMax = STM::one() / STM::eps();
00586     const int maxIters = 1000;
00587     const int termIterMax = 1;
00588     const std::string orthoType ("DGKS");
00589     const MagnitudeType orthoKappa = Teuchos::as<MagnitudeType> (-1.0);
00590     const int verbosity = Belos::Errors;
00591     const int outputStyle = Belos::General;
00592     const int outputFreq = -1;
00593     const std::string label ("Belos");
00594 
00595     RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00596     pl->set("Output Stream", outputStream,
00597       "is a reference-counted pointer to the output stream receiving\n"
00598       "all solver output.");
00599     pl->set("Lambda", lambda, "is the damping parameter.");
00600     pl->set("Rel RHS Err", relRhsErr,
00601             "estimates the error in the data defining the right-\n"
00602             "hand side.");
00603     pl->set("Rel Mat Err", relMatErr,
00604             "estimates the error in the data defining the matrix.");
00605     pl->set("Condition Limit", condMax,
00606       "bounds the estimated condition number of Abar.");
00607     pl->set("Maximum Iterations", maxIters,
00608       "allows at most the maximum number of iterations.");
00609     pl->set("Term Iter Max", termIterMax,
00610       "consecutive iterations meeting thresholds are necessary for\n"
00611        "for convergence.");
00612     pl->set("Orthogonalization", orthoType,
00613       "uses orthogonalization of either DGKS, ICGS, IMGS, or TSQR.");
00614     {
00615       OrthoManagerFactory<ScalarType, MV, OP> factory;
00616       pl->set("Orthogonalization", orthoType,
00617               "refers to the orthogonalization method to use.  Valid "
00618               "options: " + factory.validNamesString());
00619       RCP<const ParameterList> orthoParams =
00620         factory.getDefaultParameters (orthoType);
00621       pl->set ("Orthogonalization Parameters", *orthoParams,
00622                "Parameters specific to the type of orthogonalization used.");
00623     }
00624     pl->set("Orthogonalization Constant", orthoKappa,
00625       "is the threshold used by DGKS orthogonalization to determine\n"
00626       "whether or not to repeat classical Gram-Schmidt.  This parameter\n"
00627       "is ignored if \"Orthogonalization\" is not \"DGKS\".");
00628     pl->set("Verbosity", verbosity,
00629       "type(s) of solver information are outputted to the output\n"
00630       "stream.");
00631     pl->set("Output Style", outputStyle,
00632       "the style used for the solver information outputted to the\n"
00633       "output stream.");
00634     pl->set("Output Frequency", outputFreq,
00635       "is the frequency at which information is written to the\n"
00636       "output stream.");
00637     pl->set("Timer Label", label,
00638       "is the string to use as a prefix for the timer labels.");
00639     //  pl->set("Restart Timers", restartTimers_);
00640     pl->set("Block Size", 1, "Block size parameter (currently, this must always be 1).");
00641 
00642     validParams_ = pl;
00643   }
00644   return validParams_;
00645 }
00646 
00647 
00648 template<class ScalarType, class MV, class OP>
00649 void
00650 LSQRSolMgr<ScalarType,MV,OP,false>::
00651 setParameters (const Teuchos::RCP<Teuchos::ParameterList> &params)
00652 {
00653   using Teuchos::isParameterType;
00654   using Teuchos::getParameter;
00655   using Teuchos::null;
00656   using Teuchos::ParameterList;
00657   using Teuchos::parameterList;
00658   using Teuchos::RCP;
00659   using Teuchos::rcp;
00660   using Teuchos::rcp_dynamic_cast;
00661   using Teuchos::rcpFromRef;
00662   using Teuchos::Time;
00663   using Teuchos::TimeMonitor;
00664   using Teuchos::Exceptions::InvalidParameter;
00665   using Teuchos::Exceptions::InvalidParameterName;
00666   using Teuchos::Exceptions::InvalidParameterType;
00667 
00668   TEUCHOS_TEST_FOR_EXCEPTION(params.is_null(), std::invalid_argument,
00669                      "Belos::LSQRSolMgr::setParameters: "
00670                      "the input ParameterList is null.");
00671   RCP<const ParameterList> defaultParams = getValidParameters ();
00672   params->validateParametersAndSetDefaults (*defaultParams);
00673 
00674   // At this point, params is a valid parameter list with defaults
00675   // filled in.  Now we can "commit" it to our instance's parameter
00676   // list.
00677   params_ = params;
00678 
00679   // Get the damping (a.k.a. regularization) parameter lambda.
00680   lambda_ = params->get<MagnitudeType> ("Lambda");
00681 
00682   // Get the maximum number of iterations.
00683   maxIters_ = params->get<int>("Maximum Iterations");
00684 
00685   // (Re)set the timer label.
00686   {
00687     const std::string newLabel = params->get<std::string>("Timer Label");
00688 
00689     // Update parameter in our list and solver timer
00690     if (newLabel != label_) {
00691       label_ = newLabel;
00692     }
00693 
00694 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00695     std::string newSolveLabel = label_ + ": LSQRSolMgr total solve time";
00696     if (timerSolve_.is_null()) {
00697       // Ask TimeMonitor for a new timer.
00698       timerSolve_ = TimeMonitor::getNewCounter (newSolveLabel);
00699     } else {
00700       // We've already created a timer, but we may have changed its
00701       // label.  If we did change its name, then we have to forget
00702       // about the old timer and create a new one with a different
00703       // name.  This is because Teuchos::Time doesn't give you a way
00704       // to change a timer's name, once you've created it.  We assume
00705       // that if the user changed the timer's label, then the user
00706       // wants to reset the timer's results.
00707       const std::string oldSolveLabel = timerSolve_->name ();
00708 
00709       if (oldSolveLabel != newSolveLabel) {
00710         // Tell TimeMonitor to forget about the old timer.
00711         // TimeMonitor lets you clear timers by name.
00712         TimeMonitor::clearCounter (oldSolveLabel);
00713         timerSolve_ = TimeMonitor::getNewCounter (newSolveLabel);
00714       }
00715     }
00716 #endif // BELOS_TEUCHOS_TIME_MONITOR
00717   }
00718 
00719   // Check for a change in verbosity level
00720   {
00721     int newVerbosity = 0;
00722     // ParameterList gets confused sometimes about enums.  This
00723     // ensures that no matter how "Verbosity" was stored -- either an
00724     // an int, or as a Belos::MsgType enum, we will be able to extract
00725     // it.  If it was stored as some other type, we let the exception
00726     // through.
00727     try {
00728       newVerbosity = params->get<Belos::MsgType> ("Verbosity");
00729     } catch (Teuchos::Exceptions::InvalidParameterType&) {
00730       newVerbosity = params->get<int> ("Verbosity");
00731     }
00732     if (newVerbosity != verbosity_) {
00733       verbosity_ = newVerbosity;
00734     }
00735   }
00736 
00737   // (Re)set the output style.
00738   outputStyle_ = params->get<int> ("Output Style");
00739 
00740   // Get the output stream for the output manager.
00741   //
00742   // In case the output stream can't be read back in, we default to
00743   // stdout (std::cout), just to ensure reasonable behavior.
00744   {
00745     outputStream_ = params->get<RCP<std::ostream> > ("Output Stream");
00746 
00747     // We assume that a null output stream indicates that the user
00748     // doesn't want to print anything, so we replace it with a "black
00749     // hole" stream that prints nothing sent to it.  (We can't use a
00750     // null output stream, since the output manager always sends
00751     // things it wants to print to the output stream.)
00752     if (outputStream_.is_null())
00753       outputStream_ = rcp (new Teuchos::oblackholestream);
00754   }
00755 
00756   // Get the frequency of solver output.  (For example, -1 means
00757   // "never," 1 means "every iteration.")
00758   outputFreq_ = params->get<int> ("Output Frequency");
00759 
00760   // Create output manager if we need to, using the verbosity level
00761   // and output stream that we fetched above.  We do this here because
00762   // instantiating an OrthoManager using OrthoManagerFactory requires
00763   // a valid OutputManager.
00764   if (printer_.is_null()) {
00765     printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
00766   } else {
00767     printer_->setVerbosity (verbosity_);
00768     printer_->setOStream (outputStream_);
00769   }
00770 
00771   // Check if the orthogonalization changed, or if we need to
00772   // initialize it.
00773   typedef OrthoManagerFactory<ScalarType, MV, OP> factory_type;
00774   factory_type factory;
00775   bool mustMakeOrtho = false;
00776   {
00777     std::string tempOrthoType;
00778     try {
00779       tempOrthoType = params_->get<std::string> ("Orthogonalization");
00780     } catch (InvalidParameterName&) {
00781       tempOrthoType = orthoType_;
00782     }
00783     if (ortho_.is_null() || tempOrthoType != orthoType_) {
00784       mustMakeOrtho = true;
00785 
00786       // Ensure that the specified orthogonalization type is valid.
00787       if (! factory.isValidName (tempOrthoType)) {
00788         std::ostringstream os;
00789         os << "Belos::LSQRSolMgr: Invalid orthogonalization name \""
00790            << tempOrthoType << "\".  The following are valid options "
00791            << "for the \"Orthogonalization\" name parameter: ";
00792         factory.printValidNames (os);
00793         TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00794       }
00795       orthoType_ = tempOrthoType; // The name is valid, so accept it.
00796       params_->set ("Orthogonalization", orthoType_);
00797     }
00798   }
00799 
00800   // Get any parameters for the orthogonalization ("Orthogonalization
00801   // Parameters").  If not supplied, the orthogonalization manager
00802   // factory will supply default values.
00803   //
00804   // NOTE (mfh 21 Oct 2011) For the sake of backwards compatibility,
00805   // if params has an "Orthogonalization Constant" parameter and the
00806   // DGKS orthogonalization manager is to be used, the value of this
00807   // parameter will override DGKS's "depTol" parameter.
00808   //
00809   // Users must supply the orthogonalization manager parameters as a
00810   // sublist (supplying it as an RCP<ParameterList> would make the
00811   // resulting parameter list not serializable).
00812   RCP<ParameterList> orthoParams;
00813   { // The nonmember function returns an RCP<ParameterList>,
00814     // which is what we want here.
00815     using Teuchos::sublist;
00816     // Abbreviation to avoid typos.
00817     const std::string paramName ("Orthogonalization Parameters");
00818 
00819     try {
00820       orthoParams = sublist (params_, paramName, true);
00821     } catch (InvalidParameter&) {
00822       // We didn't get the parameter list from params, so get a
00823       // default parameter list from the OrthoManagerFactory.
00824       // Modify params_ so that it has the default parameter list,
00825       // and set orthoParams to ensure it's a sublist of params_
00826       // (and not just a copy of one).
00827       params_->set (paramName, factory.getDefaultParameters (orthoType_));
00828       orthoParams = sublist (params_, paramName, true);
00829     }
00830   }
00831   TEUCHOS_TEST_FOR_EXCEPTION(orthoParams.is_null(), std::logic_error,
00832                              "Failed to get orthogonalization parameters.  "
00833                              "Please report this bug to the Belos developers.");
00834 
00835   // If we need to, instantiate a new MatOrthoManager subclass
00836   // instance corresponding to the desired orthogonalization method.
00837   // We've already fetched the orthogonalization method name
00838   // (orthoType_) and its parameters (orthoParams) above.
00839   //
00840   // NOTE (mfh 21 Oct 2011) We only instantiate a new MatOrthoManager
00841   // subclass if the orthogonalization method name is different than
00842   // before.  Thus, for some orthogonalization managers, changes to
00843   // their parameters may not get propagated, if the manager type
00844   // itself didn't change.  The one exception is the "depTol"
00845   // (a.k.a. orthoKappa or "Orthogonalization Constant") parameter of
00846   // DGKS; changes to that _do_ get propagated down to the DGKS
00847   // instance.
00848   //
00849   // The most general way to fix this issue would be to supply each
00850   // orthogonalization manager class with a setParameterList() method
00851   // that takes a parameter list input, and changes the parameters as
00852   // appropriate.  A less efficient but correct way would be simply to
00853   // reinstantiate the OrthoManager every time, whether or not the
00854   // orthogonalization method name or parameters have changed.
00855   if (mustMakeOrtho) {
00856     // Create orthogonalization manager.  This requires that the
00857     // OutputManager (printer_) already be initialized.  LSQR
00858     // currently only orthogonalizes with respect to the Euclidean
00859     // inner product, so we set the inner product matrix M to null.
00860     RCP<const OP> M = null;
00861     ortho_ = factory.makeMatOrthoManager (orthoType_, M, printer_,
00862                                           label_, orthoParams);
00863   }
00864   TEUCHOS_TEST_FOR_EXCEPTION(ortho_.is_null(), std::logic_error,
00865                              "The MatOrthoManager is not yet initialized, but "
00866                              "should be by this point.  "
00867                              "Please report this bug to the Belos developers.");
00868 
00869   // Check which orthogonalization constant to use.  We only need this
00870   // if orthoType_ == "DGKS" (and we already fetched the orthoType_
00871   // parameter above).
00872   if (orthoType_ == "DGKS") {
00873     if (params->isParameter ("Orthogonalization Constant")) {
00874       orthoKappa_ = params_->get<MagnitudeType> ("Orthogonalization Constant");
00875 
00876       if (orthoKappa_ > 0 && ! ortho_.is_null()) {
00877         typedef DGKSOrthoManager<ScalarType,MV,OP> ortho_impl_type;
00878         rcp_dynamic_cast<ortho_impl_type> (ortho_)->setDepTol (orthoKappa_);
00879       }
00880     }
00881   }
00882 
00883   // Check for condition number limit, number of consecutive passed
00884   // iterations, relative RHS error, and relative matrix error.
00885   // Create the LSQR convergence test if necessary.
00886   {
00887     condMax_ = params->get<MagnitudeType> ("Condition Limit");
00888     termIterMax_ = params->get<int>("Term Iter Max");
00889     relRhsErr_ = params->get<MagnitudeType> ("Rel RHS Err");
00890     relMatErr_ = params->get<MagnitudeType> ("Rel Mat Err");
00891 
00892     // Create the LSQR convergence test if it doesn't exist yet.
00893     // Otherwise, update its parameters.
00894     if (convTest_.is_null()) {
00895       convTest_ =
00896         rcp (new LSQRStatusTest<ScalarType,MV,OP> (condMax_, termIterMax_,
00897                                                    relRhsErr_, relMatErr_));
00898     } else {
00899       convTest_->setCondLim (condMax_);
00900       convTest_->setTermIterMax (termIterMax_);
00901       convTest_->setRelRhsErr (relRhsErr_);
00902       convTest_->setRelMatErr (relMatErr_);
00903     }
00904   }
00905 
00906   // Create the status test for maximum number of iterations if
00907   // necessary.  Otherwise, update it with the new maximum iteration
00908   // count.
00909   if (maxIterTest_.is_null()) {
00910     maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
00911   } else {
00912     maxIterTest_->setMaxIters (maxIters_);
00913   }
00914 
00915   // The stopping criterion is an OR combination of the test for
00916   // maximum number of iterations, and the LSQR convergence test.
00917   // ("OR combination" means that both tests will always be evaluated,
00918   // as opposed to a SEQ combination.)
00919   typedef StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00920   // If sTest_ is not null, then maxIterTest_ and convTest_ were
00921   // already constructed on entry to this routine, and sTest_ has
00922   // their pointers.  Thus, maxIterTest_ and convTest_ have gotten any
00923   // parameter changes, so we don't need to do anything to sTest_.
00924   if (sTest_.is_null()) {
00925     sTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::OR,
00926                                          maxIterTest_,
00927                                          convTest_));
00928   }
00929 
00930   if (outputTest_.is_null()) {
00931     // Create the status test output class.
00932     // This class manages and formats the output from the status test.
00933     StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
00934     outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_,
00935                                      Passed + Failed + Undefined);
00936     // Set the solver string for the output test.
00937     const std::string solverDesc = " LSQR ";
00938     outputTest_->setSolverDesc (solverDesc);
00939   } else {
00940     // FIXME (mfh 18 Sep 2011) How do we change the output style of
00941     // outputTest_, without destroying and recreating it?
00942     outputTest_->setOutputManager (printer_);
00943     outputTest_->setChild (sTest_);
00944     outputTest_->setOutputFrequency (outputFreq_);
00945     // Since outputTest_ can only be created here, I'm assuming that
00946     // the fourth constructor argument ("printStates") was set
00947     // correctly on constrution; I don't need to reset it (and I can't
00948     // set it anyway, given StatusTestOutput's interface).
00949   }
00950 
00951   // Inform the solver manager that the current parameters were set.
00952   isSet_ = true;
00953 }
00954 
00955 
00956 template<class ScalarType, class MV, class OP>
00957 Belos::ReturnType LSQRSolMgr<ScalarType,MV,OP,false>::solve() {
00958   using Teuchos::RCP;
00959   using Teuchos::rcp;
00960 
00961   // Set the current parameters if they were not set before.  NOTE:
00962   // This may occur if the user generated the solver manager with the
00963   // default constructor, but did not set any parameters using
00964   // setParameters().
00965   if (!isSet_) {
00966     setParameters (Teuchos::parameterList (*getValidParameters()));
00967   }
00968 
00969   TEUCHOS_TEST_FOR_EXCEPTION(problem_.is_null(), LSQRSolMgrLinearProblemFailure,
00970                      "The linear problem to solve is null.");
00971   TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), LSQRSolMgrLinearProblemFailure,
00972                      "LSQRSolMgr::solve(): The linear problem is not ready, "
00973                      "as its setProblem() method has not been called.");
00974   TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs (*(problem_->getRHS ())) != 1,
00975                      LSQRSolMgrBlockSizeFailure,
00976                      "LSQRSolMgr::solve(): The current implementation of LSQR "
00977                      "only knows how to solve problems with one right-hand "
00978                      "side, but the linear problem to solve has "
00979                      << MVT::GetNumberVecs (*(problem_->getRHS ()))
00980                      << " right-hand sides.");
00981 
00982   // We've validated the LinearProblem instance above.  If any of the
00983   // StatusTests needed to be initialized using information from the
00984   // LinearProblem, now would be the time to do so.  (This is true of
00985   // GMRES, where the residual convergence test(s) to instantiate
00986   // depend on knowing whether there is a left preconditioner.  This
00987   // is why GMRES has an "isSTSet_" Boolean member datum, which tells
00988   // you whether the status tests have been instantiated and are ready
00989   // for use.
00990 
00991   // test isFlexible might go here.
00992 
00993   // Next the right-hand sides to solve are identified.  Among other things,
00994   // this enables getCurrLHSVec() to get the current initial guess vector,
00995   // and getCurrRHSVec() to get the current right-hand side (in Iter).
00996   std::vector<int> currRHSIdx(1, 0);
00997   problem_->setLSIndex(currRHSIdx);
00998 
00999   // Reset the status test.
01000   outputTest_->reset();
01001 
01002   // Don't assume convergence unless we've verified that the
01003   // convergence test passed.
01004   bool isConverged = false;
01005 
01006   // FIXME: Currently we are setting the initial guess to zero, since
01007   // the solver doesn't yet know how to handle a nonzero initial
01008   // guess.  This could be fixed by rewriting the solver to work with
01009   // the residual and a delta.
01010   //
01011   // In a least squares problem with a nonzero initial guess, the
01012   // minimzation problem involves the distance (in a norm depending on
01013   // the preconditioner) between the solution and the the initial
01014   // guess.
01015 
01017   // Solve the linear problem using LSQR
01019 
01020   // Parameter list for the LSQR iteration.
01021   Teuchos::ParameterList plist;
01022 
01023   // Use the solver manager's "Lambda" parameter to set the
01024   // iteration's "Lambda" parameter.  We know that the solver
01025   // manager's parameter list (params_) does indeed contain the
01026   // "Lambda" parameter, because solve() always ensures that
01027   // setParameters() has been called.
01028   plist.set ("Lambda", lambda_);
01029 
01030   typedef LSQRIter<ScalarType,MV,OP> iter_type;
01031   RCP<iter_type> lsqr_iter =
01032     rcp (new iter_type (problem_, printer_, outputTest_, plist));
01033 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01034   Teuchos::TimeMonitor slvtimer(*timerSolve_);
01035 #endif
01036 
01037   // Reset the number of iterations.
01038   lsqr_iter->resetNumIters();
01039   // Reset the number of calls that the status test output knows about.
01040   outputTest_->resetNumCalls();
01041   // Set the new state and initialize the solver.
01042   LSQRIterationState<ScalarType,MV> newstate;
01043   lsqr_iter->initializeLSQR(newstate);
01044   // tell lsqr_iter to iterate
01045   try {
01046     lsqr_iter->iterate();
01047 
01048     // First check for convergence.  If we didn't converge, then check
01049     // whether we reached the maximum number of iterations.  If
01050     // neither of those happened, there must have been a bug.
01051     if (convTest_->getStatus() == Belos::Passed) {
01052       isConverged = true;
01053     } else if (maxIterTest_->getStatus() == Belos::Passed) {
01054       isConverged = false;
01055     } else {
01056       TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
01057                          "LSQRSolMgr::solve(): LSQRIteration::iterate() "
01058                          "returned without either the convergence test or "
01059                          "the maximum iteration count test passing."
01060                          "Please report this bug to the Belos developers.");
01061     }
01062   } catch (const std::exception &e) {
01063     printer_->stream(Belos::Errors) << "Error! Caught std::exception in LSQRIter::iterate() at iteration "
01064                                     << lsqr_iter->getNumIters() << std::endl
01065                                     << e.what() << std::endl;
01066     throw;
01067   }
01068 
01069   // identify current linear system as solved LinearProblem
01070   problem_->setCurrLS();
01071   // print final summary
01072   sTest_->print (printer_->stream (Belos::FinalSummary));
01073 
01074   // Print timing information, if the corresponding compile-time and
01075   // run-time options are enabled.
01076 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01077   // Calling summarize() can be expensive, so don't call unless the
01078   // user wants to print out timing details.  summarize() will do all
01079   // the work even if it's passed a "black hole" output stream.
01080   if (verbosity_ & TimingDetails)
01081     Teuchos::TimeMonitor::summarize (printer_->stream (Belos::TimingDetails));
01082 #endif // BELOS_TEUCHOS_TIME_MONITOR
01083 
01084   // A posteriori solve information
01085   numIters_ = maxIterTest_->getNumIters();
01086   matCondNum_ = convTest_->getMatCondNum();
01087   matNorm_ = convTest_->getMatNorm();
01088   resNorm_ = convTest_->getResidNorm();
01089   matResNorm_ = convTest_->getLSResidNorm();
01090 
01091   if (! isConverged) {
01092     return Belos::Unconverged;
01093   } else {
01094     return Belos::Converged;
01095   }
01096 }
01097 
01098 // LSQRSolMgr requires the solver manager to return an eponymous std::string.
01099 template<class ScalarType, class MV, class OP>
01100 std::string LSQRSolMgr<ScalarType,MV,OP,false>::description () const
01101 {
01102   std::ostringstream oss;
01103   oss << "LSQRSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
01104   oss << "{";
01105   oss << "Ortho Type='"<<orthoType_<<"'";
01106   oss << ", Lambda="<< lambda_;
01107   oss << ", condition number limit="<< condMax_;
01108   oss << ", relative RHS Error="<< relRhsErr_;
01109   oss << ", relative Matrix Error="<< relMatErr_;
01110   oss << ", maximum number of iterations="<< maxIters_;
01111   oss << ", termIterMax="<<termIterMax_;
01112   oss << "}";
01113   return oss.str();
01114 }
01115 
01116 } // end Belos namespace
01117 
01118 #endif /* BELOS_LSQR_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines