Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_DefaultScaledAdjointLinearOp_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
43#ifndef THYRA_DEFAULT_SCALED_ADJOINT_LINEAR_OP_DEF_HPP
44#define THYRA_DEFAULT_SCALED_ADJOINT_LINEAR_OP_DEF_HPP
45
46
47#include "Thyra_DefaultScaledAdjointLinearOp_decl.hpp"
48#include "Thyra_VectorSpaceBase.hpp"
49#include "Thyra_AssertOp.hpp"
50
51
52namespace Thyra {
53
54
55//Constructors/initializers/accessors
56
57
58template<class Scalar>
60 const Scalar &scalar
61 ,const EOpTransp &transp
63 )
64{
65 initializeImpl(scalar,transp,Op,false);
66}
67
68
69template<class Scalar>
71 const Scalar &scalar
72 ,const EOpTransp &transp
73 ,const Teuchos::RCP<const LinearOpBase<Scalar> > &Op
74 )
75{
76 initializeImpl(scalar,transp,Op,false);
77}
78
79
80template<class Scalar>
83{
84 return getOpImpl().getNonconstObj();
85}
86
87
88template<class Scalar>
91{
92 return getOpImpl();
93}
94
95
96template<class Scalar>
98{
100 origOp_.uninitialize();
101 overallScalar_ = ST::zero();
102 overallTransp_ = NOTRANS;
103 allScalarETransp_ = Teuchos::null;
104
105}
106
107
108// Overridden from Teuchos::Describable
109
110
111template<class Scalar>
113{
114 assertInitialized();
115 std::ostringstream oss;
117 << overallScalar() << ","<<toString(overallTransp())<<","
118 << origOp_->description() << "}";
119 return oss.str();
120}
121
122
123template<class Scalar>
125 Teuchos::FancyOStream &out_arg
126 ,const Teuchos::EVerbosityLevel verbLevel
127 ) const
128{
130 using Teuchos::RCP;
132 using Teuchos::OSTab;
133 assertInitialized();
134 RCP<FancyOStream> out = rcp(&out_arg,false);
135 OSTab tab(out);
136 switch(verbLevel) {
139 *out << this->description() << std::endl;
140 break;
144 {
145 *out
147 << "rangeDim=" << this->range()->dim()
148 << ",domainDim=" << this->domain()->dim() << "}\n";
149 OSTab tab2(out);
150 *out
151 << "overallScalar="<< overallScalar() << std::endl
152 << "overallTransp="<<toString(overallTransp()) << std::endl
153 << "Constituent transformations:\n";
154 for( int i = 0; i <= my_index_; ++i ) {
155 const ScalarETransp<Scalar> &scalar_transp = (*allScalarETransp_)[my_index_-i];
156 OSTab tab3(out,i+1);
157 if(scalar_transp.scalar != ST::one() && scalar_transp.transp != NOTRANS)
158 *out << "scalar="<<scalar_transp.scalar<<",transp="<<toString(scalar_transp.transp)<<std::endl;
159 else if(scalar_transp.scalar != ST::one())
160 *out << "scalar="<<scalar_transp.scalar<<std::endl;
161 else if( scalar_transp.transp != NOTRANS )
162 *out << "transp="<<toString(scalar_transp.transp)<<std::endl;
163 else
164 *out << "no-transformation\n";
165 }
166 tab.incrTab(my_index_+2);
167 *out << "origOp = " << Teuchos::describe(*origOp_,verbLevel);
168 break;
169 }
170 default:
171 TEUCHOS_TEST_FOR_EXCEPT(true); // Should never be called!
172 }
173}
174
175
176// Overridden from LinearOpBase
177
178
179template<class Scalar>
182{
183 assertInitialized();
184 return ( real_trans(this->overallTransp())==NOTRANS
185 ? this->getOrigOp()->range() : this->getOrigOp()->domain() );
186}
187
188
189template<class Scalar>
192{
193 assertInitialized();
194 return ( real_trans(this->overallTransp())==NOTRANS
195 ? this->getOrigOp()->domain() : this->getOrigOp()->range() );
196}
197
198
199template<class Scalar>
202{
203 return Teuchos::null; // Not supported yet but could be
204}
205
206
207// Overridden from ScaledAdointLinearOpBase
208
209
210template<class Scalar>
212{
213 return overallScalar_;
214}
215
216
217template<class Scalar>
219{
220 return overallTransp_;
221}
222
223
224template<class Scalar>
227{
228 return origOp_.getNonconstObj();
229}
230
231
232template<class Scalar>
235{
236 return origOp_;
237}
238
239
240// protected
241
242
243// Overridden from LinearOpBase
244
245
246template<class Scalar>
248{
249 assertInitialized();
250 return Thyra::opSupported(
251 *this->getOrigOp(), trans_trans(this->overallTransp(), M_trans) );
252}
253
254
255template<class Scalar>
257 const EOpTransp M_trans,
259 const Ptr<MultiVectorBase<Scalar> > &Y,
260 const Scalar alpha,
261 const Scalar beta
262 ) const
263{
264 using Teuchos::as;
265 assertInitialized();
266 Thyra::apply(
267 *this->getOrigOp(), trans_trans(M_trans, this->overallTransp()),
268 X, Y, as<Scalar>(this->overallScalar()*alpha), beta
269 );
270}
271
272
273// private
274
275
276template<class Scalar>
278 const Scalar &scalar
279 ,const EOpTransp &transp
280 ,const Teuchos::RCP<const LinearOpBase<Scalar> > &Op
281 ,const bool isConst
282 )
283{
285#ifdef TEUCHOS_DEBUG
287 Op.get()==NULL, std::invalid_argument
288 ,"DefaultScaledAdjointLinearOp<"<<ST::name()<<">::initialize(scalar,transp,Op): "
289 "Error!, Op.get()==NULL is not allowed!"
290 );
291#endif // TEUCHOS_DEBUG
293 saOp = Teuchos::rcp_dynamic_cast<const DefaultScaledAdjointLinearOp<Scalar> >(Op);
294 if(saOp.get()) {
295 origOp_ = saOp->origOp_;
296 overallScalar_ = saOp->overallScalar_*scalar;
297 overallTransp_ = trans_trans(saOp->overallTransp_,transp) ;
298 my_index_ = saOp->my_index_ + 1;
299 allScalarETransp_ = saOp->allScalarETransp_;
300 }
301 else {
302 if(isConst)
303 origOp_.initialize(Op);
304 else
305 origOp_.initialize(Teuchos::rcp_const_cast<LinearOpBase<Scalar> >(Op));
306 overallScalar_ = scalar;
307 overallTransp_ = transp;
308 my_index_ = 0;
309 allScalarETransp_ = Teuchos::rcp(new allScalarETransp_t());
310 }
311 allScalarETransp_->push_back(ScalarETransp<Scalar>(scalar,transp));
312 // Set the object label
313 std::string Op_label = Op->getObjectLabel();
314 if(Op_label.length()==0)
315 Op_label = "ANYM";
316 std::ostringstream label;
317 if(scalar!=ST::one())
318 label << scalar << "*";
319 switch(transp) {
320 case NOTRANS:
321 break; // No output
322 case CONJ:
323 label << "conj";
324 break;
325 case TRANS:
326 label << "trans";
327 break;
328 case CONJTRANS:
329 label << "adj";
330 break;
331 default:
332 TEUCHOS_TEST_FOR_EXCEPT("Invalid EOpTransp value!");
333 }
334 label << "(" << Op_label << ")";
335 this->setObjectLabel(label.str());
336}
337
338
339template<class Scalar>
340typename DefaultScaledAdjointLinearOp<Scalar>::CNLOC
341DefaultScaledAdjointLinearOp<Scalar>::getOpImpl() const
342{
343 assertInitialized();
344 if( my_index_ > 0 ) {
345 const ScalarETransp<Scalar> &scalar_transp = allScalarETransp_->at(my_index_);
347 Op = Teuchos::rcp(new DefaultScaledAdjointLinearOp<Scalar>());
348 Op->origOp_ = origOp_;
349 Op->overallScalar_ = overallScalar_/scalar_transp.scalar;
350 Op->overallTransp_ = trans_trans(overallTransp_,scalar_transp.transp);
351 Op->my_index_ = my_index_ - 1;
352 Op->allScalarETransp_ = allScalarETransp_;
353 return CNLOC(
354 Teuchos::rcp_implicit_cast<LinearOpBase<Scalar> >(Op)
355 );
356 }
357 return origOp_;
358}
359
360
361} // namespace Thyra
362
363
364#endif // THYRA_DEFAULT_SCALED_ADJOINT_LINEAR_OP_DEF_HPP
virtual std::string description() const
T * get() const
Concrete decorator LinearOpBase subclass that wraps a LinearOpBase object and adds on an extra scalin...
RCP< LinearOpBase< Scalar > > getNonconstOp()
Return the non-const linear operator passed into initialize().
bool opSupportedImpl(EOpTransp M_trans) const
Return if the operation is supported on the logical linear operator.
std::string description() const
Outputs DefaultScaledAdjointLinearOp<Scalar>{this->getOrigOp().description()) along with the dimensio...
RCP< const VectorSpaceBase< Scalar > > domain() const
Return the domain space of the logical linear operator.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Prints out the original operator as well as all of the scalings and transpositions in the order that ...
void uninitialize()
Set to uninitialized and (optionally) extract the objects passed into initialize().
void applyImpl(const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha, const Scalar beta) const
Apply the linear operator (or its transpose) to a multi-vector : Y = alpha*op(M)*X + beta*Y.
RCP< const LinearOpBase< Scalar > > clone() const
RCP< const LinearOpBase< Scalar > > getOrigOp() const
RCP< const VectorSpaceBase< Scalar > > range() const
Return the range space of the logical linear operator.
void initialize(const Scalar &scalar, const EOpTransp &transp, const RCP< LinearOpBase< Scalar > > &Op)
Initialize with an operator with by defining adjoint (transpose) and scaling arguments.
RCP< const LinearOpBase< Scalar > > getOp() const
Return the const linear operator passed into initialize().
Base class for all linear operators.
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)
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...
EOpTransp trans_trans(EOpTransp trans1, EOpTransp trans2)
Combine two transpose arguments.
@ TRANS
Use the transposed operator.
@ NOTRANS
Use the non-transposed operator.
@ CONJTRANS
Use the transposed operator with complex-conjugate clements (same as TRANS for real scalar types).
@ CONJ
Use the non-transposed operator with complex-conjugate elements (same as NOTRANS for real scalar type...
TypeTo as(const TypeFrom &t)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)