Thyra Version of the Day
Loading...
Searching...
No Matches
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
52namespace Thyra {
53
54
55// Constructors/initializers/accessors
56
57
58template<class Scalar>
60{}
61
62
63template<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
79template<class Scalar>
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
95template<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
111template<class Scalar>
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
127template<class Scalar>
129{
130 lows_.uninitialize();
131 fwdSolveCriteria_ = Teuchos::null;
132 adjSolveCriteria_ = Teuchos::null;
133}
134
135
136// Overridden form InverseLinearOpBase
137
138
139template<class Scalar>
141{
142 return lows_.isConst();
143}
144
145
146template<class Scalar>
149{
150 return lows_.getNonconstObj();
151}
152
153
154template<class Scalar>
157{
158 return lows_.getConstObj();
159}
160
161
162// Overridden from LinearOpBase
163
164
165template<class Scalar>
168{
169 assertInitialized();
170 return lows_.getConstObj()->domain();
171}
172
173
174template<class Scalar>
177{
178 assertInitialized();
179 return lows_.getConstObj()->range();
180}
181
182
183template<class Scalar>
186{
187 return Teuchos::null; // Not supported yet but could be!
188}
189
190
191// Overridden from Teuchos::Describable
192
193
194template<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
209template<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) {
222 out << this->description() << std::endl;
223 break;
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
254template<class Scalar>
256{
257 if (nonnull(lows_)) {
258 return solveSupports(*lows_.getConstObj(),M_trans);
259 }
260 return false;
261}
262
263
264template<class Scalar>
266 const EOpTransp M_trans,
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
326template<class Scalar>
327template<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
359template<class Scalar>
361Thyra::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
377template<class Scalar>
379Thyra::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
virtual std::string description() const
Ptr< T > ptr() const
Exception type thrown on an catastrophic solve failure.
Concrete LinearOpBase subclass that creates an implicit LinearOpBase object using the inverse action ...
void applyImpl(const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha, const Scalar beta) const
RCP< const LinearOpWithSolveBase< Scalar > > getLows() const
RCP< const LinearOpBase< Scalar > > clone() const
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 .
DefaultInverseLinearOp()
Constructs to uninitialized (see postconditions for uninitialize()).
RCP< const VectorSpaceBase< Scalar > > domain() const
Returns this->getLows()->range() if <t>this->getLows().get()!=NULL and returns Teuchos::null otherwis...
RCP< const VectorSpaceBase< Scalar > > range() const
Returns this->getLows()->domain() if <t>this->getLows().get()!=NULL and returns Teuchos::null otherwi...
void describe(FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
RCP< LinearOpWithSolveBase< Scalar > > getNonconstLows()
bool opSupportedImpl(EOpTransp M_trans) const
Returns true only if all constituent operators support M_trans.
Base class for all linear operators that can support a high-level solve operation.
Interface for a collection of column vectors called a multi-vector.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
@ SOLVE_STATUS_CONVERGED
The requested solution criteria has likely been achieved.
EThrowOnSolveFailure
Determines what to do if inverse solve fails.
@ THROW_ON_SOLVE_FAILURE
Throw an exception if a solve fails to converge.
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
const char * toString(EConj conj)
Return a string name for a EOpTransp value. `*.
EOpTransp real_trans(EOpTransp transp)
Return NOTRANS or TRANS for real scalar valued operators and this also is used for determining struct...
@ NOTRANS
Use the non-transposed operator.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Simple struct that defines the requested solution criteria for a solve.
Simple struct for the return status from a solve.
std::string message
A simple one-line message (i.e. no newlines) returned from the solver.
ESolveStatus solveStatus
The return status of the solve.