Belos Version of the Day
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 
00061   template<class Scalar, class MV, class OP>
00062   class GmresBaseIteration : public Iteration<Scalar, MV, OP> {
00063   private:
00066     typedef GmresBase<Scalar, MV, OP> impl_type;
00067     
00068   public:
00069     typedef Scalar scalar_type;
00070     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00071     typedef MV multivector_type;
00072     typedef OP operator_type;
00073 
00093     void 
00094     iterate() 
00095     {
00096       // GmresBaseIteration is an Iteration subclass, and the status
00097       // checks (should) only call Iteration methods to get the
00098       // information they want.  We add another "status test," which
00099       // ensures that GMRES doesn't exceed the maximum number of
00100       // iterations allowed (for storage of the basis vectors).  The
00101       // status test will also perform intermediate output via the
00102       // OutputManager object.
00103       //
00104       // FIXME (mfh 14 Jan 2011): The owning solver manager object
00105       // really should set up this Iteration so that the status test
00106       // doesn't allow the GmresBase subclass to advance beyond its
00107       // capacity to advance.  For example, the status test should set
00108       // the maximum number of iterations no bigger than the amount of
00109       // vector storage that the GmresBase subclass allocates.  Thus,
00110       // the additional "status test" should be unnecessary.
00111       while (stest_->checkStatus(this) != Passed)
00112   { 
00113     if (! impl_->canAdvance())
00114       throw std::logic_error("Attempt to advance iteration when it is "
00115            "no longer possible to advance!");
00116     impl_->advance();
00117   }
00118       
00119     }
00120 
00126     void initialize();
00127 
00137     int getNumIters() const { return impl_->getNumIters() + 1; }
00138     
00151     void resetNumIters (int iter = 0) { 
00152       TEST_FOR_EXCEPTION(iter < 0, std::invalid_argument, "Belos::GmresBaseIter"
00153        "ation::resetNumIters: iter = " << iter << " is "
00154        "invalid; only nonnegative values are allowed.");
00155       impl_->backOut (iter); 
00156     }
00157 
00159     void restart () { 
00160       TEST_FOR_EXCEPTION(impl_.is_null(), std::logic_error, 
00161        "Belos::GmresBaseIteration::restart(): GmresBase "
00162        "subclass instance is null.");
00163       impl_->restart ();
00164     }
00165     
00192     Teuchos::RCP<const MV> 
00193     getNativeResiduals (std::vector<magnitude_type>* norms) const
00194     {
00195       using Teuchos::RCP;
00196 
00197       RCP<const MV> nativeResVec = impl_->currentNativeResidualVector ();
00198       const magnitude_type nativeResNorm = impl_->currentNativeResidualNorm ();
00199       if (norms != NULL)
00200   {
00201     // This should always be 1
00202     const int blockSize = getBlockSize();
00203     TEST_FOR_EXCEPTION(blockSize != 1, std::logic_error,
00204            "This implementation of Arnoldi/GMRES only "
00205            "supports a block size of 1, but the current "
00206            "block size is " << blockSize << ".");
00207     // Silence the compiler error about comparing signed and unsigned.
00208     if (norms->size() < (unsigned int) blockSize)
00209       norms->resize (blockSize);
00210     for (int k = 0; k < blockSize; ++k)
00211       (*norms)[k] = nativeResNorm;
00212   }
00213       return nativeResVec;
00214     }
00215 
00224     Teuchos::RCP<MV> getCurrentUpdate() const { 
00225       return impl_->getCurrentUpdate();
00226     }
00227 
00229     const LinearProblem<Scalar, MV, OP>& getProblem() const { return *lp_; }
00230 
00235     int getBlockSize() const { return 1; }
00236 
00248     void setBlockSize (int blockSize) {
00249       TEST_FOR_EXCEPTION(blockSize != 1, std::invalid_argument,
00250        "Belos::GmresBaseIteration::setBlockSize: blockSize = " 
00251        << blockSize << " is invalid; only blockSize = 1 is "
00252        "allowed for this iteration.");
00253     }
00254 
00256     bool isInitialized() {
00257       return ! impl_.is_null();
00258     }
00259 
00276     GmresBaseIteration (const Teuchos::RCP<LinearProblem<Scalar, MV, OP> >& problem,
00277       const Teuchos::RCP<const OrthoManager<Scalar, MV> >& ortho,
00278       const Teuchos::RCP<OutputManager<Scalar> >& printer,
00279       const Teuchos::RCP<StatusTest<Scalar, MV, OP> >& tester,
00280       const Teuchos::RCP<const Teuchos::ParameterList >& params);
00281     
00283     virtual ~GmresBaseIteration() {}
00284 
00291     void
00292     updateSolution ()
00293     {
00294       impl_->updateSolution ();
00295     }
00296 
00297   private:
00306     mutable Teuchos::RCP<GmresBase<Scalar, MV, OP> > impl_;
00307 
00309     //
00314     Teuchos::RCP<LinearProblem<Scalar, MV, OP> > lp_;
00319     Teuchos::RCP<const OrthoManager<Scalar, MV> > ortho_;
00321     Teuchos::RCP<OutputManager<Scalar> > outMan_;
00323     Teuchos::RCP<StatusTest<Scalar, MV, OP> > stest_;
00325     Teuchos::RCP<const Teuchos::ParameterList> params_;
00326 
00331     static Teuchos::RCP<LinearProblem<Scalar, MV, OP> >
00332     validatedProblem (const Teuchos::RCP<LinearProblem<Scalar, MV, OP> >& problem);
00333   };
00334 
00335   template<class Scalar, class MV, class OP>
00336   GmresBaseIteration<Scalar,MV,OP>::
00337   GmresBaseIteration (const Teuchos::RCP<LinearProblem<Scalar, MV, OP> >& problem,
00338           const Teuchos::RCP<const OrthoManager<Scalar, MV> >& ortho,
00339           const Teuchos::RCP<OutputManager<Scalar> >& printer,
00340           const Teuchos::RCP<StatusTest<Scalar, MV, OP> >& tester,
00341           const Teuchos::RCP<const Teuchos::ParameterList>& params) :
00342     lp_ (validatedProblem (problem)),
00343     ortho_ (ortho), 
00344     outMan_ (printer), 
00345     stest_ (tester), 
00346     params_ (params)
00347   {
00348     initialize();
00349   }
00350 
00351   template<class Scalar, class MV, class OP>
00352   void
00353   GmresBaseIteration<Scalar,MV,OP>::initialize ()
00354   {
00355     if (! isInitialized())
00356       {
00357   typedef GmresBaseFactory<Scalar, MV, OP> factory_type;
00358   impl_ = factory_type::create (lp_, ortho_, outMan_, params_);
00359       }
00360   }
00361 
00362   template<class Scalar, class MV, class OP>
00363   Teuchos::RCP<LinearProblem<Scalar, MV, OP> >
00364   GmresBaseIteration<Scalar,MV,OP>::
00365   validatedProblem (const Teuchos::RCP<LinearProblem<Scalar, MV, OP> >& problem)
00366   {
00367     const char prefix[] = "Belos::GmresBaseIteration constructor: ";
00368     TEST_FOR_EXCEPTION(problem.is_null(), std::invalid_argument,
00369            prefix << "The linear problem (Belos::LinearProblem "
00370            "instance) that you want me to solve is null.");
00371     if (! problem->isProblemSet())
00372       {
00373   const bool notInitialized = problem->getOperator().is_null() || 
00374     problem->getRHS().is_null() || problem->getLHS().is_null();
00375   TEST_FOR_EXCEPTION(notInitialized, std::invalid_argument,
00376          prefix << "The given linear problem (Belos::Linear"
00377          "Problem) instance is not fully initialized: "
00378          << (problem->getOperator().is_null() ? "the operator A is null; " : "")
00379          << (problem->getRHS().is_null() ? "the right-hand side B is null; " : "")
00380          << (problem->getLHS().is_null() ? "the initial guess X is null" : "")
00381          << ".");
00382   if (! problem->isProblemSet())
00383     problem->setProblem();
00384   TEST_FOR_EXCEPTION(! problem->isProblemSet(), std::logic_error,
00385          prefix << "Although the given LinearProblem instance "
00386          "has non-null operator (matrix A), right-hand side B,"
00387          " and initial guess X, and although its setProblem() "
00388          "method has been called, its isProblemSet() method "
00389          "returns false.  This should not happen.");
00390       }
00391     return problem;
00392   }
00393 
00394 } // namespace Belos
00395 
00396 #endif // __Belos_GmresBaseIteration_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines