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