Belos Package Browser (Single Doxygen Collection)  Development
BelosInnerSolver.hpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
4 // Belos: Block Linear Solvers Package
5 // Copyright 2004 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
42 #ifndef __Belos_InnerSolver_hpp
43 #define __Belos_InnerSolver_hpp
44 
46 #include <BelosLinearProblem.hpp>
47 #include <BelosMultiVecTraits.hpp>
48 #include <BelosOperatorTraits.hpp>
49 #include <Teuchos_ScalarTraits.hpp>
50 
51 namespace Belos {
52 
64  template<class Scalar, class MV, class OP>
65  Teuchos::RCP<LinearProblem<Scalar, MV, OP> >
66  problemWithNewRHS (const Teuchos::RCP<const LinearProblem<Scalar, MV, OP> >& problem,
67  const Teuchos::RCP<const MV>& B)
68  {
69  using Teuchos::is_null;
70  using Teuchos::nonnull;
71  using Teuchos::null;
72  using Teuchos::RCP;
73  using Teuchos::rcp;
74  typedef LinearProblem<Scalar, MV, OP> lp_type;
75  typedef MultiVecTraits<Scalar, MV> MVT;
76 
77  RCP<const OP> A = problem->getOperator ();
78  RCP<MV> X_orig = problem->getLHS ();
79  TEUCHOS_TEST_FOR_EXCEPTION(is_null (X_orig), std::invalid_argument,
80  "problemWithNewRHS(): The original LinearProblem's "
81  "initial guess / current approximate solution (getLHS())"
82  " is null. We need an initial guess or current approxim"
83  "ate solution in order to know the domain of the (right-"
84  "preconditioned, if applicable) operator. This is "
85  "because Belos::MultiVecTraits does not include the idea"
86  " of the domain and range of an operator, or the space "
87  "to which a vector belongs.");
88  TEUCHOS_TEST_FOR_EXCEPTION(is_null (B), std::invalid_argument,
89  "problemWithNewRHS(): the given new right-hand side B "
90  "is null.");
91  RCP<MV> X = MVT::CloneCopy (problem->getLHS ());
92 
93  RCP<lp_type> lp (new lp_type (A, X, B));
94  lp->setLeftPrec (problem->getLeftPrec ());
95  lp->setRightPrec (problem->getRightPrec ());
96  // Compute initial residual(s) and prepare the problem for solution.
97  lp->setProblem ();
98  return lp;
99  }
100 
101 
137  template<class Scalar, class MV, class OP>
138  class InnerSolver {
139  public:
140  typedef Scalar scalar_type;
141  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
142  typedef MV multivector_type;
143  typedef OP operator_type;
144 
146  virtual ~InnerSolver() {}
147 
154  virtual Teuchos::RCP<const Teuchos::ParameterList>
155  getCurrentParameters() const = 0;
156 
200  virtual InnerSolveResult
201  solve (const Teuchos::RCP<MV>& X,
202  const Teuchos::RCP<const MV>& B,
203  const magnitude_type convTol,
204  const int maxItersPerRestart,
205  const int maxNumRestarts) = 0;
206 
250  virtual InnerSolveResult
251  solve (const Teuchos::RCP<MV>& X,
252  const Teuchos::RCP<const MV>& B) = 0;
253  };
254 
255 
264  template<class Scalar, class MV, class OP>
265  class OperatorTraits<Scalar, MV, InnerSolver<Scalar, MV, OP> > {
266  public:
267  static void
269  const MV& x,
270  MV& y,
271  ETrans trans = NOTRANS)
272  {
273  using Teuchos::RCP;
274  using Teuchos::rcpFromRef;
275 
276  TEUCHOS_TEST_FOR_EXCEPTION(trans != NOTRANS, std::invalid_argument,
277  "Belos::InnerSolver is not able to solve the "
278  "transposed system.");
279  RCP<const MV> x_ptr = rcpFromRef (x);
280  RCP<MV> y_ptr = rcpFromRef (y);
281  (void) Op.solve (y_ptr, x_ptr);
282  }
283 
284  };
285 
293  template<class Scalar, class MV, class OP>
295  public:
297  (void) Scalar::this_specialization_is_not_defined();
298  (void) MV::this_specialization_is_not_defined();
299  (void) OP::this_specialization_is_not_defined();
300  }
301  };
302 
322  template<class Scalar, class MV, class OP>
324  public:
325  typedef Scalar scalar_type;
326  typedef MV multivector_type;
327  typedef OP operator_type;
330 
335  static Teuchos::RCP<OP>
337  {
338  using Teuchos::rcp;
339  using Teuchos::rcp_implicit_cast;
340  // If this class is not specialized for the given combination of
341  // (Scalar, MV, OP), the constructor of wrapper_type here will
342  // (deliberately) raise a compiler error.
343  return rcp_implicit_cast<operator_type> (rcp (new wrapper_type (solver)));
344  }
345 
356  static Teuchos::RCP<inner_solver_type>
357  getInnerSolver (const Teuchos::RCP<operator_type>& op)
358  {
359  using Teuchos::RCP;
360  using Teuchos::rcp_dynamic_cast;
361  // If this class is not specialized for the given combination of
362  // (Scalar, MV, OP), the instantiation of the wrapper_type class
363  // here will (deliberately) raise a compiler error.
364  RCP<wrapper_type> wrapper = rcp_dynamic_cast<wrapper_type> (op, true);
365  return wrapper->getInnerSolver();
366  }
367  };
368 
369 } // namespace Belos
370 
371 #endif // __Belos_InnerSolver_hpp
virtual Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const =0
Current parameters for the inner solver implementation.
Teuchos::RCP< LinearProblem< Scalar, MV, OP > > problemWithNewRHS(const Teuchos::RCP< const LinearProblem< Scalar, MV, OP > > &problem, const Teuchos::RCP< const MV > &B)
New LinearProblem with different right-hand side.
static Teuchos::RCP< inner_solver_type > getInnerSolver(const Teuchos::RCP< operator_type > &op)
Return the given wrapper&#39;s inner solver object.
static Teuchos::RCP< OP > makeInnerSolverOperator(const Teuchos::RCP< InnerSolver< Scalar, MV, OP > > &solver)
Wrap the given inner solver in a wrapper_type.
Undefined wrapper type, to check at compile time whether InnerSolverTraits has been specialized...
Represents the result of an inner solve.
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
ETrans
Whether to apply the (conjugate) transpose of an operator.
Definition: BelosTypes.hpp:81
Traits class which defines basic operations on multivectors.
virtual InnerSolveResult solve(const Teuchos::RCP< MV > &X, const Teuchos::RCP< const MV > &B, const magnitude_type convTol, const int maxItersPerRestart, const int maxNumRestarts)=0
Solve for the given right-hand side(s) B.
Wrap an InnerSolver in an OP (operator).
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitude_type
UndefinedWrapperType< Scalar, MV, OP > wrapper_type
virtual ~InnerSolver()
Virtual destructor, for correctness.
static void Apply(const InnerSolver< Scalar, MV, OP > &Op, const MV &x, MV &y, ETrans trans=NOTRANS)
Inner solver interface.
InnerSolver< scalar_type, multivector_type, operator_type > inner_solver_type
Class which defines basic traits for the operator type.