Belos Package Browser (Single Doxygen Collection) Development
BelosGmresBaseIteration.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_GmresBaseIteration_hpp
00043 #define __Belos_GmresBaseIteration_hpp
00044 
00048 
00049 #include "BelosIteration.hpp"
00050 #include "BelosGmresBaseFactory.hpp"
00051 
00052 namespace Belos {
00053 
00064   template<class Scalar, class MV, class OP>
00065   class GmresBaseIteration : public Iteration<Scalar, MV, OP> {
00066   private:
00069     typedef GmresBase<Scalar, MV, OP> impl_type;
00070     
00071   public:
00072     typedef Scalar scalar_type;
00073     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00074     typedef MV multivector_type;
00075     typedef OP operator_type;
00076 
00096     void 
00097     iterate() 
00098     {
00099       // GmresBaseIteration is an Iteration subclass, and the status
00100       // checks (should) only call Iteration methods to get the
00101       // information they want.  We add another "status test," which
00102       // ensures that GMRES doesn't exceed the maximum number of
00103       // iterations allowed (for storage of the basis vectors).  The
00104       // status test will also perform intermediate output via the
00105       // OutputManager object.
00106       while (stest_->checkStatus(this) != Passed && impl_->canAdvance()) {
00107   impl_->advance();
00108       }
00109     }
00110 
00116     void initialize();
00117 
00127     int getNumIters() const { return impl_->getNumIters() + 1; }
00128     
00141     void resetNumIters (int iter = 0) { 
00142       TEUCHOS_TEST_FOR_EXCEPTION(iter < 0, std::invalid_argument, "Belos::GmresBaseIter"
00143        "ation::resetNumIters: iter = " << iter << " is "
00144        "invalid; only nonnegative values are allowed.");
00145       impl_->backOut (iter); 
00146     }
00147 
00149     void restart () { 
00150       TEUCHOS_TEST_FOR_EXCEPTION(impl_.is_null(), std::logic_error, 
00151        "Belos::GmresBaseIteration::restart(): GmresBase "
00152        "subclass instance is null.");
00153       impl_->restart ();
00154     }
00155     
00182     Teuchos::RCP<const MV> 
00183     getNativeResiduals (std::vector<magnitude_type>* norms) const
00184     {
00185       using Teuchos::RCP;
00186 
00187       RCP<const MV> nativeResVec = impl_->currentNativeResidualVector ();
00188       const magnitude_type nativeResNorm = impl_->currentNativeResidualNorm ();
00189       if (norms != NULL)
00190   {
00191     // This should always be 1
00192     const int blockSize = getBlockSize();
00193     TEUCHOS_TEST_FOR_EXCEPTION(blockSize != 1, std::logic_error,
00194            "This implementation of Arnoldi/GMRES only "
00195            "supports a block size of 1, but the current "
00196            "block size is " << blockSize << ".");
00197     // Silence the compiler error about comparing signed and unsigned.
00198     if (norms->size() < (unsigned int) blockSize)
00199       norms->resize (blockSize);
00200     for (int k = 0; k < blockSize; ++k)
00201       (*norms)[k] = nativeResNorm;
00202   }
00203       return nativeResVec;
00204     }
00205 
00214     Teuchos::RCP<MV> getCurrentUpdate() const { 
00215       return impl_->getCurrentUpdate();
00216     }
00217 
00219     const LinearProblem<Scalar, MV, OP>& getProblem() const { return *lp_; }
00220 
00225     int getBlockSize() const { return 1; }
00226 
00238     void setBlockSize (int blockSize) {
00239       TEUCHOS_TEST_FOR_EXCEPTION(blockSize != 1, std::invalid_argument,
00240        "Belos::GmresBaseIteration::setBlockSize: blockSize = " 
00241        << blockSize << " is invalid; only blockSize = 1 is "
00242        "allowed for this iteration.");
00243     }
00244 
00246     bool isInitialized() {
00247       return ! impl_.is_null();
00248     }
00249 
00266     GmresBaseIteration (const Teuchos::RCP<LinearProblem<Scalar, MV, OP> >& problem,
00267       const Teuchos::RCP<const OrthoManager<Scalar, MV> >& ortho,
00268       const Teuchos::RCP<OutputManager<Scalar> >& printer,
00269       const Teuchos::RCP<StatusTest<Scalar, MV, OP> >& tester,
00270       const Teuchos::RCP<const Teuchos::ParameterList >& params);
00271     
00273     virtual ~GmresBaseIteration() {}
00274 
00281     void
00282     updateSolution ()
00283     {
00284       impl_->updateSolution ();
00285     }
00286 
00287   private:
00296     mutable Teuchos::RCP<GmresBase<Scalar, MV, OP> > impl_;
00297 
00299     //
00304     Teuchos::RCP<LinearProblem<Scalar, MV, OP> > lp_;
00309     Teuchos::RCP<const OrthoManager<Scalar, MV> > ortho_;
00311     Teuchos::RCP<OutputManager<Scalar> > outMan_;
00313     Teuchos::RCP<StatusTest<Scalar, MV, OP> > stest_;
00315     Teuchos::RCP<const Teuchos::ParameterList> params_;
00316 
00321     static Teuchos::RCP<LinearProblem<Scalar, MV, OP> >
00322     validatedProblem (const Teuchos::RCP<LinearProblem<Scalar, MV, OP> >& problem);
00323   };
00324 
00325   template<class Scalar, class MV, class OP>
00326   GmresBaseIteration<Scalar,MV,OP>::
00327   GmresBaseIteration (const Teuchos::RCP<LinearProblem<Scalar, MV, OP> >& problem,
00328           const Teuchos::RCP<const OrthoManager<Scalar, MV> >& ortho,
00329           const Teuchos::RCP<OutputManager<Scalar> >& printer,
00330           const Teuchos::RCP<StatusTest<Scalar, MV, OP> >& tester,
00331           const Teuchos::RCP<const Teuchos::ParameterList>& params) :
00332     lp_ (validatedProblem (problem)),
00333     ortho_ (ortho), 
00334     outMan_ (printer), 
00335     stest_ (tester), 
00336     params_ (params)
00337   {
00338     initialize();
00339   }
00340 
00341   template<class Scalar, class MV, class OP>
00342   void
00343   GmresBaseIteration<Scalar,MV,OP>::initialize ()
00344   {
00345     if (! isInitialized())
00346       {
00347   typedef GmresBaseFactory<Scalar, MV, OP> factory_type;
00348   impl_ = factory_type::create (lp_, ortho_, outMan_, params_);
00349       }
00350   }
00351 
00352   template<class Scalar, class MV, class OP>
00353   Teuchos::RCP<LinearProblem<Scalar, MV, OP> >
00354   GmresBaseIteration<Scalar,MV,OP>::
00355   validatedProblem (const Teuchos::RCP<LinearProblem<Scalar, MV, OP> >& problem)
00356   {
00357     const char prefix[] = "Belos::GmresBaseIteration constructor: ";
00358     TEUCHOS_TEST_FOR_EXCEPTION(problem.is_null(), std::invalid_argument,
00359            prefix << "The linear problem (Belos::LinearProblem "
00360            "instance) that you want me to solve is null.");
00361     if (! problem->isProblemSet())
00362       {
00363   const bool notInitialized = problem->getOperator().is_null() || 
00364     problem->getRHS().is_null() || problem->getLHS().is_null();
00365   TEUCHOS_TEST_FOR_EXCEPTION(notInitialized, std::invalid_argument,
00366          prefix << "The given linear problem (Belos::Linear"
00367          "Problem) instance is not fully initialized: "
00368          << (problem->getOperator().is_null() ? "the operator A is null; " : "")
00369          << (problem->getRHS().is_null() ? "the right-hand side B is null; " : "")
00370          << (problem->getLHS().is_null() ? "the initial guess X is null" : "")
00371          << ".");
00372   if (! problem->isProblemSet())
00373     problem->setProblem();
00374   TEUCHOS_TEST_FOR_EXCEPTION(! problem->isProblemSet(), std::logic_error,
00375          prefix << "Although the given LinearProblem instance "
00376          "has non-null operator (matrix A), right-hand side B,"
00377          " and initial guess X, and although its setProblem() "
00378          "method has been called, its isProblemSet() method "
00379          "returns false.  This should not happen.");
00380       }
00381     return problem;
00382   }
00383 
00384 } // namespace Belos
00385 
00386 #endif // __Belos_GmresBaseIteration_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines