Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_DefaultAddedLinearOp_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_ADDED_LINEAR_OP_DEF_HPP
43#define THYRA_DEFAULT_ADDED_LINEAR_OP_DEF_HPP
44
45#include "Thyra_DefaultAddedLinearOp_decl.hpp"
46#include "Thyra_DefaultScaledAdjointLinearOp.hpp"
47#include "Thyra_AssertOp.hpp"
48#include "Teuchos_Utils.hpp"
49
50
51namespace Thyra {
52
53
54// Inline members only used in this class impl
55
56
57template<class Scalar>
58inline
59void DefaultAddedLinearOp<Scalar>::assertInitialized() const
60{
61#ifdef TEUCHOS_DEBUG
62 TEUCHOS_TEST_FOR_EXCEPT( !( numOps() > 0 ) );
63#endif
64}
65
66
67template<class Scalar>
68inline
69std::string DefaultAddedLinearOp<Scalar>::getClassName() const
70{
72}
73
74
75template<class Scalar>
76inline
77Ordinal DefaultAddedLinearOp<Scalar>::getRangeDim() const
78{
79 return (numOps() > 0 ? this->range()->dim() : 0);
80}
81
82
83template<class Scalar>
84inline
85Ordinal DefaultAddedLinearOp<Scalar>::getDomainDim() const
86{
87 return (numOps() > 0 ? this->domain()->dim() : 0);
88}
89
90
91// Constructors/initializers/accessors
92
93
94template<class Scalar>
97
98
99template<class Scalar>
101 const ArrayView<const RCP<LinearOpBase<Scalar> > > &Ops )
102{
103 initialize(Ops);
104}
105
106
107template<class Scalar>
109 const ArrayView<const RCP<const LinearOpBase<Scalar> > > &Ops )
110{
111 initialize(Ops);
112}
113
114
115template<class Scalar>
117 const ArrayView<const RCP<LinearOpBase<Scalar> > > &Ops )
118{
119 const int l_numOps = Ops.size();
120 Ops_.resize(l_numOps);
121 for( int k = 0; k < l_numOps; ++k )
122 Ops_[k].initialize(Ops[k]);
123 validateOps();
124 setupDefaultObjectLabel();
125}
126
127
128template<class Scalar>
130 const ArrayView<const RCP<const LinearOpBase<Scalar> > > &Ops )
131{
132 const int l_numOps = Ops.size();
133 Ops_.resize(l_numOps);
134 for( int k = 0; k < l_numOps; ++k )
135 Ops_[k].initialize(Ops[k]);
136 validateOps();
137 setupDefaultObjectLabel();
138}
139
140
141template<class Scalar>
143{
144 Ops_.resize(0);
145 setupDefaultObjectLabel();
146}
147
148
149// Overridden form AddedLinearOpBase
150
151
152template<class Scalar>
154{
155 return Ops_.size();
156}
157
158
159template<class Scalar>
161{
162#ifdef TEUCHOS_DEBUG
163 TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= k && k < numOps() ) );
164#endif
165 return Ops_[k].isConst();
166}
167
168
169template<class Scalar>
172{
173#ifdef TEUCHOS_DEBUG
174 TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= k && k < numOps() ) );
175#endif
176 return Ops_[k].getNonconstObj();
177}
178
179
180template<class Scalar>
183{
184#ifdef TEUCHOS_DEBUG
185 TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= k && k < numOps() ) );
186#endif
187 return Ops_[k];
188}
189
190
191// Overridden from LinearOpBase
192
193
194template<class Scalar>
197{
198 if (numOps()) {
199 return getOp(0)->range();
200 }
201 return Teuchos::null;
202}
203
204
205template<class Scalar>
208{
209 if (numOps()) {
210 return getOp(numOps()-1)->domain();
211 }
212 return Teuchos::null;
213}
214
215
216template<class Scalar>
219{
220 return Teuchos::null; // Not supported yet but could be!
221}
222
223
224// Overridden from Teuchos::Describable
225
226
227template<class Scalar>
229{
230 std::ostringstream oss;
231 oss << getClassName() << "{numOps="<<numOps()
232 <<",rangeDim=" << getRangeDim()
233 << ",domainDim="<< getDomainDim() <<"}";
234 return oss.str();
235}
236
237
238template<class Scalar>
240 Teuchos::FancyOStream &out_arg
241 ,const Teuchos::EVerbosityLevel verbLevel
242 ) const
243{
244 using Teuchos::RCP;
246 using Teuchos::OSTab;
247 RCP<FancyOStream> out = rcp(&out_arg,false);
248 OSTab tab(out);
249 const int l_numOps = Ops_.size();
250 switch(verbLevel) {
253 *out << this->description() << std::endl;
254 break;
258 {
259 *out << this->description() << std::endl;
260 OSTab tab2(out);
261 *out
262 << "Constituent LinearOpBase objects for M = Op[0]*...*Op[numOps-1]:\n";
263 tab.incrTab();
264 for( int k = 0; k < l_numOps; ++k ) {
265 *out << "Op["<<k<<"] = " << Teuchos::describe(*getOp(k),verbLevel);
266 }
267 break;
268 }
269 default:
270 TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
271 }
272}
273
274
275// protected
276
277
278// Overridden from LinearOpBase
279
280
281template<class Scalar>
283{
284 bool isOpSupported = true;
285 for( int k = 0; k < static_cast<int>(Ops_.size()); ++k )
286 if(!Thyra::opSupported(*getOp(k),M_trans)) isOpSupported = false;
287 return isOpSupported;
288 // ToDo: Cache these?
289}
290
291
292template<class Scalar>
294 const EOpTransp M_trans,
296 const Ptr<MultiVectorBase<Scalar> > &Y,
297 const Scalar alpha,
298 const Scalar beta
299 ) const
300{
302#ifdef TEUCHOS_DEBUG
304 getClassName()+"::apply(...)", *this, M_trans, X, &*Y
305 );
306#endif // TEUCHOS_DEBUG
307 //
308 // Y = alpha * op(M) * X + beta*Y
309 //
310 // =>
311 //
312 // Y = beta*Y + sum(alpha*op(Op[j])*X),j=0...numOps-1)
313 //
314 const int l_numOps = Ops_.size();
315 for( int j = 0; j < l_numOps; ++j )
316 Thyra::apply(*getOp(j), M_trans, X, Y, alpha, j==0 ? beta : ST::one());
317}
318
319
320// private
321
322
323template<class Scalar>
325{
326 using Teuchos::toString;
327#ifdef TEUCHOS_DEBUG
328 try {
329 const int l_numOps = Ops_.size();
330 for( int k = 0; k < l_numOps; ++k ) {
331 TEUCHOS_TEST_FOR_EXCEPT( Ops_[k]().get() == NULL );
332 if( k > 0 ) {
334 getClassName()+"::initialize(...)"
335 ,*Ops_[0], NOTRANS, ("Ops[0]")
336 ,*Ops_[k], NOTRANS, ("Ops["+toString(k)+"]")
337 );
338 }
339 }
340 }
341 catch(...) {
342 uninitialize();
343 throw;
344 }
345#endif
346}
347
348
349template<class Scalar>
350void DefaultAddedLinearOp<Scalar>::setupDefaultObjectLabel()
351{
352 std::ostringstream label;
353 const int l_numOps = Ops_.size();
354 for( int k = 0; k < l_numOps; ++k ) {
355 std::string Op_k_label = Ops_[k]->getObjectLabel();
356 if (Op_k_label.length() == 0)
357 Op_k_label = "ANYM";
358 if (k > 0)
359 label << "+";
360 label << "("<<Op_k_label<<")";
361 }
362 this->setObjectLabel(label.str());
363 validateOps();
364}
365
366
367} // end namespace Thyra
368
369
370template<class Scalar>
372Thyra::nonconstAdd(
373 const RCP<LinearOpBase<Scalar> > &A,
374 const RCP<LinearOpBase<Scalar> > &B,
375 const std::string &label
376 )
377{
378 using Teuchos::tuple;
379 RCP<LinearOpBase<Scalar> > alo =
380 defaultAddedLinearOp<Scalar>(
381 tuple<RCP<LinearOpBase<Scalar> > >(A, B)()
382 );
383 if (label.length())
384 alo->setObjectLabel(label);
385 return alo;
386}
387
388
389template<class Scalar>
391Thyra::add(
392 const RCP<const LinearOpBase<Scalar> > &A,
393 const RCP<const LinearOpBase<Scalar> > &B,
394 const std::string &label
395 )
396{
397 using Teuchos::tuple;
398 RCP<LinearOpBase<Scalar> > alo =
399 defaultAddedLinearOp<Scalar>(
400 tuple<RCP<const LinearOpBase<Scalar> > >(A, B)()
401 );
402 if (label.length())
403 alo->setObjectLabel(label);
404 return alo;
405}
406
407
408template<class Scalar>
410Thyra::nonconstSubtract(
411 const RCP<LinearOpBase<Scalar> > &A,
412 const RCP<LinearOpBase<Scalar> > &B,
413 const std::string &label
414 )
415{
416 typedef ScalarTraits<Scalar> ST;
417 using Teuchos::tuple;
418 RCP<LinearOpBase<Scalar> > alo =
419 defaultAddedLinearOp<Scalar>(
420 tuple<RCP<LinearOpBase<Scalar> > >(
421 A, nonconstScale<Scalar>(-ST::one(),B) )()
422 );
423 if (label.length())
424 alo->setObjectLabel(label);
425 return alo;
426}
427
428
429template<class Scalar>
431Thyra::subtract(
432 const RCP<const LinearOpBase<Scalar> > &A,
433 const RCP<const LinearOpBase<Scalar> > &B,
434 const std::string &label
435 )
436{
437 typedef ScalarTraits<Scalar> ST;
438 using Teuchos::tuple;
439 RCP<LinearOpBase<Scalar> > alo =
440 defaultAddedLinearOp<Scalar>(
441 tuple<RCP<const LinearOpBase<Scalar> > >(
442 A, scale<Scalar>(-ST::one(),B)
443 )()
444 );
445 if (label.length())
446 alo->setObjectLabel(label);
447 return alo;
448}
449
450
451//
452// Explicit instantiation macro
453//
454
455#define THYRA_DEFAULT_ADDED_LINEAR_OP_INSTANT(SCALAR) \
456 template class DefaultAddedLinearOp<SCALAR >; \
457 \
458 template RCP<LinearOpBase<SCALAR > > \
459 nonconstAdd( \
460 const RCP<LinearOpBase<SCALAR > > &A, \
461 const RCP<LinearOpBase<SCALAR > > &B, \
462 const std::string &label \
463 ); \
464 \
465 template RCP<const LinearOpBase<SCALAR > > \
466 add( \
467 const RCP<const LinearOpBase<SCALAR > > &A, \
468 const RCP<const LinearOpBase<SCALAR > > &B, \
469 const std::string &label \
470 ); \
471 \
472 template RCP<LinearOpBase<SCALAR > > \
473 nonconstSubtract( \
474 const RCP<LinearOpBase<SCALAR > > &A, \
475 const RCP<LinearOpBase<SCALAR > > &B, \
476 const std::string &label \
477 ); \
478 \
479 template RCP<const LinearOpBase<SCALAR > > \
480 subtract( \
481 const RCP<const LinearOpBase<SCALAR > > &A, \
482 const RCP<const LinearOpBase<SCALAR > > &B, \
483 const std::string &label \
484 );
485
486
487#endif // THYRA_DEFAULT_ADDED_LINEAR_OP_DEF_HPP
virtual std::string description() const
Concrete composite LinearOpBase subclass that creates an implicitly added linear operator out of one ...
RCP< LinearOpBase< Scalar > > getNonconstOp(const int k)
void applyImpl(const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha, const Scalar beta) const
std::string description() const
Prints just the name DefaultAddedLinearOp along with the overall dimensions and the number of constit...
RCP< const VectorSpaceBase< Scalar > > domain() const
Returns this->getOp(this->numOps()-1).domain() if <t>this->numOps() > 0 and returns Teuchos::null oth...
RCP< const VectorSpaceBase< Scalar > > range() const
Returns this->getOp(0).range() if <t>this->numOps() > 0 and returns Teuchos::null otherwise.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Prints the details about the constituent linear operators.
DefaultAddedLinearOp()
Constructs to uninitialized.
void initialize(const ArrayView< const RCP< LinearOpBase< Scalar > > > &Ops)
Initialize given a list of non-const linear operators.
bool opSupportedImpl(EOpTransp M_trans) const
Returns true only if all constituent operators support M_trans.
RCP< const LinearOpBase< Scalar > > clone() const
RCP< const LinearOpBase< Scalar > > getOp(const int k) const
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 THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(FUNC_NAME, M, M_T, X, Y)
This is a very useful macro that should be used to validate that the spaces for the multi-vector vers...
#define THYRA_ASSERT_LINEAR_OP_PLUS_LINEAR_OP_SPACES_NAMES(FUNC_NAME, M1, M1_T, M1_N, M2, M2_T, M2_N)
Assert that a linear operator addition matches up.
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
const char * toString(EConj conj)
Return a string name for a EOpTransp value. `*.
@ NOTRANS
Use the non-transposed operator.
T_To & dyn_cast(T_From &from)
std::string toString(const T &t)