Thyra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Thyra_DefaultInverseLinearOp_def.hpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Thyra: Interfaces and Support for Abstract Numerical Algorithms
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
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 Roscoe A. Bartlett (bartlettra@ornl.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #ifndef THYRA_DEFAULT_INVERSE_LINEAR_OP_DEF_HPP
43 #define THYRA_DEFAULT_INVERSE_LINEAR_OP_DEF_HPP
44 
45 #include "Thyra_DefaultInverseLinearOp_decl.hpp"
46 #include "Thyra_MultiVectorStdOps.hpp"
47 #include "Thyra_AssertOp.hpp"
48 #include "Teuchos_Utils.hpp"
49 #include "Teuchos_TypeNameTraits.hpp"
50 
51 
52 namespace Thyra {
53 
54 
55 // Constructors/initializers/accessors
56 
57 
58 template<class Scalar>
60 {}
61 
62 
63 template<class Scalar>
66  const SolveCriteria<Scalar> *fwdSolveCriteria,
67  const EThrowOnSolveFailure throwOnFwdSolveFailure,
68  const SolveCriteria<Scalar> *adjSolveCriteria,
69  const EThrowOnSolveFailure throwOnAdjSolveFailure
70  )
71 {
72  initializeImpl(
73  lows,fwdSolveCriteria,throwOnFwdSolveFailure
74  ,adjSolveCriteria,throwOnAdjSolveFailure
75  );
76 }
77 
78 
79 template<class Scalar>
81  const Teuchos::RCP<const LinearOpWithSolveBase<Scalar> > &lows,
82  const SolveCriteria<Scalar> *fwdSolveCriteria,
83  const EThrowOnSolveFailure throwOnFwdSolveFailure,
84  const SolveCriteria<Scalar> *adjSolveCriteria,
85  const EThrowOnSolveFailure throwOnAdjSolveFailure
86  )
87 {
88  initializeImpl(
89  lows,fwdSolveCriteria,throwOnFwdSolveFailure
90  ,adjSolveCriteria,throwOnAdjSolveFailure
91  );
92 }
93 
94 
95 template<class Scalar>
98  const SolveCriteria<Scalar> *fwdSolveCriteria,
99  const EThrowOnSolveFailure throwOnFwdSolveFailure,
100  const SolveCriteria<Scalar> *adjSolveCriteria,
101  const EThrowOnSolveFailure throwOnAdjSolveFailure
102  )
103 {
104  initializeImpl(
105  lows,fwdSolveCriteria,throwOnFwdSolveFailure
106  ,adjSolveCriteria,throwOnAdjSolveFailure
107  );
108 }
109 
110 
111 template<class Scalar>
113  const Teuchos::RCP<const LinearOpWithSolveBase<Scalar> > &lows,
114  const SolveCriteria<Scalar> *fwdSolveCriteria,
115  const EThrowOnSolveFailure throwOnFwdSolveFailure,
116  const SolveCriteria<Scalar> *adjSolveCriteria,
117  const EThrowOnSolveFailure throwOnAdjSolveFailure
118  )
119 {
120  initializeImpl(
121  lows,fwdSolveCriteria,throwOnFwdSolveFailure
122  ,adjSolveCriteria,throwOnAdjSolveFailure
123  );
124 }
125 
126 
127 template<class Scalar>
129 {
130  lows_.uninitialize();
131  fwdSolveCriteria_ = Teuchos::null;
132  adjSolveCriteria_ = Teuchos::null;
133 }
134 
135 
136 // Overridden form InverseLinearOpBase
137 
138 
139 template<class Scalar>
141 {
142  return lows_.isConst();
143 }
144 
145 
146 template<class Scalar>
149 {
150  return lows_.getNonconstObj();
151 }
152 
153 
154 template<class Scalar>
157 {
158  return lows_.getConstObj();
159 }
160 
161 
162 // Overridden from LinearOpBase
163 
164 
165 template<class Scalar>
168 {
169  assertInitialized();
170  return lows_.getConstObj()->domain();
171 }
172 
173 
174 template<class Scalar>
177 {
178  assertInitialized();
179  return lows_.getConstObj()->range();
180 }
181 
182 
183 template<class Scalar>
186 {
187  return Teuchos::null; // Not supported yet but could be!
188 }
189 
190 
191 // Overridden from Teuchos::Describable
192 
193 
194 template<class Scalar>
196 {
197  assertInitialized();
198  std::ostringstream oss;
199  oss
201  << "lows="<<lows_.getConstObj()->description()
202  << ",fwdSolveCriteria="<<(fwdSolveCriteria_.get()?"...":"DEFAULT")
203  << ",adjSolveCriteria="<<(adjSolveCriteria_.get()?"...":"DEFAULT")
204  << "}";
205  return oss.str();
206 }
207 
208 
209 template<class Scalar>
212  const Teuchos::EVerbosityLevel verbLevel
213  ) const
214 {
215  using Teuchos::RCP;
216  using Teuchos::OSTab;
217  assertInitialized();
218  OSTab tab(out);
219  switch(verbLevel) {
221  case Teuchos::VERB_LOW:
222  out << this->description() << std::endl;
223  break;
225  case Teuchos::VERB_HIGH:
227  {
228  out
230  << "rangeDim=" << this->range()->dim()
231  << ",domainDim=" << this->domain()->dim() << "}:\n";
232  OSTab tab2(out);
233  out << "lows = ";
234  if(!lows_.getConstObj().get()) {
235  out << " NULL\n";
236  }
237  else {
238  out << Teuchos::describe(*lows_.getConstObj(),verbLevel);
239  }
240  break;
241  }
242  default:
243  TEUCHOS_TEST_FOR_EXCEPT(true); // Should never be called!
244  }
245 }
246 
247 
248 // protected
249 
250 
251 // Overridden from LinearOpBase
252 
253 
254 template<class Scalar>
256 {
257  if (nonnull(lows_)) {
258  return solveSupports(*lows_.getConstObj(),M_trans);
259  }
260  return false;
261 }
262 
263 
264 template<class Scalar>
266  const EOpTransp M_trans,
267  const MultiVectorBase<Scalar> &X,
268  const Ptr<MultiVectorBase<Scalar> > &Y,
269  const Scalar alpha,
270  const Scalar beta
271  ) const
272 {
274  assertInitialized();
275  // ToDo: Put in hooks for propogating verbosity level
276  //
277  // Y = alpha*op(M)*X + beta*Y
278  //
279  // =>
280  //
281  // Y = beta*Y
282  // Y += alpha*inv(op(lows))*X
283  //
285  if(beta==ST::zero()) {
286  T = Teuchos::rcpFromPtr(Y);
287  }
288  else {
289  T = createMembers(Y->range(),Y->domain()->dim());
290  scale(beta, Y);
291  }
292  //
293  const Ptr<const SolveCriteria<Scalar> > solveCriteria =
294  (
295  real_trans(M_trans)==NOTRANS
296  ? fwdSolveCriteria_.ptr()
297  : adjSolveCriteria_.ptr()
298  );
299  assign(T.ptr(), ST::zero()); // Have to initialize before solve!
300  SolveStatus<Scalar> solveStatus =
301  Thyra::solve<Scalar>(*lows_.getConstObj(), M_trans, X, T.ptr(), solveCriteria);
302 
304  nonnull(solveCriteria) && solveStatus.solveStatus!=SOLVE_STATUS_CONVERGED
305  && ( real_trans(M_trans)==NOTRANS
306  ? throwOnFwdSolveFailure_==THROW_ON_SOLVE_FAILURE
307  : throwOnAdjSolveFailure_==THROW_ON_SOLVE_FAILURE )
309  ,"Error, the LOWS object " << lows_.getConstObj()->description() << " returned an unconverged"
310  "status of " << toString(solveStatus.solveStatus) << " with the message "
311  << solveStatus.message << "."
312  );
313  //
314  if(beta==ST::zero()) {
315  scale(alpha, Y);
316  }
317  else {
318  update( alpha, *T, Y );
319  }
320 }
321 
322 
323 // private
324 
325 
326 template<class Scalar>
327 template<class LOWS>
329  const Teuchos::RCP<LOWS> &lows,
330  const SolveCriteria<Scalar> *fwdSolveCriteria,
331  const EThrowOnSolveFailure throwOnFwdSolveFailure,
332  const SolveCriteria<Scalar> *adjSolveCriteria,
333  const EThrowOnSolveFailure throwOnAdjSolveFailure
334  )
335 {
336  lows_.initialize(lows);
337  if(fwdSolveCriteria)
338  fwdSolveCriteria_ = Teuchos::rcp(new SolveCriteria<Scalar>(*fwdSolveCriteria));
339  else
340  fwdSolveCriteria_ = Teuchos::null;
341  if(adjSolveCriteria)
342  adjSolveCriteria_ = Teuchos::rcp(new SolveCriteria<Scalar>(*adjSolveCriteria));
343  else
344  adjSolveCriteria_ = Teuchos::null;
345  throwOnFwdSolveFailure_ = throwOnFwdSolveFailure;
346  throwOnAdjSolveFailure_ = throwOnAdjSolveFailure;
347  const std::string lowsLabel = lows_.getConstObj()->getObjectLabel();
348  if(lowsLabel.length())
349  this->setObjectLabel( "inv("+lowsLabel+")" );
350 }
351 
352 
353 } // end namespace Thyra
354 
355 
356 // Related non-member functions
357 
358 
359 template<class Scalar>
361 Thyra::nonconstInverse(
362  const RCP<LinearOpWithSolveBase<Scalar> > &A,
363  const Ptr<const SolveCriteria<Scalar> > &fwdSolveCriteria,
364  const EThrowOnSolveFailure throwOnFwdSolveFailure,
365  const Ptr<const SolveCriteria<Scalar> > &adjSolveCriteria,
366  const EThrowOnSolveFailure throwOnAdjSolveFailure
367  )
368 {
369  return Teuchos::rcp(
370  new DefaultInverseLinearOp<Scalar>(
371  A, fwdSolveCriteria.get(), throwOnFwdSolveFailure,
372  adjSolveCriteria.get(), throwOnAdjSolveFailure
373  )
374  );
375 }
376 
377 template<class Scalar>
379 Thyra::inverse(
380  const RCP<const LinearOpWithSolveBase<Scalar> > &A,
381  const Ptr<const SolveCriteria<Scalar> > &fwdSolveCriteria,
382  const EThrowOnSolveFailure throwOnFwdSolveFailure,
383  const Ptr<const SolveCriteria<Scalar> > &adjSolveCriteria,
384  const EThrowOnSolveFailure throwOnAdjSolveFailure
385  )
386 {
387  return Teuchos::rcp(
388  new DefaultInverseLinearOp<Scalar>(
389  A, fwdSolveCriteria.get(), throwOnFwdSolveFailure,
390  adjSolveCriteria.get(), throwOnAdjSolveFailure
391  )
392  );
393 }
394 
395 
396 //
397 // Explicit instantiation macro
398 //
399 // Must be expanded from within the Thyra namespace!
400 //
401 
402 
403 #define THYRA_DEFAULT_INVERSE_LINEAR_OP_INSTANT(SCALAR) \
404  \
405  template class DefaultInverseLinearOp<SCALAR >; \
406  \
407  template RCP<LinearOpBase<SCALAR > > \
408  nonconstInverse( \
409  const RCP<LinearOpWithSolveBase<SCALAR > > &A, \
410  const Ptr<const SolveCriteria<SCALAR > > &fwdSolveCriteria, \
411  const EThrowOnSolveFailure throwOnFwdSolveFailure, \
412  const Ptr<const SolveCriteria<SCALAR > > &adjSolveCriteria, \
413  const EThrowOnSolveFailure throwOnAdjSolveFailure \
414  ); \
415  \
416  template RCP<LinearOpBase<SCALAR > > \
417  inverse( \
418  const RCP<const LinearOpWithSolveBase<SCALAR > > &A, \
419  const Ptr<const SolveCriteria<SCALAR > > &fwdSolveCriteria, \
420  const EThrowOnSolveFailure throwOnFwdSolveFailure, \
421  const Ptr<const SolveCriteria<SCALAR > > &adjSolveCriteria, \
422  const EThrowOnSolveFailure throwOnAdjSolveFailure \
423  );
424 
425 
426 #endif // THYRA_DEFAULT_INVERSE_LINEAR_OP_DEF_HPP
Base class for all linear operators that can support a high-level solve operation.
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
basic_OSTab< char > OSTab
RCP< const LinearOpBase< Scalar > > clone() const
DefaultInverseLinearOp()
Constructs to uninitialized (see postconditions for uninitialize()).
bool opSupportedImpl(EOpTransp M_trans) const
Returns true only if all constituent operators support M_trans.
void describe(FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Use the non-transposed operator.
EOpTransp real_trans(EOpTransp transp)
Return NOTRANS or TRANS for real scalar valued operators and this also is used for determining struct...
RCP< const VectorSpaceBase< Scalar > > range() const
Returns this-&gt;getLows()-&gt;domain() if &lt;t&gt;this-&gt;getLows().get()!=NULL and returns Teuchos::null otherwi...
void applyImpl(const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha, const Scalar beta) const
Throw an exception if a solve fails to converge.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Interface for a collection of column vectors called a multi-vector.
Ptr< T > ptr() const
virtual std::string description() const
TEUCHOSCORE_LIB_DLL_EXPORT std::string toString(const EVerbosityLevel verbLevel)
Simple struct for the return status from a solve.
void initialize(const RCP< LinearOpWithSolveBase< Scalar > > &lows, const SolveCriteria< Scalar > *fwdSolveCriteria=NULL, const EThrowOnSolveFailure throwOnFwdSolveFailure=THROW_ON_SOLVE_FAILURE, const SolveCriteria< Scalar > *adjSolveCriteria=NULL, const EThrowOnSolveFailure throwOnAdjSolveFailure=THROW_ON_SOLVE_FAILURE)
Initialize given a non-const LinearOpWithSolveBase object and an optional .
RCP< const VectorSpaceBase< Scalar > > domain() const
Returns this-&gt;getLows()-&gt;range() if &lt;t&gt;this-&gt;getLows().get()!=NULL and returns Teuchos::null otherwis...
bool nonnull(const boost::shared_ptr< T > &p)
EThrowOnSolveFailure
Determines what to do if inverse solve fails.
The requested solution criteria has likely been achieved.
Exception type thrown on an catastrophic solve failure.
Simple struct that defines the requested solution criteria for a solve.
RCP< const LinearOpWithSolveBase< Scalar > > getLows() const
Concrete LinearOpBase subclass that creates an implicit LinearOpBase object using the inverse action ...
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
RCP< LinearOpWithSolveBase< Scalar > > getNonconstLows()