Belos Package Browser (Single Doxygen Collection) Development
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 
00049 #include "BelosConfigDefs.hpp"
00050 #include "BelosTypes.hpp"
00051 
00052 #include "BelosLinearProblem.hpp"
00053 #include "BelosSolverManager.hpp"
00054 
00055 #include "BelosLSQRIteration.hpp"
00056 #include "BelosLSQRIter.hpp"
00057 #include "BelosOrthoManagerFactory.hpp"
00058 #include "BelosStatusTestMaxIters.hpp"
00059 #include "BelosLSQRStatusTest.hpp"
00060 #include "BelosStatusTestCombo.hpp"
00061 #include "BelosStatusTestOutputFactory.hpp"
00062 #include "BelosOutputManager.hpp"
00063 #include "Teuchos_as.hpp"
00064 #include "Teuchos_BLAS.hpp"
00065 #include "Teuchos_LAPACK.hpp"
00066 
00067 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00068 #include "Teuchos_TimeMonitor.hpp"
00069 #endif
00070 
00081 namespace Belos {
00082 
00083   
00085 
00086 
00093 class LSQRSolMgrLinearProblemFailure : public BelosError {
00094 public:
00095   LSQRSolMgrLinearProblemFailure(const std::string& what_arg)
00096     : BelosError(what_arg)
00097   {}
00098 };
00099 
00106 class LSQRSolMgrOrthoFailure : public BelosError {
00107 public:
00108   LSQRSolMgrOrthoFailure(const std::string& what_arg)
00109     : BelosError(what_arg)
00110   {}
00111 };
00112 
00119 class LSQRSolMgrBlockSizeFailure : public BelosError {
00120 public:
00121   LSQRSolMgrBlockSizeFailure(const std::string& what_arg) 
00122   : BelosError(what_arg)
00123   {}
00124 };
00125 
00170 template<class ScalarType, class MV, class OP>
00171 class LSQRSolMgr : public SolverManager<ScalarType,MV,OP> {
00172   
00173 private:
00174   typedef MultiVecTraits<ScalarType,MV> MVT;
00175   typedef OperatorTraits<ScalarType,MV,OP> OPT;
00176   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00177   typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00178   typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00179   
00180 public:
00181   
00183 
00184   
00191   LSQRSolMgr();
00192   
00215   // This LSQR implementation only supports block size 1.  Like CG, LSQR is a short
00216   // recurrence method that, in finite precision arithmetic and without reorthogonalization, 
00217   // does not have the "n" step convergence property.  Without either blocks or 
00218   // reorthogonalization, there is nothing to "Orthogonalize."
00219 
00220   LSQRSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00221      const Teuchos::RCP<Teuchos::ParameterList> &pl );
00222   
00224   virtual ~LSQRSolMgr() {};
00226   
00228 
00229 
00232   const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00233     return *problem_;
00234   }
00235   
00238   Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00239   
00242   Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
00243   
00249   Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00250     return Teuchos::tuple(timerSolve_);
00251   }
00252   
00254   int getNumIters() const {
00255     return numIters_;
00256   }
00257 
00262   MagnitudeType getMatCondNum () const {
00263     return matCondNum_; 
00264   }
00265 
00270   MagnitudeType getMatNorm () const {
00271     return matNorm_;
00272   }
00273 
00282   MagnitudeType getResNorm () const {
00283     return resNorm_;
00284   }
00285 
00287   MagnitudeType getMatResNorm () const {
00288     return matResNorm_;
00289   }
00290   
00299   bool isLOADetected() const { return false; }
00300   
00302   
00304 
00305    
00307   void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; }
00308   
00310   void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
00311     
00313    
00315 
00316 
00320   void reset( const ResetType type ) { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
00322  
00324 
00325     
00344   ReturnType solve();
00345     
00347     
00350     
00352   std::string description() const;
00353     
00355     
00356 private:
00357 
00358   // Linear problem.
00359   Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00360   
00361   // Output manager.
00362   Teuchos::RCP<OutputManager<ScalarType> > printer_;
00363   Teuchos::RCP<std::ostream> outputStream_;
00364   
00365   // Status test.
00366   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00367   Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00368   Teuchos::RCP<LSQRStatusTest<ScalarType,MV,OP> > convTest_;
00369   Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00370   
00371   // Orthogonalization manager.
00372   Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_; 
00373     
00374   // Current parameter list.
00375   Teuchos::RCP<Teuchos::ParameterList> params_;
00381   mutable Teuchos::RCP<const Teuchos::ParameterList> validParams_;
00382 
00383   // Current solver input parameters
00384   MagnitudeType lambda_;
00385   MagnitudeType relRhsErr_;
00386   MagnitudeType relMatErr_;
00387   MagnitudeType condMax_;
00388   int maxIters_, termIterMax_;
00389   std::string orthoType_; 
00390   MagnitudeType orthoKappa_;
00391   int verbosity_, outputStyle_, outputFreq_;
00392 
00393   // Terminal solver state values
00394   int numIters_;
00395   MagnitudeType matCondNum_;
00396   MagnitudeType matNorm_;
00397   MagnitudeType resNorm_;
00398   MagnitudeType matResNorm_;
00399     
00400   // Timers.
00401   std::string label_;
00402   Teuchos::RCP<Teuchos::Time> timerSolve_;
00403 
00404   // Internal state variables.
00405   bool isSet_;
00406   bool loaDetected_;
00407 };
00408 
00409 // Empty Constructor
00410 template<class ScalarType, class MV, class OP>
00411 LSQRSolMgr<ScalarType,MV,OP>::LSQRSolMgr() :
00412   isSet_(false),
00413   loaDetected_(false)
00414 {}
00415 
00416 
00417 // Basic Constructor
00418 template<class ScalarType, class MV, class OP>
00419 LSQRSolMgr<ScalarType,MV,OP>::
00420 LSQRSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00421       const Teuchos::RCP<Teuchos::ParameterList> &pl) : 
00422   problem_(problem),
00423   isSet_(false),
00424   loaDetected_(false)
00425 {
00426   // The linear problem to solve is allowed to be null here.  The user
00427   // must then set a nonnull linear problem (by calling setProblem())
00428   // before calling solve().
00429   // 
00430   // Similarly, users are allowed to set a null parameter list here,
00431   // but they must first set a nonnull parameter list (by calling
00432   // setParameters()) before calling solve().
00433   if (! is_null (pl)) {
00434     setParameters (pl);  
00435   }
00436 }
00437 
00438 
00439 template<class ScalarType, class MV, class OP>
00440 Teuchos::RCP<const Teuchos::ParameterList> 
00441 LSQRSolMgr<ScalarType,MV,OP>::getValidParameters() const
00442 {
00443   using Teuchos::ParameterList;
00444   using Teuchos::parameterList;
00445   using Teuchos::RCP;
00446   using Teuchos::rcp;
00447   using Teuchos::rcpFromRef;
00448   typedef Teuchos::ScalarTraits<MagnitudeType> STM;
00449 
00450   // Set all the valid parameters and their default values.
00451   if (is_null (validParams_)) {
00452     // We use Teuchos::as just in case MagnitudeType doesn't have a
00453     // constructor that takes an int.  Otherwise, we could just write
00454     // "MagnitudeType(10)".
00455     const MagnitudeType ten = Teuchos::as<MagnitudeType> (10);
00456     const MagnitudeType sqrtEps = STM::squareroot (STM::eps());
00457 
00458     const MagnitudeType lambda = STM::zero();
00459     RCP<std::ostream> outputStream = rcpFromRef (std::cout);
00460     const MagnitudeType relRhsErr = ten * sqrtEps;
00461     const MagnitudeType relMatErr = ten * sqrtEps;
00462     const MagnitudeType condMax = STM::one() / STM::eps();
00463     const int maxIters = 1000;
00464     const int termIterMax = 1;
00465     const std::string orthoType ("DGKS");
00466     const MagnitudeType orthoKappa = Teuchos::as<MagnitudeType> (-1.0);
00467     const int verbosity = Belos::Errors;
00468     const int outputStyle = Belos::General;
00469     const int outputFreq = -1;
00470     const std::string label ("Belos");
00471 
00472     RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00473     pl->set("Output Stream", outputStream,
00474       "is a reference-counted pointer to the output stream receiving\n"
00475       "all solver output.");
00476     pl->set("Lambda", lambda, "is the damping parameter.");
00477     pl->set("Rel RHS Err", relRhsErr,
00478       "estimates the error in the data defining the right-\n"
00479       "hand side.");
00480     pl->set("Rel Mat Err", relMatErr,
00481       "estimates the error in the data defining the matrix.");
00482     pl->set("Condition Limit", condMax,
00483       "bounds the estimated condition number of Abar.");
00484     pl->set("Maximum Iterations", maxIters,
00485       "allows at most the maximum number of iterations.");
00486     pl->set("Term Iter Max", termIterMax,
00487       "consecutive iterations meeting thresholds are necessary for\n"
00488        "for convergence.");
00489     pl->set("Orthogonalization", orthoType,
00490       "uses orthogonalization of either DGKS, ICGS, IMGS, or TSQR.");
00491     {
00492       OrthoManagerFactory<ScalarType, MV, OP> factory;
00493       pl->set("Orthogonalization", orthoType,
00494         "refers to the orthogonalization method to use.  Valid "
00495         "options: " + factory.validNamesString());
00496       RCP<const ParameterList> orthoParams = 
00497   factory.getDefaultParameters (orthoType);
00498       pl->set ("Orthogonalization Parameters", *orthoParams, 
00499          "Parameters specific to the type of orthogonalization used.");
00500     }
00501     pl->set("Orthogonalization Constant", orthoKappa,
00502       "is the threshold used by DGKS orthogonalization to determine\n"
00503       "whether or not to repeat classical Gram-Schmidt.  This parameter\n"
00504       "is ignored if \"Orthogonalization\" is not \"DGKS\".");
00505     pl->set("Verbosity", verbosity,
00506       "type(s) of solver information are outputted to the output\n"
00507       "stream.");
00508     pl->set("Output Style", outputStyle,
00509       "the style used for the solver information outputted to the\n"
00510       "output stream.");
00511     pl->set("Output Frequency", outputFreq,
00512       "is the frequency at which information is written to the\n"
00513       "output stream.");  
00514     pl->set("Timer Label", label,
00515       "is the string to use as a prefix for the timer labels.");
00516     //  pl->set("Restart Timers", restartTimers_);
00517     pl->set("Block Size", 1, "Block size parameter (currently, this must always be 1).");
00518 
00519     validParams_ = pl;
00520   }
00521   return validParams_;
00522 }
00523 
00524 
00525 template<class ScalarType, class MV, class OP>
00526 void 
00527 LSQRSolMgr<ScalarType,MV,OP>::
00528 setParameters (const Teuchos::RCP<Teuchos::ParameterList> &params)
00529 {
00530   using Teuchos::isParameterType;
00531   using Teuchos::getParameter;
00532   using Teuchos::null;
00533   using Teuchos::ParameterList;
00534   using Teuchos::parameterList;
00535   using Teuchos::RCP;
00536   using Teuchos::rcp;
00537   using Teuchos::rcp_dynamic_cast;
00538   using Teuchos::rcpFromRef;
00539   using Teuchos::Time;
00540   using Teuchos::TimeMonitor;
00541   using Teuchos::Exceptions::InvalidParameter;
00542   using Teuchos::Exceptions::InvalidParameterName;
00543   using Teuchos::Exceptions::InvalidParameterType;
00544 
00545   TEUCHOS_TEST_FOR_EXCEPTION(params.is_null(), std::invalid_argument,
00546          "Belos::LSQRSolMgr::setParameters: "
00547          "the input ParameterList is null.");
00548   RCP<const ParameterList> defaultParams = getValidParameters ();
00549   params->validateParametersAndSetDefaults (*defaultParams);
00550 
00551   // At this point, params is a valid parameter list with defaults
00552   // filled in.  Now we can "commit" it to our instance's parameter
00553   // list.
00554   params_ = params;
00555 
00556   // Get the damping (a.k.a. regularization) parameter lambda.
00557   lambda_ = params->get<MagnitudeType> ("Lambda");
00558 
00559   // Get the maximum number of iterations.
00560   maxIters_ = params->get<int>("Maximum Iterations");
00561 
00562   // (Re)set the timer label.
00563   {
00564     const std::string newLabel = params->get<std::string>("Timer Label");
00565 
00566     // Update parameter in our list and solver timer
00567     if (newLabel != label_) {
00568       label_ = newLabel;
00569     }
00570 
00571 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00572     std::string newSolveLabel = label_ + ": LSQRSolMgr total solve time";
00573     if (timerSolve_.is_null()) {
00574       // Ask TimeMonitor for a new timer.
00575       timerSolve_ = TimeMonitor::getNewCounter (newSolveLabel);
00576     } else {
00577       // We've already created a timer, but we may have changed its
00578       // label.  If we did change its name, then we have to forget
00579       // about the old timer and create a new one with a different
00580       // name.  This is because Teuchos::Time doesn't give you a way
00581       // to change a timer's name, once you've created it.  We assume
00582       // that if the user changed the timer's label, then the user
00583       // wants to reset the timer's results.
00584       const std::string oldSolveLabel = timerSolve_->name ();
00585 
00586       if (oldSolveLabel != newSolveLabel) {
00587   // Tell TimeMonitor to forget about the old timer.
00588   // TimeMonitor lets you clear timers by name.
00589   TimeMonitor::clearCounter (oldSolveLabel);
00590   timerSolve_ = TimeMonitor::getNewCounter (newSolveLabel);
00591       }
00592     }
00593 #endif // BELOS_TEUCHOS_TIME_MONITOR
00594   }
00595 
00596   // Check for a change in verbosity level
00597   {
00598     int newVerbosity = 0;
00599     // ParameterList gets confused sometimes about enums.  This
00600     // ensures that no matter how "Verbosity" was stored -- either an
00601     // an int, or as a Belos::MsgType enum, we will be able to extract
00602     // it.  If it was stored as some other type, we let the exception
00603     // through.
00604     try {
00605       newVerbosity = params->get<Belos::MsgType> ("Verbosity");
00606     } catch (Teuchos::Exceptions::InvalidParameterType&) {
00607       newVerbosity = params->get<int> ("Verbosity");
00608     }
00609     if (newVerbosity != verbosity_) {
00610       verbosity_ = newVerbosity;
00611     }
00612   }
00613 
00614   // (Re)set the output style.
00615   outputStyle_ = params->get<int> ("Output Style");
00616 
00617   // Get the output stream for the output manager.
00618   //
00619   // In case the output stream can't be read back in, we default to
00620   // stdout (std::cout), just to ensure reasonable behavior.
00621   {
00622     outputStream_ = params->get<RCP<std::ostream> > ("Output Stream");
00623 
00624     // We assume that a null output stream indicates that the user
00625     // doesn't want to print anything, so we replace it with a "black
00626     // hole" stream that prints nothing sent to it.  (We can't use a
00627     // null output stream, since the output manager always sends
00628     // things it wants to print to the output stream.)
00629     if (outputStream_.is_null())
00630       outputStream_ = rcp (new Teuchos::oblackholestream);
00631   }
00632 
00633   // Get the frequency of solver output.  (For example, -1 means
00634   // "never," 1 means "every iteration.")
00635   outputFreq_ = params->get<int> ("Output Frequency");
00636 
00637   // Create output manager if we need to, using the verbosity level
00638   // and output stream that we fetched above.  We do this here because
00639   // instantiating an OrthoManager using OrthoManagerFactory requires
00640   // a valid OutputManager.
00641   if (printer_.is_null()) {
00642     printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
00643   } else {
00644     printer_->setVerbosity (verbosity_);
00645     printer_->setOStream (outputStream_);
00646   }
00647 
00648   // Check if the orthogonalization changed, or if we need to
00649   // initialize it.
00650   typedef OrthoManagerFactory<ScalarType, MV, OP> factory_type;
00651   factory_type factory;
00652   bool mustMakeOrtho = false;
00653   {
00654     std::string tempOrthoType;
00655     try {
00656       tempOrthoType = params_->get<std::string> ("Orthogonalization");
00657     } catch (InvalidParameterName&) {
00658       tempOrthoType = orthoType_;
00659     }
00660     if (ortho_.is_null() || tempOrthoType != orthoType_) {
00661       mustMakeOrtho = true;
00662 
00663       // Ensure that the specified orthogonalization type is valid.
00664       if (! factory.isValidName (tempOrthoType)) {
00665   std::ostringstream os;
00666   os << "Belos::LSQRSolMgr: Invalid orthogonalization name \"" 
00667      << tempOrthoType << "\".  The following are valid options "
00668      << "for the \"Orthogonalization\" name parameter: ";
00669   factory.printValidNames (os);
00670   TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00671       }
00672       orthoType_ = tempOrthoType; // The name is valid, so accept it.
00673       params_->set ("Orthogonalization", orthoType_);
00674     }
00675   }
00676 
00677   // Get any parameters for the orthogonalization ("Orthogonalization
00678   // Parameters").  If not supplied, the orthogonalization manager
00679   // factory will supply default values.
00680   //
00681   // NOTE (mfh 21 Oct 2011) For the sake of backwards compatibility,
00682   // if params has an "Orthogonalization Constant" parameter and the
00683   // DGKS orthogonalization manager is to be used, the value of this
00684   // parameter will override DGKS's "depTol" parameter.
00685   //
00686   // Users must supply the orthogonalization manager parameters as a
00687   // sublist (supplying it as an RCP<ParameterList> would make the
00688   // resulting parameter list not serializable).
00689   RCP<ParameterList> orthoParams;
00690   { // The nonmember function returns an RCP<ParameterList>, 
00691     // which is what we want here.
00692     using Teuchos::sublist;
00693     // Abbreviation to avoid typos.
00694     const std::string paramName ("Orthogonalization Parameters");
00695 
00696     try {
00697       orthoParams = sublist (params_, paramName, true);
00698     } catch (InvalidParameter&) {
00699       // We didn't get the parameter list from params, so get a
00700       // default parameter list from the OrthoManagerFactory.
00701       // Modify params_ so that it has the default parameter list,
00702       // and set orthoParams to ensure it's a sublist of params_
00703       // (and not just a copy of one).
00704       params_->set (paramName, factory.getDefaultParameters (orthoType_));
00705       orthoParams = sublist (params_, paramName, true);
00706     }
00707   }
00708   TEUCHOS_TEST_FOR_EXCEPTION(orthoParams.is_null(), std::logic_error, 
00709            "Failed to get orthogonalization parameters.  "
00710            "Please report this bug to the Belos developers.");
00711 
00712   // If we need to, instantiate a new MatOrthoManager subclass
00713   // instance corresponding to the desired orthogonalization method.
00714   // We've already fetched the orthogonalization method name
00715   // (orthoType_) and its parameters (orthoParams) above.
00716   //
00717   // NOTE (mfh 21 Oct 2011) We only instantiate a new MatOrthoManager
00718   // subclass if the orthogonalization method name is different than
00719   // before.  Thus, for some orthogonalization managers, changes to
00720   // their parameters may not get propagated, if the manager type
00721   // itself didn't change.  The one exception is the "depTol"
00722   // (a.k.a. orthoKappa or "Orthogonalization Constant") parameter of
00723   // DGKS; changes to that _do_ get propagated down to the DGKS
00724   // instance.
00725   //
00726   // The most general way to fix this issue would be to supply each
00727   // orthogonalization manager class with a setParameterList() method
00728   // that takes a parameter list input, and changes the parameters as
00729   // appropriate.  A less efficient but correct way would be simply to
00730   // reinstantiate the OrthoManager every time, whether or not the
00731   // orthogonalization method name or parameters have changed.
00732   if (mustMakeOrtho) {
00733     // Create orthogonalization manager.  This requires that the
00734     // OutputManager (printer_) already be initialized.  LSQR
00735     // currently only orthogonalizes with respect to the Euclidean
00736     // inner product, so we set the inner product matrix M to null.
00737     RCP<const OP> M = null;
00738     ortho_ = factory.makeMatOrthoManager (orthoType_, M, printer_, 
00739             label_, orthoParams);
00740   }
00741   TEUCHOS_TEST_FOR_EXCEPTION(ortho_.is_null(), std::logic_error, 
00742            "The MatOrthoManager is not yet initialized, but "
00743            "should be by this point.  "
00744            "Please report this bug to the Belos developers.");
00745 
00746   // Check which orthogonalization constant to use.  We only need this
00747   // if orthoType_ == "DGKS" (and we already fetched the orthoType_
00748   // parameter above).
00749   if (orthoType_ == "DGKS") {
00750     if (params->isParameter ("Orthogonalization Constant")) {
00751       orthoKappa_ = params_->get<MagnitudeType> ("Orthogonalization Constant");
00752 
00753       if (orthoKappa_ > 0 && ! ortho_.is_null()) {
00754   typedef DGKSOrthoManager<ScalarType,MV,OP> ortho_impl_type;
00755   rcp_dynamic_cast<ortho_impl_type> (ortho_)->setDepTol (orthoKappa_);
00756       }
00757     }
00758   }
00759 
00760   // Check for condition number limit, number of consecutive passed
00761   // iterations, relative RHS error, and relative matrix error.
00762   // Create the LSQR convergence test if necessary.
00763   {
00764     condMax_ = params->get<MagnitudeType> ("Condition Limit");
00765     termIterMax_ = params->get<int>("Term Iter Max");
00766     relRhsErr_ = params->get<MagnitudeType> ("Rel RHS Err");
00767     relMatErr_ = params->get<MagnitudeType> ("Rel Mat Err");
00768 
00769     // Create the LSQR convergence test if it doesn't exist yet. 
00770     // Otherwise, update its parameters.
00771     if (convTest_.is_null()) {
00772       convTest_ = 
00773   rcp (new LSQRStatusTest<ScalarType,MV,OP> (condMax_, termIterMax_, 
00774                relRhsErr_, relMatErr_));
00775     } else {
00776       convTest_->setCondLim (condMax_);
00777       convTest_->setTermIterMax (termIterMax_);
00778       convTest_->setRelRhsErr (relRhsErr_);
00779       convTest_->setRelMatErr (relMatErr_);
00780     }
00781   }
00782 
00783   // Create the status test for maximum number of iterations if
00784   // necessary.  Otherwise, update it with the new maximum iteration
00785   // count.
00786   if (maxIterTest_.is_null()) {
00787     maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
00788   } else {
00789     maxIterTest_->setMaxIters (maxIters_);
00790   }
00791 
00792   // The stopping criterion is an OR combination of the test for
00793   // maximum number of iterations, and the LSQR convergence test.
00794   // ("OR combination" means that both tests will always be evaluated,
00795   // as opposed to a SEQ combination.)
00796   typedef StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00797   // If sTest_ is not null, then maxIterTest_ and convTest_ were
00798   // already constructed on entry to this routine, and sTest_ has
00799   // their pointers.  Thus, maxIterTest_ and convTest_ have gotten any
00800   // parameter changes, so we don't need to do anything to sTest_.
00801   if (sTest_.is_null()) {
00802     sTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::OR, 
00803            maxIterTest_, 
00804            convTest_));
00805   }
00806   
00807   if (outputTest_.is_null()) {
00808     // Create the status test output class.
00809     // This class manages and formats the output from the status test.
00810     StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
00811     outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_, 
00812              Passed + Failed + Undefined);
00813     // Set the solver string for the output test.
00814     const std::string solverDesc = " LSQR ";
00815     outputTest_->setSolverDesc (solverDesc);
00816   } else {
00817     // FIXME (mfh 18 Sep 2011) How do we change the output style of
00818     // outputTest_, without destroying and recreating it?
00819     outputTest_->setOutputManager (printer_);
00820     outputTest_->setChild (sTest_);
00821     outputTest_->setOutputFrequency (outputFreq_);
00822     // Since outputTest_ can only be created here, I'm assuming that
00823     // the fourth constructor argument ("printStates") was set
00824     // correctly on constrution; I don't need to reset it (and I can't
00825     // set it anyway, given StatusTestOutput's interface).
00826   }
00827 
00828   // Inform the solver manager that the current parameters were set.
00829   isSet_ = true;
00830 }
00831 
00832 
00833 template<class ScalarType, class MV, class OP>
00834 Belos::ReturnType LSQRSolMgr<ScalarType,MV,OP>::solve() {
00835   using Teuchos::RCP;
00836   using Teuchos::rcp;
00837 
00838   // Set the current parameters if they were not set before.  NOTE:
00839   // This may occur if the user generated the solver manager with the
00840   // default constructor, but did not set any parameters using
00841   // setParameters().
00842   if (!isSet_) {
00843     setParameters (Teuchos::parameterList (*getValidParameters()));
00844   }
00845 
00846   TEUCHOS_TEST_FOR_EXCEPTION(problem_.is_null(), LSQRSolMgrLinearProblemFailure,
00847          "The linear problem to solve is null.");
00848   TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), LSQRSolMgrLinearProblemFailure,
00849                      "LSQRSolMgr::solve(): The linear problem is not ready, "
00850          "as its setProblem() method has not been called.");
00851   TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs (*(problem_->getRHS ())) != 1, 
00852          LSQRSolMgrBlockSizeFailure, 
00853          "LSQRSolMgr::solve(): The current implementation of LSQR "
00854          "only knows how to solve problems with one right-hand "
00855          "side, but the linear problem to solve has " 
00856          << MVT::GetNumberVecs (*(problem_->getRHS ())) 
00857          << " right-hand sides.");
00858 
00859   // We've validated the LinearProblem instance above.  If any of the
00860   // StatusTests needed to be initialized using information from the
00861   // LinearProblem, now would be the time to do so.  (This is true of
00862   // GMRES, where the residual convergence test(s) to instantiate
00863   // depend on knowing whether there is a left preconditioner.  This
00864   // is why GMRES has an "isSTSet_" Boolean member datum, which tells
00865   // you whether the status tests have been instantiated and are ready
00866   // for use.
00867 
00868   // test isFlexible might go here.
00869 
00870   // Next the right-hand sides to solve are identified.  Among other things,
00871   // this enables getCurrLHSVec() to get the current initial guess vector,
00872   // and getCurrRHSVec() to get the current right-hand side (in Iter).
00873   std::vector<int> currRHSIdx(1, 0); 
00874   problem_->setLSIndex(currRHSIdx);
00875 
00876   // Reset the status test.  
00877   outputTest_->reset();
00878 
00879   // Don't assume convergence unless we've verified that the
00880   // convergence test passed.
00881   bool isConverged = false;
00882 
00883   // FIXME: Currently we are setting the initial guess to zero, since
00884   // the solver doesn't yet know how to handle a nonzero initial
00885   // guess.  This could be fixed by rewriting the solver to work with
00886   // the residual and a delta.
00887   //
00888   // In a least squares problem with a nonzero initial guess, the
00889   // minimzation problem involves the distance (in a norm depending on
00890   // the preconditioner) between the solution and the the initial
00891   // guess.
00892 
00894   // Solve the linear problem using LSQR
00896 
00897   // Parameter list for the LSQR iteration.
00898   Teuchos::ParameterList plist;
00899 
00900   // Use the solver manager's "Lambda" parameter to set the
00901   // iteration's "Lambda" parameter.  We know that the solver
00902   // manager's parameter list (params_) does indeed contain the
00903   // "Lambda" parameter, because solve() always ensures that
00904   // setParameters() has been called.
00905   plist.set ("Lambda", lambda_);
00906 
00907   typedef LSQRIter<ScalarType,MV,OP> iter_type;
00908   RCP<iter_type> lsqr_iter = 
00909     rcp (new iter_type (problem_, printer_, outputTest_, plist));
00910 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00911   Teuchos::TimeMonitor slvtimer(*timerSolve_);
00912 #endif
00913 
00914   // Reset the number of iterations.
00915   lsqr_iter->resetNumIters();
00916   // Reset the number of calls that the status test output knows about.
00917   outputTest_->resetNumCalls();
00918   // Set the new state and initialize the solver.
00919   LSQRIterationState<ScalarType,MV> newstate;
00920   lsqr_iter->initializeLSQR(newstate);
00921   // tell lsqr_iter to iterate
00922   try {
00923     lsqr_iter->iterate();    
00924 
00925     // First check for convergence.  If we didn't converge, then check
00926     // whether we reached the maximum number of iterations.  If
00927     // neither of those happened, there must have been a bug.
00928     if (convTest_->getStatus() == Belos::Passed) {
00929       isConverged = true;
00930     } else if (maxIterTest_->getStatus() == Belos::Passed) {
00931       isConverged = false;
00932     } else {
00933       TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
00934        "LSQRSolMgr::solve(): LSQRIteration::iterate() "
00935        "returned without either the convergence test or "
00936        "the maximum iteration count test passing."
00937        "Please report this bug to the Belos developers.");
00938     }
00939   } catch (const std::exception &e) {
00940     printer_->stream(Belos::Errors) << "Error! Caught std::exception in LSQRIter::iterate() at iteration " 
00941             << lsqr_iter->getNumIters() << std::endl 
00942             << e.what() << std::endl;
00943     throw;
00944   }
00945       
00946   // identify current linear system as solved LinearProblem 
00947   problem_->setCurrLS();
00948   // print final summary
00949   sTest_->print (printer_->stream (Belos::FinalSummary));
00950 
00951   // Print timing information, if the corresponding compile-time and
00952   // run-time options are enabled.
00953 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00954   // Calling summarize() can be expensive, so don't call unless the
00955   // user wants to print out timing details.  summarize() will do all
00956   // the work even if it's passed a "black hole" output stream.
00957   if (verbosity_ & TimingDetails)
00958     Teuchos::TimeMonitor::summarize (printer_->stream (Belos::TimingDetails));
00959 #endif // BELOS_TEUCHOS_TIME_MONITOR
00960 
00961   // A posteriori solve information
00962   numIters_ = maxIterTest_->getNumIters();
00963   matCondNum_ = convTest_->getMatCondNum();
00964   matNorm_ = convTest_->getMatNorm();
00965   resNorm_ = convTest_->getResidNorm();
00966   matResNorm_ = convTest_->getLSResidNorm();
00967  
00968   if (! isConverged) {
00969     return Belos::Unconverged;
00970   } else {
00971     return Belos::Converged;
00972   }
00973 }
00974 
00975 // LSQRSolMgr requires the solver manager to return an eponymous std::string.
00976 template<class ScalarType, class MV, class OP>
00977 std::string LSQRSolMgr<ScalarType,MV,OP>::description() const
00978 {
00979   std::ostringstream oss;
00980   oss << "LSQRSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
00981   oss << "{";
00982   oss << "Ortho Type='"<<orthoType_<<"'";
00983   oss << ", Lambda="<< lambda_;
00984   oss << ", condition number limit="<< condMax_;
00985   oss << ", relative RHS Error="<< relRhsErr_;
00986   oss << ", relative Matrix Error="<< relMatErr_;
00987   oss << ", maximum number of iterations="<< maxIters_;
00988   oss << ", termIterMax="<<termIterMax_;
00989   oss << "}";
00990   return oss.str();
00991 }
00992   
00993 } // end Belos namespace
00994 
00995 #endif /* BELOS_LSQR_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines