Belos Version of the Day
BelosInnerSolver.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_InnerSolver_hpp
00043 #define __Belos_InnerSolver_hpp
00044 
00045 #include <BelosInnerSolveResult.hpp>
00046 #include <BelosLinearProblem.hpp>
00047 #include <BelosMultiVecTraits.hpp>
00048 #include <BelosOperatorTraits.hpp>
00049 #include <Teuchos_ScalarTraits.hpp>
00050 
00051 namespace Belos {
00052 
00064   template<class Scalar, class MV, class OP>
00065   Teuchos::RCP<LinearProblem<Scalar, MV, OP> >
00066   problemWithNewRHS (const Teuchos::RCP<const LinearProblem<Scalar, MV, OP> >& problem,
00067          const Teuchos::RCP<const MV>& B)
00068   {
00069     using Teuchos::is_null;
00070     using Teuchos::nonnull;
00071     using Teuchos::null;
00072     using Teuchos::RCP;
00073     using Teuchos::rcp;
00074     typedef LinearProblem<Scalar, MV, OP> lp_type;
00075     typedef MultiVecTraits<Scalar, MV> MVT;
00076 
00077     RCP<const OP> A = problem->getOperator (); 
00078     RCP<MV> X_orig = problem->getLHS ();
00079     TEST_FOR_EXCEPTION(is_null (X_orig), std::invalid_argument,
00080            "problemWithNewRHS(): The original LinearProblem's "
00081            "initial guess / current approximate solution (getLHS())"
00082            " is null.  We need an initial guess or current approxim"
00083            "ate solution in order to know the domain of the (right-"
00084            "preconditioned, if applicable) operator.  This is "
00085            "because Belos::MultiVecTraits does not include the idea"
00086            " of the domain and range of an operator, or the space "
00087            "to which a vector belongs.");
00088     TEST_FOR_EXCEPTION(is_null (B), std::invalid_argument,
00089            "problemWithNewRHS(): the given new right-hand side B "
00090            "is null.");
00091     RCP<MV> X = MVT::CloneCopy (problem->getLHS ());
00092 
00093     RCP<lp_type> lp (new lp_type (A, X, B));
00094     lp->setLeftPrec (problem->getLeftPrec ());
00095     lp->setRightPrec (problem->getRightPrec ());
00096     // Compute initial residual(s) and prepare the problem for solution.
00097     lp->setProblem ();
00098     return lp;
00099   }
00100 
00101 
00134   template<class Scalar, class MV, class OP>
00135   class InnerSolver {
00136   public:
00137     typedef Scalar scalar_type;
00138     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00139     typedef MV multivector_type;
00140     typedef OP operator_type;
00141 
00143     virtual ~InnerSolver() {}
00144 
00151     virtual Teuchos::RCP<const Teuchos::ParameterList> 
00152     getCurrentParameters() const = 0;
00153 
00197     virtual InnerSolveResult
00198     solve (const Teuchos::RCP<MV>& X,
00199      const Teuchos::RCP<const MV>& B,
00200      const magnitude_type convTol,
00201      const int maxItersPerRestart,
00202      const int maxNumRestarts) = 0;
00203 
00247     virtual InnerSolveResult
00248     solve (const Teuchos::RCP<MV>& X,
00249      const Teuchos::RCP<const MV>& B) = 0;
00250   };
00251 
00252 
00261   template<class Scalar, class MV, class OP>
00262   class OperatorTraits<Scalar, MV, InnerSolver<Scalar, MV, OP> > {
00263   public:
00264     static void
00265     Apply (const InnerSolver<Scalar, MV, OP>& Op,
00266      const MV& x,
00267      MV& y,
00268      ETrans trans = NOTRANS)
00269     {
00270       using Teuchos::RCP;
00271       using Teuchos::rcpFromRef;
00272 
00273       TEST_FOR_EXCEPTION(trans != NOTRANS, std::invalid_argument,
00274        "Belos::InnerSolver is not able to solve the "
00275        "transposed system.");
00276       RCP<const MV> x_ptr = rcpFromRef (x);
00277       RCP<MV> y_ptr = rcpFromRef (y);
00278       (void) Op.solve (y_ptr, x_ptr);
00279     }
00280 
00281   };
00282 
00290   template<class Scalar, class MV, class OP>
00291   class UndefinedWrapperType {
00292   public:
00293     UndefinedWrapperType() {
00294       (void) Scalar::this_specialization_is_not_defined();
00295       (void) MV::this_specialization_is_not_defined();
00296       (void) OP::this_specialization_is_not_defined();
00297     }
00298   };
00299 
00319   template<class Scalar, class MV, class OP>
00320   class InnerSolverTraits {
00321   public:
00322     typedef Scalar scalar_type;
00323     typedef MV multivector_type;
00324     typedef OP operator_type;
00325     typedef InnerSolver<scalar_type, multivector_type, operator_type> inner_solver_type;
00326     typedef UndefinedWrapperType<Scalar, MV, OP> wrapper_type;
00327 
00332     static Teuchos::RCP<OP>
00333     makeInnerSolverOperator (const Teuchos::RCP<InnerSolver<Scalar, MV, OP> >& solver)
00334     {
00335       using Teuchos::rcp;
00336       using Teuchos::rcp_implicit_cast;
00337       // If this class is not specialized for the given combination of
00338       // (Scalar, MV, OP), the constructor of wrapper_type here will
00339       // (deliberately) raise a compiler error.
00340       return rcp_implicit_cast<operator_type> (rcp (new wrapper_type (solver)));
00341     }
00342 
00353     static Teuchos::RCP<inner_solver_type>
00354     getInnerSolver (const Teuchos::RCP<operator_type>& op)
00355     {
00356       using Teuchos::RCP;
00357       using Teuchos::rcp_dynamic_cast;
00358       // If this class is not specialized for the given combination of
00359       // (Scalar, MV, OP), the instantiation of the wrapper_type class
00360       // here will (deliberately) raise a compiler error.
00361       RCP<wrapper_type> wrapper = rcp_dynamic_cast<wrapper_type> (op, true);
00362       return wrapper->getInnerSolver();
00363     }
00364   };
00365 
00366 } // namespace Belos
00367 
00368 #endif // __Belos_InnerSolver_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines