Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_DefaultMultipliedLinearOp_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_MULTIPLIED_LINEAR_OP_DEF_HPP
43#define THYRA_DEFAULT_MULTIPLIED_LINEAR_OP_DEF_HPP
44
45#include "Thyra_DefaultMultipliedLinearOp_decl.hpp"
46#include "Thyra_AssertOp.hpp"
47
48
49namespace Thyra {
50
51
52// Inline members only used in implementation
53
54
55template<class Scalar>
56inline
57void DefaultMultipliedLinearOp<Scalar>::assertInitialized() const
58{
59#ifdef TEUCHOS_DEBUG
60 TEUCHOS_TEST_FOR_EXCEPT( !( numOps() > 0 ) );
61#endif
62}
63
64
65template<class Scalar>
66inline
67std::string DefaultMultipliedLinearOp<Scalar>::getClassName() const
68{
70}
71
72
73template<class Scalar>
74inline
75Ordinal DefaultMultipliedLinearOp<Scalar>::getRangeDim() const
76{
77 return (numOps() > 0 ? this->range()->dim() : 0);
78}
79
80
81template<class Scalar>
82inline
83Ordinal DefaultMultipliedLinearOp<Scalar>::getDomainDim() const
84{
85 return (numOps() > 0 ? this->domain()->dim() : 0);
86}
87
88
89// Constructors/initializers/accessors
90
91
92template<class Scalar>
94{}
95
96
97template<class Scalar>
99 const ArrayView<const RCP<LinearOpBase<Scalar> > > &Ops )
100{
101 const int nOps = Ops.size();
102 Ops_.resize(nOps);
103 for( int k = 0; k < nOps; ++k )
104 Ops_[k].initialize(Ops[k]);
105 validateOps();
106 setupDefaultObjectLabel();
107}
108
109
110template<class Scalar>
112 const ArrayView<const RCP<const LinearOpBase<Scalar> > > &Ops )
113{
114 const int nOps = Ops.size();
115 Ops_.resize(nOps);
116 for( int k = 0; k < nOps; ++k )
117 Ops_[k].initialize(Ops[k]);
118 validateOps();
119 setupDefaultObjectLabel();
120}
121
122
123template<class Scalar>
125{
126 Ops_.resize(0);
127 setupDefaultObjectLabel();
128}
129
130
131// Overridden form MultipliedLinearOpBase
132
133
134template<class Scalar>
136{
137 return Ops_.size();
138}
139
140
141template<class Scalar>
143{
144#ifdef TEUCHOS_DEBUG
145 TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= k && k < numOps() ) );
146#endif
147 return Ops_[k].isConst();
148}
149
150
151template<class Scalar>
154{
155#ifdef TEUCHOS_DEBUG
156 TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= k && k < numOps() ) );
157#endif
158 return Ops_[k].getNonconstObj();
159}
160
161
162template<class Scalar>
165{
166#ifdef TEUCHOS_DEBUG
167 TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= k && k < numOps() ) );
168#endif
169 return Ops_[k];
170}
171
172
173// Overridden from LinearOpBase
174
175
176template<class Scalar>
179{
180 if (numOps()) {
181 return getOp(0)->range();
182 }
183 return Teuchos::null;
184}
185
186
187template<class Scalar>
190{
191 if (numOps()) {
192 return getOp(numOps()-1)->domain();
193 }
194 return Teuchos::null;
195}
196
197
198template<class Scalar>
201{
202 return Teuchos::null; // Not supported yet but could be!
203}
204
205
206// Overridden from Teuchos::Describable
207
208
209template<class Scalar>
211{
212 std::ostringstream oss;
213 oss << getClassName() << "{numOps="<<numOps()
214 <<",rangeDim=" << getRangeDim()
215 << ",domainDim="<< getDomainDim() <<"}";
216 return oss.str();
217}
218
219template<class Scalar>
221 Teuchos::FancyOStream &out_arg,
222 const Teuchos::EVerbosityLevel verbLevel
223 ) const
224{
226 using Teuchos::OSTab;
227 RCP<FancyOStream> out = rcp(&out_arg,false);
228 OSTab tab(out);
229 const int nOps = Ops_.size();
230 switch(verbLevel) {
233 *out << this->description() << std::endl;
234 break;
238 {
239 *out << this->description() << std::endl;
240 OSTab tab2(out);
241 *out
242 << "Constituent LinearOpBase objects for M = Op[0]*...*Op[numOps-1]:\n";
243 OSTab tab3(out);
244 for( int k = 0; k < nOps; ++k ) {
245 *out << "Op["<<k<<"] = " << Teuchos::describe(*getOp(k),verbLevel);
246 }
247 break;
248 }
249 default:
250 TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
251 }
252}
253
254
255// protected
256
257
258// Overridden from LinearOpBase
259
260
261template<class Scalar>
263{
264 bool overallOpSupported = true;
265 for( int k = 0; k < static_cast<int>(Ops_.size()); ++k )
266 if(!Thyra::opSupported(*getOp(k),M_trans)) overallOpSupported = false;
267 return overallOpSupported;
268 // ToDo: Cache these?
269}
270
271
272template<class Scalar>
274 const EOpTransp M_trans,
276 const Ptr<MultiVectorBase<Scalar> > &Y,
277 const Scalar alpha,
278 const Scalar beta
279 ) const
280{
281 using Teuchos::rcpFromPtr;
282 using Teuchos::rcpFromRef;
283#ifdef TEUCHOS_DEBUG
285 getClassName()+"::apply(...)", *this, M_trans, X, &*Y
286 );
287#endif // TEUCHOS_DEBUG
288 const int nOps = Ops_.size();
289 const Ordinal m = X.domain()->dim();
290 if( real_trans(M_trans)==NOTRANS ) {
291 //
292 // Y = alpha * M * X + beta*Y
293 // =>
294 // Y = alpha * op(Op[0]) * op(Op[1]) * ... * op(Op[numOps-1]) * X + beta*Y
295 //
296 RCP<MultiVectorBase<Scalar> > T_kp1, T_k; // Temporary propagated between loops
297 for( int k = nOps-1; k >= 0; --k ) {
300 if(k==0) Y_k = rcpFromPtr(Y); else Y_k = T_k = createMembers(getOp(k)->range(), m);
301 if(k==nOps-1) X_k = rcpFromRef(X); else X_k = T_kp1;
302 if( k > 0 )
303 Thyra::apply(*getOp(k), M_trans, *X_k, Y_k.ptr());
304 else
305 Thyra::apply(*getOp(k), M_trans, *X_k, Y_k.ptr(), alpha, beta);
306 T_kp1 = T_k;
307 }
308 }
309 else {
310 //
311 // Y = alpha * M' * X + beta*Y
312 // =>
313 // Y = alpha * Op[numOps-1]' * Op[1]' * ... * Op[0]' * X + beta * Y
314 //
315 RCP<MultiVectorBase<Scalar> > T_km1, T_k; // Temporary propagated between loops
316 for( int k = 0; k <= nOps-1; ++k ) {
319 if(k==nOps-1) Y_k = rcpFromPtr(Y); else Y_k = T_k = createMembers(getOp(k)->domain(), m);
320 if(k==0) X_k = rcpFromRef(X); else X_k = T_km1;
321 if( k < nOps-1 )
322 Thyra::apply(*getOp(k), M_trans, *X_k, Y_k.ptr());
323 else
324 Thyra::apply(*getOp(k), M_trans, *X_k, Y_k.ptr(), alpha, beta);
325 T_km1 = T_k;
326 }
327 }
328}
329
330
331// private
332
333
334template<class Scalar>
336{
337 using Teuchos::toString;
338#ifdef TEUCHOS_DEBUG
339 try {
340 const int nOps = Ops_.size();
341 for( int k = 0; k < nOps; ++k ) {
342 TEUCHOS_TEST_FOR_EXCEPT( Ops_[k]().get() == NULL );
343 if( k < nOps-1 ) {
345 getClassName()+"::initialize(...)"
346 ,*Ops_[k],NOTRANS,("Ops["+toString(k)+"]")
347 ,*Ops_[k+1],NOTRANS,("Ops["+toString(k+1)+"]")
348 );
349 }
350 }
351 }
352 catch(...) {
353 uninitialize();
354 throw;
355 }
356#endif
357}
358
359
360template<class Scalar>
361void DefaultMultipliedLinearOp<Scalar>::setupDefaultObjectLabel()
362{
363 std::ostringstream label;
364 const int nOps = Ops_.size();
365 for( int k = 0; k < nOps; ++k ) {
366 std::string Op_k_label = Ops_[k]->getObjectLabel();
367 if (Op_k_label.length() == 0)
368 Op_k_label = "ANYM";
369 if (k > 0)
370 label << "*";
371 label << "("<<Op_k_label<<")";
372 }
373 this->setObjectLabel(label.str());
374 validateOps();
375}
376
377
378} // end namespace Thyra
379
380
381//
382// Nonmember implementations
383//
384
385template<class Scalar>
387Thyra::nonconstMultiply(
388 const RCP<LinearOpBase<Scalar> > &A,
389 const RCP<LinearOpBase<Scalar> > &B,
390 const std::string &M_label
391 )
392{
393 using Teuchos::tuple;
394 RCP<DefaultMultipliedLinearOp<Scalar> > multOp =
395 defaultMultipliedLinearOp<Scalar>(tuple(A, B)());
396 if(M_label.length())
397 multOp->setObjectLabel(M_label);
398 return multOp;
399}
400
401
402template<class Scalar>
404Thyra::multiply(
405 const RCP<const LinearOpBase<Scalar> > &A,
406 const RCP<const LinearOpBase<Scalar> > &B,
407 const std::string &M_label
408 )
409{
410 using Teuchos::tuple;
411 RCP<DefaultMultipliedLinearOp<Scalar> > multOp =
412 defaultMultipliedLinearOp<Scalar>(tuple(A, B)());
413 if(M_label.length())
414 multOp->setObjectLabel(M_label);
415 return multOp;
416}
417
418
419template<class Scalar>
421Thyra::multiply(
422 const RCP<const LinearOpBase<Scalar> > &A,
423 const RCP<const LinearOpBase<Scalar> > &B,
424 const RCP<const LinearOpBase<Scalar> > &C,
425 const std::string &M_label
426 )
427{
428 using Teuchos::tuple;
429 RCP<DefaultMultipliedLinearOp<Scalar> > multOp =
430 defaultMultipliedLinearOp<Scalar>(tuple(A, B, C)());
431 if(M_label.length())
432 multOp->setObjectLabel(M_label);
433 return multOp;
434}
435
436
437//
438// Explicit instantiation macro
439//
440// Must be expanded from within the Thyra namespace!
441//
442
443
444#define THYRA_DEFAULT_MULTIPLIED_LINEAR_OP_INSTANT(SCALAR) \
445 \
446 template class DefaultMultipliedLinearOp<SCALAR >; \
447 \
448 template RCP<LinearOpBase<SCALAR > > \
449 nonconstMultiply( \
450 const RCP<LinearOpBase<SCALAR > > &A, \
451 const RCP<LinearOpBase<SCALAR > > &B, \
452 const std::string &M_label \
453 ); \
454 \
455 template RCP<const LinearOpBase<SCALAR > > \
456 multiply( \
457 const RCP<const LinearOpBase<SCALAR > > &A, \
458 const RCP<const LinearOpBase<SCALAR > > &B, \
459 const std::string &M_label \
460 ); \
461 \
462 template RCP<const LinearOpBase<SCALAR > > \
463 multiply( \
464 const RCP<const LinearOpBase<SCALAR > > &A, \
465 const RCP<const LinearOpBase<SCALAR > > &B, \
466 const RCP<const LinearOpBase<SCALAR > > &C, \
467 const std::string &M_label \
468 ); \
469
470
471#endif // THYRA_DEFAULT_MULTIPLIED_LINEAR_OP_DEF_HPP
virtual std::string description() const
Ptr< T > ptr() const
Concrete composite LinearOpBase subclass that creates an implicitly multiplied linear operator out of...
RCP< const LinearOpBase< Scalar > > getOp(const int k) const
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Prints the details about the constituent linear operators.
RCP< const VectorSpaceBase< Scalar > > range() const
Returns this->getOp(0).range() if <t>this->numOps() > 0 and returns Teuchos::null otherwise.
void initialize(const ArrayView< const RCP< LinearOpBase< Scalar > > > &Ops)
Initialize given a list of non-const linear operators.
RCP< const VectorSpaceBase< Scalar > > domain() const
Returns this->getOp(this->numOps()-1).domain() if <t>this->numOps() > 0 and returns Teuchos::null oth...
std::string description() const
Prints just the name DefaultMultipliedLinearOp along with the overall dimensions and the number of co...
RCP< const LinearOpBase< Scalar > > clone() const
RCP< LinearOpBase< Scalar > > getNonconstOp(const int k)
bool opSupportedImpl(EOpTransp M_trans) const
Returns true only if all constituent operators support M_trans.
void applyImpl(const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha, const Scalar beta) const
Base class for all linear operators.
virtual RCP< const VectorSpaceBase< Scalar > > domain() const =0
Return a smart pointer for the domain space for this operator.
Interface for a collection of column vectors called a multi-vector.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
#define THYRA_ASSERT_LINEAR_OP_TIMES_LINEAR_OP_SPACES_NAMES(FUNC_NAME, M1, M1_T, M1_N, M2, M2_T, M2_N)
Assert that a linear operator multiplication matches up.
#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...
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. `*.
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.
std::string toString(const T &t)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)