Thyra  Version of the Day
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>
65  const Teuchos::RCP<LinearOpWithSolveBase<Scalar> > &lows,
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>
97  const Teuchos::RCP<LinearOpWithSolveBase<Scalar> > &lows,
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>
147 Teuchos::RCP<LinearOpWithSolveBase<Scalar> >
149 {
150  return lows_.getNonconstObj();
151 }
152 
153 
154 template<class Scalar>
155 Teuchos::RCP<const LinearOpWithSolveBase<Scalar> >
157 {
158  return lows_.getConstObj();
159 }
160 
161 
162 // Overridden from LinearOpBase
163 
164 
165 template<class Scalar>
166 Teuchos::RCP< const VectorSpaceBase<Scalar> >
168 {
169  assertInitialized();
170  return lows_.getConstObj()->domain();
171 }
172 
173 
174 template<class Scalar>
175 Teuchos::RCP< const VectorSpaceBase<Scalar> >
177 {
178  assertInitialized();
179  return lows_.getConstObj()->range();
180 }
181 
182 
183 template<class Scalar>
184 Teuchos::RCP<const LinearOpBase<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
200  << Teuchos::Describable::description() << "{"
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>
211  Teuchos::FancyOStream &out,
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) {
220  case Teuchos::VERB_DEFAULT:
221  case Teuchos::VERB_LOW:
222  out << this->description() << std::endl;
223  break;
224  case Teuchos::VERB_MEDIUM:
225  case Teuchos::VERB_HIGH:
226  case Teuchos::VERB_EXTREME:
227  {
228  out
229  << Teuchos::Describable::description() << "{"
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 {
273  typedef Teuchos::ScalarTraits<Scalar> ST;
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  //
284  Teuchos::RCP<MultiVectorBase<Scalar> > T;
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 
303  TEUCHOS_TEST_FOR_EXCEPTION(
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>
360 Teuchos::RCP<Thyra::LinearOpBase<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>
378 Teuchos::RCP<Thyra::LinearOpBase<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.
RCP< const VectorSpaceBase< Scalar > > domain() const
Returns this->getLows()->range() if <t>this->getLows().get()!=NULL and returns Teuchos::null otherwis...
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
DefaultInverseLinearOp()
Constructs to uninitialized (see postconditions for uninitialize()).
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...
bool opSupportedImpl(EOpTransp M_trans) const
Returns true only if all constituent operators support M_trans.
RCP< const LinearOpBase< Scalar > > clone() const
RCP< const LinearOpWithSolveBase< Scalar > > getLows() const
RCP< const VectorSpaceBase< Scalar > > range() const
Returns this->getLows()->domain() if <t>this->getLows().get()!=NULL and returns Teuchos::null otherwi...
Throw an exception if a solve fails to converge.
Interface for a collection of column vectors called a multi-vector.
void describe(FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
void applyImpl(const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha, const Scalar beta) const
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 .
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.
const char * toString(EConj conj)
Return a string name for a EOpTransp value. `*.
Simple struct that defines the requested solution criteria for a solve.
Concrete LinearOpBase subclass that creates an implicit LinearOpBase object using the inverse action ...
RCP< LinearOpWithSolveBase< Scalar > > getNonconstLows()