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 
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 "BelosDGKSOrthoManager.hpp"
00058 #include "BelosICGSOrthoManager.hpp"
00059 #include "BelosIMGSOrthoManager.hpp"
00060 #include "BelosStatusTestMaxIters.hpp"
00061 #include "BelosLSQRStatusTest.hpp"
00062 #include "BelosStatusTestCombo.hpp"
00063 #include "BelosStatusTestOutputFactory.hpp"
00064 #include "BelosOutputManager.hpp"
00065 #include "Teuchos_BLAS.hpp"
00066 #include "Teuchos_LAPACK.hpp"
00067 
00068 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00069 #include "Teuchos_TimeMonitor.hpp"
00070 #endif
00071 
00082 namespace Belos {
00083 
00084   
00086 
00087 
00094 class LSQRSolMgrLinearProblemFailure : public BelosError {public:
00095     LSQRSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
00096   {}};
00097 
00104 class LSQRSolMgrOrthoFailure : public BelosError {public:
00105     LSQRSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00106   {}};
00107 
00114 class LSQRSolMgrBlockSizeFailure : public BelosError {public:
00115     LSQRSolMgrBlockSizeFailure(const std::string& what_arg) : BelosError(what_arg)
00116   {}};
00117 
00118 template<class ScalarType, class MV, class OP>
00119 class LSQRSolMgr : public SolverManager<ScalarType,MV,OP> {
00120   
00121 private:
00122   typedef MultiVecTraits<ScalarType,MV> MVT;
00123   typedef OperatorTraits<ScalarType,MV,OP> OPT;
00124   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00125   typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00126   typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00127   
00128 public:
00129   
00131 
00132   
00139   LSQRSolMgr();
00140   
00161   LSQRSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00162      const Teuchos::RCP<Teuchos::ParameterList> &pl );
00163   
00165   virtual ~LSQRSolMgr() {};
00167   
00169 
00170 
00173   const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00174     return *problem_;
00175   }
00176   
00179   Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00180   
00183   Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
00184   
00190   Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00191     return Teuchos::tuple(timerSolve_);
00192   }
00193   
00195   int getNumIters() const {
00196     return numIters_;
00197   }
00198   
00201   bool isLOADetected() const { return false; }
00202   
00204   
00206 
00207    
00209   void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; }
00210   
00212   void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
00213     
00215    
00217 
00218 
00222   void reset( const ResetType type ) { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
00224  
00226 
00227     
00246   ReturnType solve();
00247     
00249     
00252     
00254   std::string description() const;
00255     
00257     
00258 private:
00259 
00260 
00261 
00262   // Method to convert std::string to enumerated type for residual.
00263   Belos::ScaleType convertStringToScaleType( std::string& scaleType ) {
00264     if (scaleType == "Norm of Initial Residual") {
00265       return Belos::NormOfInitRes;
00266     } else if (scaleType == "Norm of Preconditioned Initial Residual") {
00267       return Belos::NormOfPrecInitRes;
00268     } else if (scaleType == "Norm of RHS") {
00269       return Belos::NormOfRHS;
00270     } else if (scaleType == "None") {
00271       return Belos::None;
00272     } else
00273       TEST_FOR_EXCEPTION( true ,std::logic_error,
00274         "Belos::BlockGmresSolMgr(): Invalid residual scaling type.");
00275   }
00276 
00277   // Method for checking current status test against defined linear problem.
00278   //bool checkStatusTest(); ???only needed for gmres ???
00279 
00280   // Linear problem.
00281   Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00282   
00283   // Output manager.
00284   Teuchos::RCP<OutputManager<ScalarType> > printer_;
00285   Teuchos::RCP<std::ostream> outputStream_;
00286   
00287   // Status test.
00288   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00289   Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00290   Teuchos::RCP<LSQRStatusTest<ScalarType,MV,OP> > convTest_;
00291   Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00292   
00293   // Orthogonalization manager.
00294   Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_; 
00295     
00296   // Current parameter list.
00297   Teuchos::RCP<Teuchos::ParameterList> params_;
00298 
00299     
00300   // Default solver values.
00301   static const ScalarType lambda_default_;  // MagnitueType, lsqr only
00302   static const MagnitudeType relRhsErr_default_; // lsqr only
00303   static const MagnitudeType relMatErr_default_; // lsqr only 
00304   static const MagnitudeType condlim_default_;   // lsqr only
00305   static const int maxIters_default_; 
00306   static const int termIterMax_default_;         // lsqr only
00307   static const std::string orthoType_default_;
00308   static const MagnitudeType orthoKappa_default_;
00309   static const int verbosity_default_;
00310   static const int outputStyle_default_;
00311   static const int outputFreq_default_;
00312   static const Teuchos::RCP<std::ostream> outputStream_default_;
00313   static const std::string label_default_;
00314   
00315   // Current solver values.
00316   ScalarType lambda_;
00317   MagnitudeType relRhsErr_, relMatErr_, condlim_;
00318   int maxIters_, termIterMax_;
00319   int numIters_;
00320   std::string orthoType_; 
00321   MagnitudeType orthoKappa_;
00322   int verbosity_, outputStyle_, outputFreq_;
00323     
00324   // Timers.
00325   std::string label_;
00326   Teuchos::RCP<Teuchos::Time> timerSolve_;
00327 
00328   // Internal state variables.
00329   bool isSet_;
00330   bool isSTSet_;
00331   bool loaDetected_;
00332 
00333 };
00334 
00335 
00336 // Default solver values.
00337 template<class ScalarType, class MV, class OP>
00338 const ScalarType LSQRSolMgr<ScalarType,MV,OP>::lambda_default_ = 0.0;
00339 
00340 template<class ScalarType, class MV, class OP>
00341 const typename LSQRSolMgr<ScalarType,MV,OP>::MagnitudeType LSQRSolMgr<ScalarType,MV,OP>::relRhsErr_default_ = 0.0;
00342 
00343 template<class ScalarType, class MV, class OP>
00344 const typename LSQRSolMgr<ScalarType,MV,OP>::MagnitudeType LSQRSolMgr<ScalarType,MV,OP>::relMatErr_default_ = 0.0;
00345 
00346 template<class ScalarType, class MV, class OP>
00347 const typename LSQRSolMgr<ScalarType,MV,OP>::MagnitudeType LSQRSolMgr<ScalarType,MV,OP>::condlim_default_ = 0.0;
00348 
00349 template<class ScalarType, class MV, class OP>
00350 const int LSQRSolMgr<ScalarType,MV,OP>::maxIters_default_ = 1000;
00351 
00352 template<class ScalarType, class MV, class OP>
00353 const int LSQRSolMgr<ScalarType,MV,OP>::termIterMax_default_ = 1;
00354 
00355 template<class ScalarType, class MV, class OP>
00356 const std::string LSQRSolMgr<ScalarType,MV,OP>::orthoType_default_ = "DGKS";
00357 
00358 template<class ScalarType, class MV, class OP>
00359 const typename LSQRSolMgr<ScalarType,MV,OP>::MagnitudeType LSQRSolMgr<ScalarType,MV,OP>::orthoKappa_default_ = -1.0;
00360 
00361 template<class ScalarType, class MV, class OP>
00362 const int LSQRSolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors;
00363 
00364 template<class ScalarType, class MV, class OP>
00365 const int LSQRSolMgr<ScalarType,MV,OP>::outputStyle_default_ = Belos::General;
00366 
00367 template<class ScalarType, class MV, class OP>
00368 const int LSQRSolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1;
00369 
00370 template<class ScalarType, class MV, class OP>
00371 const Teuchos::RCP<std::ostream> LSQRSolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);
00372 
00373 template<class ScalarType, class MV, class OP>
00374 const std::string LSQRSolMgr<ScalarType,MV,OP>::label_default_ = "Belos";
00375 
00376 // Empty Constructor
00377 template<class ScalarType, class MV, class OP>
00378 LSQRSolMgr<ScalarType,MV,OP>::LSQRSolMgr() :
00379   outputStream_(outputStream_default_),
00380   lambda_(lambda_default_),
00381   relRhsErr_(relRhsErr_default_),
00382   relMatErr_(relMatErr_default_),
00383   condlim_(condlim_default_),
00384   maxIters_(maxIters_default_),
00385   termIterMax_(termIterMax_default_),
00386   orthoType_(orthoType_default_),
00387   orthoKappa_(orthoKappa_default_),
00388   verbosity_(verbosity_default_),
00389   outputStyle_(outputStyle_default_),
00390   outputFreq_(outputFreq_default_),
00391   label_(label_default_),
00392   isSet_(false),
00393   isSTSet_(false),
00394   loaDetected_(false)
00395 {}
00396 
00397 
00398 // Basic Constructor
00399 template<class ScalarType, class MV, class OP>
00400 LSQRSolMgr<ScalarType,MV,OP>::LSQRSolMgr( 
00401   const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00402   const Teuchos::RCP<Teuchos::ParameterList> &pl ) : 
00403   problem_(problem),
00404   outputStream_(outputStream_default_),
00405   lambda_(lambda_default_),
00406   relRhsErr_(relRhsErr_default_),
00407   relMatErr_(relMatErr_default_),
00408   condlim_(condlim_default_),
00409   maxIters_(maxIters_default_),
00410   termIterMax_(termIterMax_default_),
00411   orthoType_(orthoType_default_),
00412   orthoKappa_(orthoKappa_default_),
00413   verbosity_(verbosity_default_),
00414   outputStyle_(outputStyle_default_),
00415   outputFreq_(outputFreq_default_),
00416   label_(label_default_),
00417   isSet_(false),
00418   isSTSet_(false),
00419   loaDetected_(false)
00420 {
00421   TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
00422 
00423   // If the parameter list pointer is null, 
00424   // then set the current parameters to the default parameter list.
00425   if ( !is_null(pl) ) {
00426     setParameters( pl );  
00427   }
00428 }
00429 
00430 
00431 template<class ScalarType, class MV, class OP>
00432 Teuchos::RCP<const Teuchos::ParameterList> 
00433 LSQRSolMgr<ScalarType,MV,OP>::getValidParameters() const
00434 {
00435   static Teuchos::RCP<const Teuchos::ParameterList> validPL;
00436   
00437   // Set all the valid parameters and their default values.
00438   if(is_null(validPL)) {
00439     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00440     pl->set("Output Stream", outputStream_default_,
00441       "is a reference-counted pointer to the output stream receiving\n"
00442       "all solver output.");
00443     pl->set("Lambda", lambda_default_, "is the damping parameter.");
00444     pl->set("Rel RHS Err", relRhsErr_default_,
00445       "estimates the error in the data defining the rignt-\n"
00446       "hand side.");
00447     pl->set("Rel Mat Err", relMatErr_default_,
00448       "estimates the error in the data defining the matrix.");
00449     pl->set("Condition limit", condlim_default_,
00450       "bounds the estimated condition number of Abar.");
00451     pl->set("Maximum Iterations", maxIters_default_,
00452       "allows at most the maximum number of iterations.");
00453     pl->set("Term Iter Max", termIterMax_default_,
00454       "consecutive iterations meeting thresholds are necessary for\n"
00455        "for convergence.");
00456     pl->set("Orthogonalization", orthoType_default_,
00457       "uses orthogonalization of either DGKS, ICGS or IMGS.");
00458     pl->set("Orthogonalization Constant",orthoKappa_default_,
00459       "is the threshold used by DGKS orthogonalization to determine\n"
00460       "whether or not to repeat classical Gram-Schmidt.");
00461     pl->set("Verbosity", verbosity_default_,
00462       "type(s) of solver information are outputted to the output\n"
00463       "stream.");
00464     pl->set("Output Style", outputStyle_default_,
00465       "the style used for the solver information outputted to the\n"
00466       "output stream.");
00467     pl->set("Output Frequency", outputFreq_default_,
00468       "is the frequency at which information is written to the\n"
00469       "output stream.");  
00470     pl->set("Timer Label", label_default_,
00471       "is the string to use as a prefix for the timer labels.");
00472     //  pl->set("Restart Timers", restartTimers_);
00473     validPL = pl;
00474   }
00475   return validPL;
00476 }
00477 
00478 
00479 template<class ScalarType, class MV, class OP>
00480 void LSQRSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
00481 {
00482   // Create the internal parameter list if ones doesn't already exist.
00483   if (params_ == Teuchos::null) {
00484     params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
00485   }
00486   else {
00487     params->validateParameters(*getValidParameters());
00488   }
00489 
00490   // Check for damping value lambda
00491   if (params->isParameter("Lambda")) {
00492     lambda_ = params->get("Lambda",lambda_default_);
00493 
00494     // Update parameter in our list and in status test.
00495     params_->set("Lambda",lambda_);
00496   }
00497 
00498   // Check for maximum number of iterations
00499   if (params->isParameter("Maximum Iterations")) {
00500     maxIters_ = params->get("Maximum Iterations",maxIters_default_);
00501 
00502     // Update parameter in our list and in status test.
00503     params_->set("Maximum Iterations", maxIters_);
00504     if (maxIterTest_!=Teuchos::null)
00505       maxIterTest_->setMaxIters( maxIters_ );
00506   }
00507 
00508   // Check to see if the timer label changed.
00509   if (params->isParameter("Timer Label")) {
00510     std::string tempLabel = params->get("Timer Label", label_default_);
00511 
00512     // Update parameter in our list and solver timer
00513     if (tempLabel != label_) {
00514       label_ = tempLabel;
00515       params_->set("Timer Label", label_);
00516       std::string solveLabel = label_ + ": LSQRSolMgr total solve time";
00517 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00518       timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00519 #endif
00520     }
00521   }
00522   // Flexible LSQR undiscovered
00523 
00524   // Check if the orthogonalization changed.
00525   if (params->isParameter("Orthogonalization")) {
00526     std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
00527     TEST_FOR_EXCEPTION( tempOrthoType != "DGKS" && tempOrthoType != "ICGS" && tempOrthoType != "IMGS", 
00528                         std::invalid_argument,
00529       "LSQRSolMgr: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
00530     if (tempOrthoType != orthoType_) {
00531       orthoType_ = tempOrthoType;
00532       // Create orthogonalization manager
00533       if (orthoType_=="DGKS") {
00534   if (orthoKappa_ <= 0) {
00535     ortho_ = Teuchos::rcp( new Belos::DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00536   }
00537   else {
00538     ortho_ = Teuchos::rcp( new Belos::DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00539     Teuchos::rcp_dynamic_cast<Belos::DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00540   }
00541       }
00542       else if (orthoType_=="ICGS") {
00543   ortho_ = Teuchos::rcp( new Belos::ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00544       } 
00545       else if (orthoType_=="IMGS") {
00546   ortho_ = Teuchos::rcp( new Belos::IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00547       } 
00548     }  
00549   }
00550 
00551   // Check which orthogonalization constant to use.
00552   if (params->isParameter("Orthogonalization Constant")) {
00553     orthoKappa_ = params->get("Orthogonalization Constant",orthoKappa_default_);
00554 
00555     // Update parameter in our list.
00556     params_->set("Orthogonalization Constant",orthoKappa_);
00557     if (orthoType_=="DGKS") {
00558       if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
00559   Teuchos::rcp_dynamic_cast<Belos::DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00560       }
00561     } 
00562   }
00563 
00564   // Check for a change in verbosity level
00565   if (params->isParameter("Verbosity")) {
00566     if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
00567       verbosity_ = params->get("Verbosity", verbosity_default_);
00568     } else {
00569       verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
00570     }
00571 
00572     // Update parameter in our list.
00573     params_->set("Verbosity", verbosity_);
00574     if (printer_ != Teuchos::null)
00575       printer_->setVerbosity(verbosity_);
00576   }
00577 
00578   // Check for a change in output style
00579   if (params->isParameter("Output Style")) {
00580     if (Teuchos::isParameterType<int>(*params,"Output Style")) {
00581       outputStyle_ = params->get("Output Style", outputStyle_default_);
00582     } else {
00583       outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
00584     }
00585 
00586     // Update parameter in our list.
00587     params_->set("Output Style", outputStyle_);
00588     outputTest_ = Teuchos::null;
00589   }
00590 
00591   // output stream
00592   if (params->isParameter("Output Stream")) {
00593     outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
00594 
00595     // Update parameter in our list.
00596     params_->set("Output Stream", outputStream_);
00597     if (printer_ != Teuchos::null)
00598       printer_->setOStream( outputStream_ );
00599   }
00600 
00601   // frequency level
00602   if (verbosity_ & Belos::StatusTestDetails) {
00603     if (params->isParameter("Output Frequency")) {
00604       outputFreq_ = params->get("Output Frequency", outputFreq_default_);
00605     }
00606 
00607     // Update parameter in out list and output status test.
00608     params_->set("Output Frequency", outputFreq_);
00609     if (outputTest_ != Teuchos::null)
00610       outputTest_->setOutputFrequency( outputFreq_ );
00611   }
00612 
00613   // Create output manager if we need to.
00614   if (printer_ == Teuchos::null) {
00615     printer_ = Teuchos::rcp( new Belos::OutputManager<ScalarType>(verbosity_, outputStream_) );
00616   }  
00617   
00618   // Check for convergence tolerance
00619   if (params->isParameter("Condition Limit")) {
00620     condlim_ = params->get("Condition Limit",condlim_default_);
00621 
00622     // Update parameter in our list and residual tests.
00623     params_->set("Condition Limit", condlim_);
00624     if (convTest_ != Teuchos::null)
00625       convTest_->setCondLim( condlim_ );
00626   }
00627 
00628   // dmd: Residual scaling would go here.  What are explicit and implicit residual scaling? 
00629   // dmd: No Flexible LSQR -> No Implicit residual -> No residual scaling?
00630 
00631   // Check for number of consecutive passed iterations
00632   if (params->isParameter("Term Iter Max")) {
00633     termIterMax_ = params->get("Term Iter Max", termIterMax_default_);
00634 
00635     // Update parameter in our list and residual tests
00636     params_->set("Term Iter Max", termIterMax_);
00637     if (convTest_ != Teuchos::null)
00638       convTest_->setTermIterMax( termIterMax_ );
00639   }
00640 
00641   // Check for relative RHS error
00642   if (params->isParameter("Rel RHS Err")) {
00643     relRhsErr_ = params->get("Rel RHS Err",relRhsErr_default_);
00644 
00645     // Update parameter in our list and residual tests
00646     params_->set("Rel RHS Err", relRhsErr_);
00647     if (convTest_ != Teuchos::null)
00648       convTest_->setRelRhsErr( relRhsErr_ );
00649   }
00650 
00651   // Check for relative matrix error
00652   if (params->isParameter("Rel Mat Err")) {
00653     relMatErr_ = params->get("Rel Mat Err",relMatErr_default_);
00654 
00655     // Update parameter in our list and residual tests
00656     params_->set("Rel Mat Err", relMatErr_);
00657     if (convTest_ != Teuchos::null)
00658       convTest_->setRelMatErr( relMatErr_ );
00659   }
00660 
00661   // Create status tests if we need to.
00662 
00663   // Basic test checks maximum iterations and native residual.
00664   if (maxIterTest_ == Teuchos::null)
00665     maxIterTest_ = Teuchos::rcp( new Belos::StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
00666 
00667   // Implicit residual test on the native residual to check for convergence.
00668   if (convTest_ == Teuchos::null)
00669     convTest_ = Teuchos::rcp( new LSQRStatusTest<ScalarType,MV,OP>(condlim_, termIterMax_, relRhsErr_, relMatErr_) );
00670   
00671   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00672 
00673   if (sTest_ == Teuchos::null)
00674     sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
00675   
00676   if (outputTest_ == Teuchos::null) {
00677 
00678     // Create the status test output class.
00679     // This class manages and formats the output from the status test.
00680     Belos::StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
00681     outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Belos::Passed+Belos::Failed+Belos::Undefined );
00682 
00683     // Set the solver string for the output test
00684     std::string solverDesc = " LSQR ";
00685     outputTest_->setSolverDesc( solverDesc );
00686 
00687   }
00688   // dmd: divergence from BlockGmresSolMgr end here 
00689 
00690   // Create orthogonalization manager if we need to.
00691   if (ortho_ == Teuchos::null) {
00692     if (orthoType_=="DGKS") {
00693       if (orthoKappa_ <= 0) {
00694   ortho_ = Teuchos::rcp( new Belos::DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00695       }
00696       else {
00697   ortho_ = Teuchos::rcp( new Belos::DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00698   Teuchos::rcp_dynamic_cast<Belos::DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00699       }
00700     }
00701     else if (orthoType_=="ICGS") {
00702       ortho_ = Teuchos::rcp( new Belos::ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00703     } 
00704     else if (orthoType_=="IMGS") {
00705       ortho_ = Teuchos::rcp( new Belos::IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00706     } 
00707     else {
00708       TEST_FOR_EXCEPTION(orthoType_!="ICGS"&&orthoType_!="DGKS"&&orthoType_!="IMGS",std::logic_error,
00709        "LSQRSolMgr(): Invalid orthogonalization type.");
00710     }  
00711   }
00712 
00713   // Create the timer if we need to.
00714   if (timerSolve_ == Teuchos::null) {
00715     std::string solveLabel = label_ + ": LSQRSolMgr total solve time";
00716 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00717     timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00718 #endif
00719   }
00720 
00721   // Inform the solver manager that the current parameters were set.
00722   isSet_ = true;
00723 }
00724 
00725 
00726 // Check the status test versus the defined linear problem
00727 //bool LSQRSolMgr<ScalarType,MV,OP>::checkStatusTest()
00728 
00729     
00730 // solve()
00731 template<class ScalarType, class MV, class OP>
00732 Belos::ReturnType LSQRSolMgr<ScalarType,MV,OP>::solve() {
00733 
00734   // Set the current parameters if they were not set before.
00735   // NOTE:  This may occur if the user generated the solver manager with the default constructor
00736   // but did not set any parameters using setParameters().
00737   if (!isSet_) {
00738     setParameters(Teuchos::parameterList(*getValidParameters()));
00739   }
00740 
00741   //Teuchos::BLAS<int,ScalarType> blas;
00742   //Teuchos::LAPACK<int,ScalarType> lapack;
00743   
00744   TEST_FOR_EXCEPTION(!problem_->isProblemSet(),LSQRSolMgrLinearProblemFailure,
00745                      "LSQRSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
00746   TEST_FOR_EXCEPTION(MVT::GetNumberVecs( *(problem_->getRHS()) ) != 1, LSQRSolMgrBlockSizeFailure, "LSQRSolMgr::solve(): Incorrect number of RHS vectors, should be exactly 1.");
00747 
00748   // test isFlexible might go here.
00749   // isSTSet abbreviates is Status Test set.
00750   // Inform the linear problem of the current linear system to solve.
00751   std::vector<int> currRHSIdx(1, 0);
00752   problem_->setLSIndex(currRHSIdx);
00753   // Parameter list
00754   Teuchos::ParameterList plist;
00755   plist.set("Lambda", lambda_);
00756   // Reset the status test.  
00757   outputTest_->reset();
00758   // Assume convergence is achieved, then let any failed convergence set this to false.
00759   bool isConverged = true;  
00760 
00761   // LSQR requires a zero initial guess,  ensure vanishing lhs
00762   Teuchos::RCP<MV> x = problem_->getLHS();
00763   // Set x to the zero vector
00764   MVT::MvInit(*x, SCT::zero());
00765   problem_->updateSolution(x);
00766 
00768   // LSQR solver
00769   Teuchos::RCP<LSQRIter<ScalarType,MV,OP> > lsqr_iter = Teuchos::rcp( new LSQRIter<ScalarType,MV,OP>(problem_, printer_, outputTest_, plist) );
00770 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00771   Teuchos::TimeMonitor slvtimer(*timerSolve_);
00772 #endif
00773 
00774 
00775   // Reset the number of iterations.
00776   lsqr_iter->resetNumIters();
00777   // Reset the number of calls that the status test output knows about.
00778   outputTest_->resetNumCalls();
00779   // Set the new state and initialize the solver.
00780   LSQRIterationState<ScalarType,MV> newstate;
00781   lsqr_iter->initializeLSQR(newstate);
00782   // tell lsqr_iter to iterate
00783   try {
00784     lsqr_iter->iterate();    
00785 
00787     //
00788     // check convergence first
00789     //
00791     if ( convTest_->getStatus() == Belos::Passed ) {
00792     }
00793     else if ( maxIterTest_->getStatus() == Belos::Passed ) {
00794       // we don't have convergence
00795       isConverged = false;
00796     }
00798     //
00799     // we returned from iterate(), but none of our status tests Passed.
00800     // something is wrong, and it is probably our fault.
00801     //
00803     else {
00804       TEST_FOR_EXCEPTION(true,std::logic_error,
00805        "LSQRSolMgr::solve(): Invalid return from LSQRIteration::iterate().");
00806     }
00807   }
00808   catch (const std::exception &e) {
00809     printer_->stream(Belos::Errors) << "Error! Caught std::exception in LSQRIter::iterate() at iteration " 
00810             << lsqr_iter->getNumIters() << std::endl 
00811             << e.what() << std::endl;
00812     throw;
00813   }
00814       
00815   // Inform the linear problem that we are finished with this linear system.
00816   problem_->setCurrLS();
00817   // print final summary
00818   sTest_->print( printer_->stream(Belos::FinalSummary) );
00819 
00820   // Print timing information, if the corresponding compile-time and
00821   // run-time options are enabled.
00822 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00823   // Calling summarize() can be expensive, so don't call unless the
00824   // user wants to print out timing details.  summarize() will do all
00825   // the work even if it's passed a "black hole" output stream.
00826   if (verbosity_ & TimingDetails)
00827     Teuchos::TimeMonitor::summarize( printer_->stream(Belos::TimingDetails) );
00828 #endif // BELOS_TEUCHOS_TIME_MONITOR
00829 
00830   // get iteration information for this solve
00831   numIters_ = maxIterTest_->getNumIters();
00832  
00833   if (!isConverged) {
00834     return Belos::Unconverged; // return from LSQRSolMgr::solve() 
00835   }
00836   return Belos::Converged; // return from LSQRSolMgr::solve() 
00837 }
00838 
00839 //  This method requires the solver manager to return a std::string that describes itself.
00840 template<class ScalarType, class MV, class OP>
00841 std::string LSQRSolMgr<ScalarType,MV,OP>::description() const
00842 {
00843   std::ostringstream oss;
00844   oss << "LSQRSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
00845   oss << "{";
00846   oss << "Ortho Type='"<<orthoType_<<"'";
00847   oss << ", Lambda="<< lambda_;
00848   oss << ", condition number limit="<< condlim_;
00849   oss << ", relative RHS Error="<< relRhsErr_;
00850   oss << ", relative Matrix Error="<< relMatErr_;
00851   oss << ", maximum number of iterations="<< maxIters_;
00852   oss << ", termIterMax?="<<termIterMax_;
00853   oss << "}";
00854   return oss.str();
00855 }
00856   
00857 } // end Belos namespace
00858 
00859 #endif /* BELOS_LSQR_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines