Belos Version of the Day
BelosMinresSolMgr.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_MINRES_SOLMGR_HPP
00043 #define BELOS_MINRES_SOLMGR_HPP
00044 
00047 
00048 #include "BelosConfigDefs.hpp"
00049 #include "BelosTypes.hpp"
00050 
00051 #include "BelosLinearProblem.hpp"
00052 #include "BelosSolverManager.hpp"
00053 
00054 #include "BelosMinresIter.hpp"
00055 #include "BelosStatusTestMaxIters.hpp"
00056 #include "BelosStatusTestGenResNorm.hpp"
00057 #include "BelosStatusTestCombo.hpp"
00058 #include "BelosStatusTestOutputFactory.hpp"
00059 #include "BelosOutputManager.hpp"
00060 #include "Teuchos_BLAS.hpp"
00061 #include "Teuchos_LAPACK.hpp"
00062 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00063 #include "Teuchos_TimeMonitor.hpp"
00064 #endif
00065 
00066 #include "Teuchos_StandardParameterEntryValidators.hpp"
00067 // Teuchos::ScalarTraits<int> doesn't define rmax(), alas, so we get
00068 // INT_MAX from here.
00069 #include <climits> 
00070 
00071 namespace Belos {
00072   
00074 
00075 
00085   //
00089   class MinresSolMgrLinearProblemFailure : public BelosError {
00090   public:
00091     MinresSolMgrLinearProblemFailure (const std::string& what_arg) : 
00092       BelosError(what_arg)
00093     {}
00094   };
00095 
00114   template<class ScalarType, class MV, class OP>
00115   class MinresSolMgr : public SolverManager<ScalarType,MV,OP> {
00116     
00117   private:
00118     typedef MultiVecTraits<ScalarType,MV> MVT;
00119     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00120     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00121     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00122     typedef Teuchos::ScalarTraits< MagnitudeType > MT;
00123     
00124   public:
00125 
00137     static Teuchos::RCP<const Teuchos::ParameterList> defaultParameters();
00138 
00140 
00141 
00150     MinresSolMgr();
00151 
00178     MinresSolMgr (const Teuchos::RCP<LinearProblem< ScalarType, MV, OP> > &problem,
00179                   const Teuchos::RCP<Teuchos::ParameterList> &params);
00180     
00182     virtual ~MinresSolMgr() {};
00184     
00186 
00187 
00189     const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00190       return *problem_;
00191     }
00192 
00194     Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const {
00195       if (defaultParams_.is_null()) {
00196   defaultParams_ = defaultParameters ();
00197       }
00198       return defaultParams_;
00199     }
00200     
00202     Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { 
00203       return params_; 
00204     }
00205     
00215     Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00216       return Teuchos::tuple (timerSolve_);
00217     }
00218 
00224     MagnitudeType achievedTol() const {
00225       return achievedTol_;
00226     }
00227 
00229     int getNumIters() const {
00230       return numIters_;
00231     }
00232 
00238     bool isLOADetected() const { return false; }
00239  
00241     
00243 
00244    
00245     void 
00246     setProblem (const Teuchos::RCP<LinearProblem<ScalarType, MV, OP> > &problem) 
00247     { 
00248       validateProblem (problem);
00249       problem_ = problem; 
00250     }
00251    
00252     void 
00253     setParameters (const Teuchos::RCP<Teuchos::ParameterList>& params);
00254     
00256    
00258 
00259 
00260     void 
00261     reset (const ResetType type) 
00262     { 
00263       if ((type & Belos::Problem) && ! problem_.is_null()) {
00264   problem_->setProblem (); 
00265       }
00266     }
00268  
00270 
00271     
00289     ReturnType solve();
00290     
00292     
00295     
00296     std::string description() const;
00297     
00299     
00300   private:
00302     Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00303     
00305     Teuchos::RCP<OutputManager<ScalarType> > printer_;
00306     Teuchos::RCP<std::ostream> outputStream_;
00307 
00314     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00315 
00319     Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00320 
00324     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
00325 
00329     Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > impConvTest_;
00330 
00334     Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > expConvTest_;
00335 
00340     Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00341 
00345     mutable Teuchos::RCP<const Teuchos::ParameterList> defaultParams_;
00346 
00348     Teuchos::RCP<Teuchos::ParameterList> params_;
00349 
00351     MagnitudeType convtol_;
00352 
00354     MagnitudeType achievedTol_;
00355 
00357     int maxIters_;
00358 
00360     int numIters_;
00361 
00363     int blockSize_;
00364 
00366     int verbosity_;
00367 
00369     int outputStyle_;
00370 
00372     int outputFreq_;
00373     
00375     std::string label_;
00376 
00378     Teuchos::RCP<Teuchos::Time> timerSolve_;
00379 
00381     bool parametersSet_;
00382 
00388     static void
00389     validateProblem (const Teuchos::RCP<LinearProblem<ScalarType, MV, OP> >& problem);
00390   };
00391 
00392 
00393   template<class ScalarType, class MV, class OP>
00394   Teuchos::RCP<const Teuchos::ParameterList> 
00395   MinresSolMgr<ScalarType, MV, OP>::defaultParameters()
00396   {
00397     using Teuchos::ParameterList;
00398     using Teuchos::parameterList;
00399     using Teuchos::RCP;
00400     using Teuchos::rcp;
00401     using Teuchos::rcpFromRef;
00402     using Teuchos::EnhancedNumberValidator;
00403     typedef MagnitudeType MT;
00404     typedef Teuchos::ScalarTraits<MT> MST;
00405     typedef Teuchos::ScalarTraits<int> IST;
00406 
00407     // List of parameters accepted by MINRES, and their default values.
00408     RCP<ParameterList> pl = parameterList ("MINRES");
00409   
00410     pl->set ("Convergence Tolerance", MST::squareroot (MST::eps()),
00411        "Relative residual tolerance that needs to be achieved by "
00412        "the iterative solver, in order for the linear system to be "
00413        "declared converged.",
00414        rcp (new EnhancedNumberValidator<MT> (MST::zero(), MST::rmax())));
00415     pl->set ("Maximum Iterations", static_cast<int>(1000),
00416        "Maximum number of iterations allowed for each right-hand "
00417        "side solved.",
00418        rcp (new EnhancedNumberValidator<int> (0, INT_MAX)));
00419     pl->set ("Block Size", static_cast<int>(1),
00420        "Number of vectors in each block.  WARNING: The current "
00421        "implementation of MINRES only accepts a block size of 1, "
00422        "since it can only solve for 1 right-hand side at a time.",
00423        rcp (new EnhancedNumberValidator<int> (1, 1)));
00424     pl->set ("Verbosity", (int) Belos::Errors,
00425        "The type(s) of solver information that should "
00426        "be written to the output stream.");
00427     pl->set ("Output Style", (int) Belos::General,
00428        "What style is used for the solver information written "
00429        "to the output stream.");
00430     pl->set ("Output Frequency", static_cast<int>(-1),
00431        "How often (in terms of number of iterations) intermediate "
00432        "convergence information should be written to the output stream."
00433        "  -1 means never.");
00434     pl->set ("Output Stream", rcpFromRef(std::cout),
00435        "A reference-counted pointer to the output stream where all "
00436        "solver output is sent.  The output stream defaults to stdout.");
00437     pl->set ("Timer Label", std::string("Belos"),
00438        "The string to use as a prefix for the timer labels.");
00439     return pl;
00440   }
00441 
00442   //
00443   // Empty Constructor
00444   //
00445   template<class ScalarType, class MV, class OP>
00446   MinresSolMgr<ScalarType,MV,OP>::MinresSolMgr () :
00447     numIters_ (0),
00448     parametersSet_ (false)
00449   {}
00450 
00451   //
00452   // Primary constructor (use this one)
00453   //
00454   template<class ScalarType, class MV, class OP>
00455   MinresSolMgr<ScalarType, MV, OP>::
00456   MinresSolMgr (const Teuchos::RCP<LinearProblem<ScalarType, MV, OP> > &problem,
00457     const Teuchos::RCP<Teuchos::ParameterList>& params) :
00458     problem_ (problem),
00459     numIters_ (0),
00460     parametersSet_ (false)
00461   {
00462     TEUCHOS_TEST_FOR_EXCEPTION(problem_.is_null(), std::invalid_argument, 
00463              "MinresSolMgr: The version of the constructor "
00464              "that takes a LinearProblem to solve was given a "
00465              "null LinearProblem.");
00466     setParameters (params);
00467   }
00468 
00469   template<class ScalarType, class MV, class OP>
00470   void
00471   MinresSolMgr<ScalarType, MV, OP>::
00472   validateProblem (const Teuchos::RCP<LinearProblem<ScalarType, MV, OP> >& problem) 
00473   {
00474     TEUCHOS_TEST_FOR_EXCEPTION(problem.is_null(), 
00475       MinresSolMgrLinearProblemFailure,
00476       "MINRES requires that you have provided a nonnull LinearProblem to the "
00477       "solver manager, before you call the solve() method.");
00478     TEUCHOS_TEST_FOR_EXCEPTION(problem->getOperator().is_null(),
00479       MinresSolMgrLinearProblemFailure,
00480       "MINRES requires a LinearProblem object with a non-null operator (the "
00481       "matrix A).");
00482     TEUCHOS_TEST_FOR_EXCEPTION(problem->getRHS().is_null(),
00483       MinresSolMgrLinearProblemFailure,
00484       "MINRES requires a LinearProblem object with a non-null right-hand side.");
00485     TEUCHOS_TEST_FOR_EXCEPTION( ! problem->isProblemSet(),
00486       MinresSolMgrLinearProblemFailure,
00487       "MINRES requires that before you give it a LinearProblem to solve, you "
00488       "must first call the linear problem's setProblem() method.");
00489   }
00490 
00491   template<class ScalarType, class MV, class OP>
00492   void 
00493   MinresSolMgr<ScalarType, MV, OP>::
00494   setParameters (const Teuchos::RCP<Teuchos::ParameterList>& params)
00495   {
00496     using Teuchos::ParameterList;
00497     using Teuchos::parameterList;
00498     using Teuchos::RCP;
00499     using Teuchos::rcp;
00500     using Teuchos::rcpFromRef;
00501     using Teuchos::null;
00502     using Teuchos::is_null;
00503     using std::string;
00504     using std::ostream;
00505     using std::endl;
00506 
00507     RCP<const ParameterList> defaults = getValidParameters ();
00508 
00509     RCP<ParameterList> pl;
00510     if (params.is_null()) {
00511       // We don't need to validate the default parameter values.
00512       pl = parameterList (*defaults);
00513     } else {
00514       pl = params;
00515       pl->validateParametersAndSetDefaults (*defaults);
00516     }
00517 
00518     //
00519     // Read parameters from the parameter list.  We have already
00520     // populated it with defaults.
00521     //
00522     blockSize_ = pl->get<int> ("Block Size");
00523     verbosity_ = pl->get<int> ("Verbosity");
00524     outputStyle_ = pl->get<int> ("Output Style");
00525     outputFreq_ = pl->get<int>("Output Frequency");
00526     outputStream_ = pl->get<RCP<std::ostream> > ("Output Stream");
00527     convtol_ = pl->get<MagnitudeType> ("Convergence Tolerance");
00528     maxIters_ = pl->get<int> ("Maximum Iterations");
00529     //
00530     // All done reading parameters from the parameter list.
00531     // Now we know it's valid and we can store it.
00532     //
00533     params_ = pl;
00534 
00535     // Change the timer label, and create the timer if necessary.
00536     const string newLabel = pl->get<string> ("Timer Label");
00537     {
00538       if (newLabel != label_ || timerSolve_.is_null()) {
00539   label_ = newLabel;
00540 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00541   const string solveLabel = label_ + ": MinresSolMgr total solve time";
00542   // Unregister the old timer before creating a new one.
00543   if (! timerSolve_.is_null()) {
00544     Teuchos::TimeMonitor::clearCounter (label_);
00545     timerSolve_ = Teuchos::null;
00546   }
00547   timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
00548 #endif // BELOS_TEUCHOS_TIME_MONITOR
00549       }
00550     }
00551 
00552     // Create output manager, if necessary; otherwise, set its parameters.
00553     bool recreatedPrinter = false;
00554     if (printer_.is_null()) {
00555       printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
00556       recreatedPrinter = true;
00557     } else {
00558       // Set the output stream's verbosity level.
00559       printer_->setVerbosity (verbosity_);
00560       // Tell the output manager about the new output stream.
00561       printer_->setOStream (outputStream_);
00562     }
00563 
00564     //
00565     // Set up the convergence tests
00566     //
00567     typedef StatusTestGenResNorm<ScalarType, MV, OP> res_norm_type;
00568     typedef StatusTestCombo<ScalarType, MV, OP> combo_type;
00569 
00570     // Do we need to allocate at least one of the implicit or explicit
00571     // residual norm convergence tests?
00572     const bool allocatedConvergenceTests = 
00573       impConvTest_.is_null() || expConvTest_.is_null();
00574 
00575     // Allocate or set the tolerance of the implicit residual norm
00576     // convergence test.
00577     if (impConvTest_.is_null()) {
00578       impConvTest_ = rcp (new res_norm_type (convtol_));
00579       impConvTest_->defineResForm (res_norm_type::Implicit, TwoNorm);
00580       // TODO (mfh 03 Nov 2011) Allow users to define the type of
00581       // scaling (or a custom scaling factor).
00582       impConvTest_->defineScaleForm (NormOfInitRes, TwoNorm);
00583     } else {
00584       impConvTest_->setTolerance (convtol_);
00585     }
00586 
00587     // Allocate or set the tolerance of the explicit residual norm
00588     // convergence test.
00589     if (expConvTest_.is_null()) {
00590       expConvTest_ = rcp (new res_norm_type (convtol_));
00591       expConvTest_->defineResForm (res_norm_type::Explicit, TwoNorm);
00592       // TODO (mfh 03 Nov 2011) Allow users to define the type of
00593       // scaling (or a custom scaling factor).
00594       expConvTest_->defineScaleForm (NormOfInitRes, TwoNorm);
00595     } else {
00596       expConvTest_->setTolerance (convtol_);
00597     }
00598 
00599     // Whether we need to recreate the full status test.  We only need
00600     // to do that if at least one of convTest_ or maxIterTest_ had to
00601     // be reallocated.
00602     bool needToRecreateFullStatusTest = sTest_.is_null();
00603 
00604     // Residual status test is a combo of the implicit and explicit
00605     // convergence tests.
00606     if (convTest_.is_null() || allocatedConvergenceTests) {
00607       convTest_ = rcp (new combo_type (combo_type::SEQ, impConvTest_, expConvTest_));
00608       needToRecreateFullStatusTest = true;
00609     }
00610 
00611     // Maximum number of iterations status test.  It tells the solver to
00612     // stop iteration, if the maximum number of iterations has been
00613     // exceeded.  Initialize it if we haven't yet done so, otherwise
00614     // tell it the new maximum number of iterations.
00615     if (maxIterTest_.is_null()) {
00616       maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
00617       needToRecreateFullStatusTest = true;
00618     } else {
00619       maxIterTest_->setMaxIters (maxIters_);
00620     }
00621 
00622     // Create the full status test if we need to.  
00623     //
00624     // The full status test: the maximum number of iterations have
00625     // been reached, OR the residual has converged.
00626     //
00627     // "If we need to" means either that the status test was never
00628     // created before, or that its two component tests had to be
00629     // reallocated.
00630     if (needToRecreateFullStatusTest) {
00631       sTest_ = rcp (new combo_type (combo_type::OR, maxIterTest_, convTest_));
00632     }
00633   
00634     // If necessary, create the status test output class.  This class
00635     // manages and formats the output from the status test.  We have
00636     // to recreate the output test if we had to (re)allocate either
00637     // printer_ or sTest_.
00638     if (outputTest_.is_null() || needToRecreateFullStatusTest || recreatedPrinter) {
00639       StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
00640       outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_,
00641                Passed+Failed+Undefined);
00642     } else {
00643       outputTest_->setOutputFrequency (outputFreq_);
00644     }
00645     // Set the solver string for the output test.
00646     // StatusTestOutputFactory has no constructor argument for this.
00647     outputTest_->setSolverDesc (std::string (" MINRES "));
00648 
00649     // Inform the solver manager that the current parameters were set.
00650     parametersSet_ = true;
00651 
00652     if (verbosity_ & Debug) {
00653       using std::endl;
00654 
00655       std::ostream& dbg = printer_->stream (Debug);
00656       dbg << "MINRES parameters:" << endl << params_ << endl;
00657     }
00658   }
00659 
00660 
00661   template<class ScalarType, class MV, class OP>
00662   ReturnType MinresSolMgr<ScalarType,MV,OP>::solve() 
00663   {
00664     using Teuchos::RCP;
00665     using Teuchos::rcp;
00666     using Teuchos::rcp_const_cast;
00667     using std::endl;
00668 
00669     if (! parametersSet_) {
00670       setParameters (params_);
00671     } 
00672     std::ostream& dbg = printer_->stream (Debug);
00673 
00674 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00675     Teuchos::TimeMonitor solveTimerMonitor (*timerSolve_);
00676 #endif // BELOS_TEUCHOS_TIME_MONITOR
00677 
00678     // We need a problem to solve, else we can't solve it.
00679     validateProblem (problem_);
00680 
00681     // Reset the status test for this solve.  
00682     outputTest_->reset();
00683 
00684     // The linear problem has this many right-hand sides to solve.
00685     // MINRES can solve only one at a time, so we solve for each
00686     // right-hand side in succession.
00687     const int numRHS2Solve = MVT::GetNumberVecs (*(problem_->getRHS()));
00688 
00689     // Create MINRES iteration object.  Pass along the solver
00690     // manager's parameters, which have already been validated.
00691     typedef MinresIter<ScalarType, MV, OP> iter_type;
00692     RCP<iter_type> minres_iter =
00693       rcp (new iter_type (problem_, printer_, outputTest_, *params_));
00694 
00695     // The index/indices of the right-hand sides for which MINRES did
00696     // _not_ converge.  Hopefully this is empty after the for loop
00697     // below!  If it is not empty, at least one right-hand side did
00698     // not converge.
00699     std::vector<int> notConverged;
00700     std::vector<int> currentIndices(1);
00701 
00702     numIters_ = 0;
00703 
00704     // Solve for each right-hand side in turn.
00705     for (int currentRHS = 0; currentRHS < numRHS2Solve; ++currentRHS) {
00706       // Inform the linear problem of the right-hand side(s) currently
00707       // being solved.  MINRES only knows how to solve linear problems
00708       // with one right-hand side, so we only include one index, which
00709       // is the index of the current right-hand side.
00710       currentIndices[0] = currentRHS;
00711       problem_->setLSIndex (currentIndices);
00712 
00713       dbg << "-- Current right-hand side index being solved: " 
00714     << currentRHS << endl;
00715 
00716       // Reset the number of iterations.
00717       minres_iter->resetNumIters();
00718       // Reset the number of calls that the status test output knows about.
00719       outputTest_->resetNumCalls();
00720       // Set the new state and initialize the solver.
00721       MinresIterationState<ScalarType, MV> newstate;
00722 
00723       // Get the residual vector for the current linear system 
00724       // (that is, for the current right-hand side).
00725       newstate.Y = MVT::CloneViewNonConst (*(rcp_const_cast<MV> (problem_->getInitResVec())), currentIndices);
00726       minres_iter->initializeMinres (newstate);
00727 
00728       // Attempt to solve for the solution corresponding to the
00729       // current right-hand side.
00730       while (true) {
00731   try {
00732     minres_iter->iterate();
00733 
00734     // First check for convergence
00735     if (convTest_->getStatus() == Passed) {
00736       dbg << "---- Converged after " << maxIterTest_->getNumIters() 
00737     << " iterations" << endl;
00738       break;
00739     }
00740     // Now check for max # of iterations
00741     else if (maxIterTest_->getStatus() == Passed) {
00742       dbg << "---- Did not converge after " << maxIterTest_->getNumIters() 
00743     << " iterations" << endl;
00744       // This right-hand side didn't converge!
00745       notConverged.push_back (currentRHS);
00746       break;
00747     } else {
00748       // If we get here, we returned from iterate(), but none of
00749       // our status tests Passed.  Something is wrong, and it is
00750       // probably our fault.
00751       TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
00752               "Belos::MinresSolMgr::solve(): iterations neither converged, "
00753         "nor reached the maximum number of iterations " << maxIters_ 
00754         << ".  That means something went wrong.");
00755     }
00756   } catch (const std::exception &e) {
00757     printer_->stream (Errors) 
00758       << "Error! Caught std::exception in MinresIter::iterate() at "
00759       << "iteration " << minres_iter->getNumIters() << endl
00760       << e.what() << endl;
00761     throw e;
00762   }
00763       }
00764   
00765       // Inform the linear problem that we are finished with the
00766       // current right-hand side.  It may or may not have converged,
00767       // but we don't try again if the first time didn't work.
00768       problem_->setCurrLS();
00769     
00770       // Get iteration information for this solve: total number of
00771       // iterations for all right-hand sides.
00772       numIters_ += maxIterTest_->getNumIters();
00773     }
00774     
00775     // Print final summary of the solution process
00776     sTest_->print (printer_->stream (FinalSummary));
00777 
00778     // Print timing information, if the corresponding compile-time and
00779     // run-time options are enabled.
00780 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00781     // Calling summarize() can be expensive, so don't call unless the
00782     // user wants to print out timing details.  summarize() will do all
00783     // the work even if it's passed a "black hole" output stream.
00784     if (verbosity_ & TimingDetails) {
00785       Teuchos::TimeMonitor::summarize (printer_->stream (TimingDetails));
00786     }
00787 #endif // BELOS_TEUCHOS_TIME_MONITOR
00788 
00789     // Save the convergence test value ("achieved tolerance") for this
00790     // solve.  This solver always has two residual norm status tests:
00791     // an explicit and an implicit test.  The master convergence test
00792     // convTest_ is a SEQ combo of the implicit resp. explicit tests.
00793     // If the implicit test never passes, then the explicit test won't
00794     // ever be executed.  This manifests as
00795     // expConvTest_->getTestValue()->size() < 1.  We deal with this
00796     // case by using the values returned by
00797     // impConvTest_->getTestValue().
00798     {
00799       const std::vector<MagnitudeType>* pTestValues = expConvTest_->getTestValue();
00800       if (pTestValues == NULL || pTestValues->size() < 1) {
00801   pTestValues = impConvTest_->getTestValue();
00802       }
00803       TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
00804         "Belos::MinresSolMgr::solve(): The implicit convergence test's getTestValue() "
00805   "method returned NULL.  Please report this bug to the Belos developers.");
00806       TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
00807         "Belos::MinresSolMgr::solve(): The implicit convergence test's getTestValue() "
00808         "method returned a vector of length zero.  Please report this bug to the "
00809         "Belos developers.");
00810 
00811       // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
00812       // achieved tolerances for all vectors in the current solve(), or
00813       // just for the vectors from the last deflation?
00814       achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
00815     }
00816  
00817     if (notConverged.size() > 0) {
00818       return Unconverged;
00819     } else {
00820       return Converged;
00821     }
00822   }
00823 
00824   //  This method requires the solver manager to return a std::string that describes itself.
00825   template<class ScalarType, class MV, class OP>
00826   std::string MinresSolMgr<ScalarType,MV,OP>::description() const
00827   {
00828     std::ostringstream oss;
00829     oss << "Belos::MinresSolMgr< " 
00830   << Teuchos::ScalarTraits<ScalarType>::name()
00831   <<", MV, OP >";
00832     // oss << "{";
00833     // oss << "Block Size=" << blockSize_;
00834     // oss << "}";
00835     return oss.str();
00836   }
00837   
00838 } // end Belos namespace
00839 
00840 #endif /* BELOS_MINRES_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines