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