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   isSet_(false),
00506   loaDetected_(false)
00507 {}
00508 
00509 
00510 // Basic Constructor
00511 template<class ScalarType, class MV, class OP>
00512 LSQRSolMgr<ScalarType,MV,OP,false>::
00513 LSQRSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00514             const Teuchos::RCP<Teuchos::ParameterList> &pl) :
00515   problem_(problem),
00516   isSet_(false),
00517   loaDetected_(false)
00518 {
00519   // The linear problem to solve is allowed to be null here.  The user
00520   // must then set a nonnull linear problem (by calling setProblem())
00521   // before calling solve().
00522   //
00523   // Similarly, users are allowed to set a null parameter list here,
00524   // but they must first set a nonnull parameter list (by calling
00525   // setParameters()) before calling solve().
00526   if (! is_null (pl)) {
00527     setParameters (pl);
00528   }
00529 }
00530 
00531 
00532 template<class ScalarType, class MV, class OP>
00533 Teuchos::RCP<const Teuchos::ParameterList>
00534 LSQRSolMgr<ScalarType,MV,OP,false>::getValidParameters() const
00535 {
00536   using Teuchos::ParameterList;
00537   using Teuchos::parameterList;
00538   using Teuchos::RCP;
00539   using Teuchos::rcp;
00540   using Teuchos::rcpFromRef;
00541   typedef Teuchos::ScalarTraits<MagnitudeType> STM;
00542 
00543   // Set all the valid parameters and their default values.
00544   if (is_null (validParams_)) {
00545     // We use Teuchos::as just in case MagnitudeType doesn't have a
00546     // constructor that takes an int.  Otherwise, we could just write
00547     // "MagnitudeType(10)".
00548     const MagnitudeType ten = Teuchos::as<MagnitudeType> (10);
00549     const MagnitudeType sqrtEps = STM::squareroot (STM::eps());
00550 
00551     const MagnitudeType lambda = STM::zero();
00552     RCP<std::ostream> outputStream = rcpFromRef (std::cout);
00553     const MagnitudeType relRhsErr = ten * sqrtEps;
00554     const MagnitudeType relMatErr = ten * sqrtEps;
00555     const MagnitudeType condMax = STM::one() / STM::eps();
00556     const int maxIters = 1000;
00557     const int termIterMax = 1;
00558     const std::string orthoType ("DGKS");
00559     const MagnitudeType orthoKappa = Teuchos::as<MagnitudeType> (-1.0);
00560     const int verbosity = Belos::Errors;
00561     const int outputStyle = Belos::General;
00562     const int outputFreq = -1;
00563     const std::string label ("Belos");
00564 
00565     RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00566     pl->set("Output Stream", outputStream,
00567       "is a reference-counted pointer to the output stream receiving\n"
00568       "all solver output.");
00569     pl->set("Lambda", lambda, "is the damping parameter.");
00570     pl->set("Rel RHS Err", relRhsErr,
00571             "estimates the error in the data defining the right-\n"
00572             "hand side.");
00573     pl->set("Rel Mat Err", relMatErr,
00574             "estimates the error in the data defining the matrix.");
00575     pl->set("Condition Limit", condMax,
00576       "bounds the estimated condition number of Abar.");
00577     pl->set("Maximum Iterations", maxIters,
00578       "allows at most the maximum number of iterations.");
00579     pl->set("Term Iter Max", termIterMax,
00580       "consecutive iterations meeting thresholds are necessary for\n"
00581        "for convergence.");
00582     pl->set("Orthogonalization", orthoType,
00583       "uses orthogonalization of either DGKS, ICGS, IMGS, or TSQR.");
00584     {
00585       OrthoManagerFactory<ScalarType, MV, OP> factory;
00586       pl->set("Orthogonalization", orthoType,
00587               "refers to the orthogonalization method to use.  Valid "
00588               "options: " + factory.validNamesString());
00589       RCP<const ParameterList> orthoParams =
00590         factory.getDefaultParameters (orthoType);
00591       pl->set ("Orthogonalization Parameters", *orthoParams,
00592                "Parameters specific to the type of orthogonalization used.");
00593     }
00594     pl->set("Orthogonalization Constant", orthoKappa,
00595       "is the threshold used by DGKS orthogonalization to determine\n"
00596       "whether or not to repeat classical Gram-Schmidt.  This parameter\n"
00597       "is ignored if \"Orthogonalization\" is not \"DGKS\".");
00598     pl->set("Verbosity", verbosity,
00599       "type(s) of solver information are outputted to the output\n"
00600       "stream.");
00601     pl->set("Output Style", outputStyle,
00602       "the style used for the solver information outputted to the\n"
00603       "output stream.");
00604     pl->set("Output Frequency", outputFreq,
00605       "is the frequency at which information is written to the\n"
00606       "output stream.");
00607     pl->set("Timer Label", label,
00608       "is the string to use as a prefix for the timer labels.");
00609     //  pl->set("Restart Timers", restartTimers_);
00610     pl->set("Block Size", 1, "Block size parameter (currently, this must always be 1).");
00611 
00612     validParams_ = pl;
00613   }
00614   return validParams_;
00615 }
00616 
00617 
00618 template<class ScalarType, class MV, class OP>
00619 void
00620 LSQRSolMgr<ScalarType,MV,OP,false>::
00621 setParameters (const Teuchos::RCP<Teuchos::ParameterList> &params)
00622 {
00623   using Teuchos::isParameterType;
00624   using Teuchos::getParameter;
00625   using Teuchos::null;
00626   using Teuchos::ParameterList;
00627   using Teuchos::parameterList;
00628   using Teuchos::RCP;
00629   using Teuchos::rcp;
00630   using Teuchos::rcp_dynamic_cast;
00631   using Teuchos::rcpFromRef;
00632   using Teuchos::Time;
00633   using Teuchos::TimeMonitor;
00634   using Teuchos::Exceptions::InvalidParameter;
00635   using Teuchos::Exceptions::InvalidParameterName;
00636   using Teuchos::Exceptions::InvalidParameterType;
00637 
00638   TEUCHOS_TEST_FOR_EXCEPTION(params.is_null(), std::invalid_argument,
00639                      "Belos::LSQRSolMgr::setParameters: "
00640                      "the input ParameterList is null.");
00641   RCP<const ParameterList> defaultParams = getValidParameters ();
00642   params->validateParametersAndSetDefaults (*defaultParams);
00643 
00644   // At this point, params is a valid parameter list with defaults
00645   // filled in.  Now we can "commit" it to our instance's parameter
00646   // list.
00647   params_ = params;
00648 
00649   // Get the damping (a.k.a. regularization) parameter lambda.
00650   lambda_ = params->get<MagnitudeType> ("Lambda");
00651 
00652   // Get the maximum number of iterations.
00653   maxIters_ = params->get<int>("Maximum Iterations");
00654 
00655   // (Re)set the timer label.
00656   {
00657     const std::string newLabel = params->get<std::string>("Timer Label");
00658 
00659     // Update parameter in our list and solver timer
00660     if (newLabel != label_) {
00661       label_ = newLabel;
00662     }
00663 
00664 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00665     std::string newSolveLabel = label_ + ": LSQRSolMgr total solve time";
00666     if (timerSolve_.is_null()) {
00667       // Ask TimeMonitor for a new timer.
00668       timerSolve_ = TimeMonitor::getNewCounter (newSolveLabel);
00669     } else {
00670       // We've already created a timer, but we may have changed its
00671       // label.  If we did change its name, then we have to forget
00672       // about the old timer and create a new one with a different
00673       // name.  This is because Teuchos::Time doesn't give you a way
00674       // to change a timer's name, once you've created it.  We assume
00675       // that if the user changed the timer's label, then the user
00676       // wants to reset the timer's results.
00677       const std::string oldSolveLabel = timerSolve_->name ();
00678 
00679       if (oldSolveLabel != newSolveLabel) {
00680         // Tell TimeMonitor to forget about the old timer.
00681         // TimeMonitor lets you clear timers by name.
00682         TimeMonitor::clearCounter (oldSolveLabel);
00683         timerSolve_ = TimeMonitor::getNewCounter (newSolveLabel);
00684       }
00685     }
00686 #endif // BELOS_TEUCHOS_TIME_MONITOR
00687   }
00688 
00689   // Check for a change in verbosity level
00690   {
00691     int newVerbosity = 0;
00692     // ParameterList gets confused sometimes about enums.  This
00693     // ensures that no matter how "Verbosity" was stored -- either an
00694     // an int, or as a Belos::MsgType enum, we will be able to extract
00695     // it.  If it was stored as some other type, we let the exception
00696     // through.
00697     try {
00698       newVerbosity = params->get<Belos::MsgType> ("Verbosity");
00699     } catch (Teuchos::Exceptions::InvalidParameterType&) {
00700       newVerbosity = params->get<int> ("Verbosity");
00701     }
00702     if (newVerbosity != verbosity_) {
00703       verbosity_ = newVerbosity;
00704     }
00705   }
00706 
00707   // (Re)set the output style.
00708   outputStyle_ = params->get<int> ("Output Style");
00709 
00710   // Get the output stream for the output manager.
00711   //
00712   // In case the output stream can't be read back in, we default to
00713   // stdout (std::cout), just to ensure reasonable behavior.
00714   {
00715     outputStream_ = params->get<RCP<std::ostream> > ("Output Stream");
00716 
00717     // We assume that a null output stream indicates that the user
00718     // doesn't want to print anything, so we replace it with a "black
00719     // hole" stream that prints nothing sent to it.  (We can't use a
00720     // null output stream, since the output manager always sends
00721     // things it wants to print to the output stream.)
00722     if (outputStream_.is_null())
00723       outputStream_ = rcp (new Teuchos::oblackholestream);
00724   }
00725 
00726   // Get the frequency of solver output.  (For example, -1 means
00727   // "never," 1 means "every iteration.")
00728   outputFreq_ = params->get<int> ("Output Frequency");
00729 
00730   // Create output manager if we need to, using the verbosity level
00731   // and output stream that we fetched above.  We do this here because
00732   // instantiating an OrthoManager using OrthoManagerFactory requires
00733   // a valid OutputManager.
00734   if (printer_.is_null()) {
00735     printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
00736   } else {
00737     printer_->setVerbosity (verbosity_);
00738     printer_->setOStream (outputStream_);
00739   }
00740 
00741   // Check if the orthogonalization changed, or if we need to
00742   // initialize it.
00743   typedef OrthoManagerFactory<ScalarType, MV, OP> factory_type;
00744   factory_type factory;
00745   bool mustMakeOrtho = false;
00746   {
00747     std::string tempOrthoType;
00748     try {
00749       tempOrthoType = params_->get<std::string> ("Orthogonalization");
00750     } catch (InvalidParameterName&) {
00751       tempOrthoType = orthoType_;
00752     }
00753     if (ortho_.is_null() || tempOrthoType != orthoType_) {
00754       mustMakeOrtho = true;
00755 
00756       // Ensure that the specified orthogonalization type is valid.
00757       if (! factory.isValidName (tempOrthoType)) {
00758         std::ostringstream os;
00759         os << "Belos::LSQRSolMgr: Invalid orthogonalization name \""
00760            << tempOrthoType << "\".  The following are valid options "
00761            << "for the \"Orthogonalization\" name parameter: ";
00762         factory.printValidNames (os);
00763         TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00764       }
00765       orthoType_ = tempOrthoType; // The name is valid, so accept it.
00766       params_->set ("Orthogonalization", orthoType_);
00767     }
00768   }
00769 
00770   // Get any parameters for the orthogonalization ("Orthogonalization
00771   // Parameters").  If not supplied, the orthogonalization manager
00772   // factory will supply default values.
00773   //
00774   // NOTE (mfh 21 Oct 2011) For the sake of backwards compatibility,
00775   // if params has an "Orthogonalization Constant" parameter and the
00776   // DGKS orthogonalization manager is to be used, the value of this
00777   // parameter will override DGKS's "depTol" parameter.
00778   //
00779   // Users must supply the orthogonalization manager parameters as a
00780   // sublist (supplying it as an RCP<ParameterList> would make the
00781   // resulting parameter list not serializable).
00782   RCP<ParameterList> orthoParams;
00783   { // The nonmember function returns an RCP<ParameterList>,
00784     // which is what we want here.
00785     using Teuchos::sublist;
00786     // Abbreviation to avoid typos.
00787     const std::string paramName ("Orthogonalization Parameters");
00788 
00789     try {
00790       orthoParams = sublist (params_, paramName, true);
00791     } catch (InvalidParameter&) {
00792       // We didn't get the parameter list from params, so get a
00793       // default parameter list from the OrthoManagerFactory.
00794       // Modify params_ so that it has the default parameter list,
00795       // and set orthoParams to ensure it's a sublist of params_
00796       // (and not just a copy of one).
00797       params_->set (paramName, factory.getDefaultParameters (orthoType_));
00798       orthoParams = sublist (params_, paramName, true);
00799     }
00800   }
00801   TEUCHOS_TEST_FOR_EXCEPTION(orthoParams.is_null(), std::logic_error,
00802                              "Failed to get orthogonalization parameters.  "
00803                              "Please report this bug to the Belos developers.");
00804 
00805   // If we need to, instantiate a new MatOrthoManager subclass
00806   // instance corresponding to the desired orthogonalization method.
00807   // We've already fetched the orthogonalization method name
00808   // (orthoType_) and its parameters (orthoParams) above.
00809   //
00810   // NOTE (mfh 21 Oct 2011) We only instantiate a new MatOrthoManager
00811   // subclass if the orthogonalization method name is different than
00812   // before.  Thus, for some orthogonalization managers, changes to
00813   // their parameters may not get propagated, if the manager type
00814   // itself didn't change.  The one exception is the "depTol"
00815   // (a.k.a. orthoKappa or "Orthogonalization Constant") parameter of
00816   // DGKS; changes to that _do_ get propagated down to the DGKS
00817   // instance.
00818   //
00819   // The most general way to fix this issue would be to supply each
00820   // orthogonalization manager class with a setParameterList() method
00821   // that takes a parameter list input, and changes the parameters as
00822   // appropriate.  A less efficient but correct way would be simply to
00823   // reinstantiate the OrthoManager every time, whether or not the
00824   // orthogonalization method name or parameters have changed.
00825   if (mustMakeOrtho) {
00826     // Create orthogonalization manager.  This requires that the
00827     // OutputManager (printer_) already be initialized.  LSQR
00828     // currently only orthogonalizes with respect to the Euclidean
00829     // inner product, so we set the inner product matrix M to null.
00830     RCP<const OP> M = null;
00831     ortho_ = factory.makeMatOrthoManager (orthoType_, M, printer_,
00832                                           label_, orthoParams);
00833   }
00834   TEUCHOS_TEST_FOR_EXCEPTION(ortho_.is_null(), std::logic_error,
00835                              "The MatOrthoManager is not yet initialized, but "
00836                              "should be by this point.  "
00837                              "Please report this bug to the Belos developers.");
00838 
00839   // Check which orthogonalization constant to use.  We only need this
00840   // if orthoType_ == "DGKS" (and we already fetched the orthoType_
00841   // parameter above).
00842   if (orthoType_ == "DGKS") {
00843     if (params->isParameter ("Orthogonalization Constant")) {
00844       orthoKappa_ = params_->get<MagnitudeType> ("Orthogonalization Constant");
00845 
00846       if (orthoKappa_ > 0 && ! ortho_.is_null()) {
00847         typedef DGKSOrthoManager<ScalarType,MV,OP> ortho_impl_type;
00848         rcp_dynamic_cast<ortho_impl_type> (ortho_)->setDepTol (orthoKappa_);
00849       }
00850     }
00851   }
00852 
00853   // Check for condition number limit, number of consecutive passed
00854   // iterations, relative RHS error, and relative matrix error.
00855   // Create the LSQR convergence test if necessary.
00856   {
00857     condMax_ = params->get<MagnitudeType> ("Condition Limit");
00858     termIterMax_ = params->get<int>("Term Iter Max");
00859     relRhsErr_ = params->get<MagnitudeType> ("Rel RHS Err");
00860     relMatErr_ = params->get<MagnitudeType> ("Rel Mat Err");
00861 
00862     // Create the LSQR convergence test if it doesn't exist yet.
00863     // Otherwise, update its parameters.
00864     if (convTest_.is_null()) {
00865       convTest_ =
00866         rcp (new LSQRStatusTest<ScalarType,MV,OP> (condMax_, termIterMax_,
00867                                                    relRhsErr_, relMatErr_));
00868     } else {
00869       convTest_->setCondLim (condMax_);
00870       convTest_->setTermIterMax (termIterMax_);
00871       convTest_->setRelRhsErr (relRhsErr_);
00872       convTest_->setRelMatErr (relMatErr_);
00873     }
00874   }
00875 
00876   // Create the status test for maximum number of iterations if
00877   // necessary.  Otherwise, update it with the new maximum iteration
00878   // count.
00879   if (maxIterTest_.is_null()) {
00880     maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
00881   } else {
00882     maxIterTest_->setMaxIters (maxIters_);
00883   }
00884 
00885   // The stopping criterion is an OR combination of the test for
00886   // maximum number of iterations, and the LSQR convergence test.
00887   // ("OR combination" means that both tests will always be evaluated,
00888   // as opposed to a SEQ combination.)
00889   typedef StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00890   // If sTest_ is not null, then maxIterTest_ and convTest_ were
00891   // already constructed on entry to this routine, and sTest_ has
00892   // their pointers.  Thus, maxIterTest_ and convTest_ have gotten any
00893   // parameter changes, so we don't need to do anything to sTest_.
00894   if (sTest_.is_null()) {
00895     sTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::OR,
00896                                          maxIterTest_,
00897                                          convTest_));
00898   }
00899 
00900   if (outputTest_.is_null()) {
00901     // Create the status test output class.
00902     // This class manages and formats the output from the status test.
00903     StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
00904     outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_,
00905                                      Passed + Failed + Undefined);
00906     // Set the solver string for the output test.
00907     const std::string solverDesc = " LSQR ";
00908     outputTest_->setSolverDesc (solverDesc);
00909   } else {
00910     // FIXME (mfh 18 Sep 2011) How do we change the output style of
00911     // outputTest_, without destroying and recreating it?
00912     outputTest_->setOutputManager (printer_);
00913     outputTest_->setChild (sTest_);
00914     outputTest_->setOutputFrequency (outputFreq_);
00915     // Since outputTest_ can only be created here, I'm assuming that
00916     // the fourth constructor argument ("printStates") was set
00917     // correctly on constrution; I don't need to reset it (and I can't
00918     // set it anyway, given StatusTestOutput's interface).
00919   }
00920 
00921   // Inform the solver manager that the current parameters were set.
00922   isSet_ = true;
00923 }
00924 
00925 
00926 template<class ScalarType, class MV, class OP>
00927 Belos::ReturnType LSQRSolMgr<ScalarType,MV,OP,false>::solve() {
00928   using Teuchos::RCP;
00929   using Teuchos::rcp;
00930 
00931   // Set the current parameters if they were not set before.  NOTE:
00932   // This may occur if the user generated the solver manager with the
00933   // default constructor, but did not set any parameters using
00934   // setParameters().
00935   if (!isSet_) {
00936     setParameters (Teuchos::parameterList (*getValidParameters()));
00937   }
00938 
00939   TEUCHOS_TEST_FOR_EXCEPTION(problem_.is_null(), LSQRSolMgrLinearProblemFailure,
00940                      "The linear problem to solve is null.");
00941   TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), LSQRSolMgrLinearProblemFailure,
00942                      "LSQRSolMgr::solve(): The linear problem is not ready, "
00943                      "as its setProblem() method has not been called.");
00944   TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs (*(problem_->getRHS ())) != 1,
00945                      LSQRSolMgrBlockSizeFailure,
00946                      "LSQRSolMgr::solve(): The current implementation of LSQR "
00947                      "only knows how to solve problems with one right-hand "
00948                      "side, but the linear problem to solve has "
00949                      << MVT::GetNumberVecs (*(problem_->getRHS ()))
00950                      << " right-hand sides.");
00951 
00952   // We've validated the LinearProblem instance above.  If any of the
00953   // StatusTests needed to be initialized using information from the
00954   // LinearProblem, now would be the time to do so.  (This is true of
00955   // GMRES, where the residual convergence test(s) to instantiate
00956   // depend on knowing whether there is a left preconditioner.  This
00957   // is why GMRES has an "isSTSet_" Boolean member datum, which tells
00958   // you whether the status tests have been instantiated and are ready
00959   // for use.
00960 
00961   // test isFlexible might go here.
00962 
00963   // Next the right-hand sides to solve are identified.  Among other things,
00964   // this enables getCurrLHSVec() to get the current initial guess vector,
00965   // and getCurrRHSVec() to get the current right-hand side (in Iter).
00966   std::vector<int> currRHSIdx(1, 0);
00967   problem_->setLSIndex(currRHSIdx);
00968 
00969   // Reset the status test.
00970   outputTest_->reset();
00971 
00972   // Don't assume convergence unless we've verified that the
00973   // convergence test passed.
00974   bool isConverged = false;
00975 
00976   // FIXME: Currently we are setting the initial guess to zero, since
00977   // the solver doesn't yet know how to handle a nonzero initial
00978   // guess.  This could be fixed by rewriting the solver to work with
00979   // the residual and a delta.
00980   //
00981   // In a least squares problem with a nonzero initial guess, the
00982   // minimzation problem involves the distance (in a norm depending on
00983   // the preconditioner) between the solution and the the initial
00984   // guess.
00985 
00987   // Solve the linear problem using LSQR
00989 
00990   // Parameter list for the LSQR iteration.
00991   Teuchos::ParameterList plist;
00992 
00993   // Use the solver manager's "Lambda" parameter to set the
00994   // iteration's "Lambda" parameter.  We know that the solver
00995   // manager's parameter list (params_) does indeed contain the
00996   // "Lambda" parameter, because solve() always ensures that
00997   // setParameters() has been called.
00998   plist.set ("Lambda", lambda_);
00999 
01000   typedef LSQRIter<ScalarType,MV,OP> iter_type;
01001   RCP<iter_type> lsqr_iter =
01002     rcp (new iter_type (problem_, printer_, outputTest_, plist));
01003 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01004   Teuchos::TimeMonitor slvtimer(*timerSolve_);
01005 #endif
01006 
01007   // Reset the number of iterations.
01008   lsqr_iter->resetNumIters();
01009   // Reset the number of calls that the status test output knows about.
01010   outputTest_->resetNumCalls();
01011   // Set the new state and initialize the solver.
01012   LSQRIterationState<ScalarType,MV> newstate;
01013   lsqr_iter->initializeLSQR(newstate);
01014   // tell lsqr_iter to iterate
01015   try {
01016     lsqr_iter->iterate();
01017 
01018     // First check for convergence.  If we didn't converge, then check
01019     // whether we reached the maximum number of iterations.  If
01020     // neither of those happened, there must have been a bug.
01021     if (convTest_->getStatus() == Belos::Passed) {
01022       isConverged = true;
01023     } else if (maxIterTest_->getStatus() == Belos::Passed) {
01024       isConverged = false;
01025     } else {
01026       TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
01027                          "LSQRSolMgr::solve(): LSQRIteration::iterate() "
01028                          "returned without either the convergence test or "
01029                          "the maximum iteration count test passing."
01030                          "Please report this bug to the Belos developers.");
01031     }
01032   } catch (const std::exception &e) {
01033     printer_->stream(Belos::Errors) << "Error! Caught std::exception in LSQRIter::iterate() at iteration "
01034                                     << lsqr_iter->getNumIters() << std::endl
01035                                     << e.what() << std::endl;
01036     throw;
01037   }
01038 
01039   // identify current linear system as solved LinearProblem
01040   problem_->setCurrLS();
01041   // print final summary
01042   sTest_->print (printer_->stream (Belos::FinalSummary));
01043 
01044   // Print timing information, if the corresponding compile-time and
01045   // run-time options are enabled.
01046 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01047   // Calling summarize() can be expensive, so don't call unless the
01048   // user wants to print out timing details.  summarize() will do all
01049   // the work even if it's passed a "black hole" output stream.
01050   if (verbosity_ & TimingDetails)
01051     Teuchos::TimeMonitor::summarize (printer_->stream (Belos::TimingDetails));
01052 #endif // BELOS_TEUCHOS_TIME_MONITOR
01053 
01054   // A posteriori solve information
01055   numIters_ = maxIterTest_->getNumIters();
01056   matCondNum_ = convTest_->getMatCondNum();
01057   matNorm_ = convTest_->getMatNorm();
01058   resNorm_ = convTest_->getResidNorm();
01059   matResNorm_ = convTest_->getLSResidNorm();
01060 
01061   if (! isConverged) {
01062     return Belos::Unconverged;
01063   } else {
01064     return Belos::Converged;
01065   }
01066 }
01067 
01068 // LSQRSolMgr requires the solver manager to return an eponymous std::string.
01069 template<class ScalarType, class MV, class OP>
01070 std::string LSQRSolMgr<ScalarType,MV,OP,false>::description () const
01071 {
01072   std::ostringstream oss;
01073   oss << "LSQRSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
01074   oss << "{";
01075   oss << "Ortho Type='"<<orthoType_<<"'";
01076   oss << ", Lambda="<< lambda_;
01077   oss << ", condition number limit="<< condMax_;
01078   oss << ", relative RHS Error="<< relRhsErr_;
01079   oss << ", relative Matrix Error="<< relMatErr_;
01080   oss << ", maximum number of iterations="<< maxIters_;
01081   oss << ", termIterMax="<<termIterMax_;
01082   oss << "}";
01083   return oss.str();
01084 }
01085 
01086 } // end Belos namespace
01087 
01088 #endif /* BELOS_LSQR_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines