Thyra  Version of the Day
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 
49 namespace Thyra {
50 
51 
52 // Inline members only used in implementation
53 
54 
55 template<class Scalar>
56 inline
57 void DefaultMultipliedLinearOp<Scalar>::assertInitialized() const
58 {
59 #ifdef TEUCHOS_DEBUG
60  TEUCHOS_TEST_FOR_EXCEPT( !( numOps() > 0 ) );
61 #endif
62 }
63 
64 
65 template<class Scalar>
66 inline
67 std::string DefaultMultipliedLinearOp<Scalar>::getClassName() const
68 {
69  return Teuchos::Describable::description();
70 }
71 
72 
73 template<class Scalar>
74 inline
75 Ordinal DefaultMultipliedLinearOp<Scalar>::getRangeDim() const
76 {
77  return (numOps() > 0 ? this->range()->dim() : 0);
78 }
79 
80 
81 template<class Scalar>
82 inline
83 Ordinal DefaultMultipliedLinearOp<Scalar>::getDomainDim() const
84 {
85  return (numOps() > 0 ? this->domain()->dim() : 0);
86 }
87 
88 
89 // Constructors/initializers/accessors
90 
91 
92 template<class Scalar>
93 DefaultMultipliedLinearOp<Scalar>::DefaultMultipliedLinearOp()
94 {}
95 
96 
97 template<class Scalar>
98 void DefaultMultipliedLinearOp<Scalar>::initialize(
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 
110 template<class Scalar>
111 void DefaultMultipliedLinearOp<Scalar>::initialize(
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 
123 template<class Scalar>
124 void DefaultMultipliedLinearOp<Scalar>::uninitialize()
125 {
126  Ops_.resize(0);
127  setupDefaultObjectLabel();
128 }
129 
130 
131 // Overridden form MultipliedLinearOpBase
132 
133 
134 template<class Scalar>
135 int DefaultMultipliedLinearOp<Scalar>::numOps() const
136 {
137  return Ops_.size();
138 }
139 
140 
141 template<class Scalar>
142 bool DefaultMultipliedLinearOp<Scalar>::opIsConst(const int k) const
143 {
144 #ifdef TEUCHOS_DEBUG
145  TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= k && k < numOps() ) );
146 #endif
147  return Ops_[k].isConst();
148 }
149 
150 
151 template<class Scalar>
152 RCP<LinearOpBase<Scalar> >
153 DefaultMultipliedLinearOp<Scalar>::getNonconstOp(const int k)
154 {
155 #ifdef TEUCHOS_DEBUG
156  TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= k && k < numOps() ) );
157 #endif
158  return Ops_[k].getNonconstObj();
159 }
160 
161 
162 template<class Scalar>
163 RCP<const LinearOpBase<Scalar> >
164 DefaultMultipliedLinearOp<Scalar>::getOp(const int k) const
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 
176 template<class Scalar>
177 RCP< const VectorSpaceBase<Scalar> >
178 DefaultMultipliedLinearOp<Scalar>::range() const
179 {
180  if (numOps()) {
181  return getOp(0)->range();
182  }
183  return Teuchos::null;
184 }
185 
186 
187 template<class Scalar>
188 RCP< const VectorSpaceBase<Scalar> >
189 DefaultMultipliedLinearOp<Scalar>::domain() const
190 {
191  if (numOps()) {
192  return getOp(numOps()-1)->domain();
193  }
194  return Teuchos::null;
195 }
196 
197 
198 template<class Scalar>
199 RCP<const LinearOpBase<Scalar> >
200 DefaultMultipliedLinearOp<Scalar>::clone() const
201 {
202  return Teuchos::null; // Not supported yet but could be!
203 }
204 
205 
206 // Overridden from Teuchos::Describable
207 
208 
209 template<class Scalar>
210 std::string DefaultMultipliedLinearOp<Scalar>::description() const
211 {
212  std::ostringstream oss;
213  oss << getClassName() << "{numOps="<<numOps()
214  <<",rangeDim=" << getRangeDim()
215  << ",domainDim="<< getDomainDim() <<"}";
216  return oss.str();
217 }
218 
219 template<class Scalar>
220 void DefaultMultipliedLinearOp<Scalar>::describe(
221  Teuchos::FancyOStream &out_arg,
222  const Teuchos::EVerbosityLevel verbLevel
223  ) const
224 {
225  using Teuchos::FancyOStream;
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) {
231  case Teuchos::VERB_DEFAULT:
232  case Teuchos::VERB_LOW:
233  *out << this->description() << std::endl;
234  break;
235  case Teuchos::VERB_MEDIUM:
236  case Teuchos::VERB_HIGH:
237  case Teuchos::VERB_EXTREME:
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 
261 template<class Scalar>
262 bool DefaultMultipliedLinearOp<Scalar>::opSupportedImpl(EOpTransp M_trans) const
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 
272 template<class Scalar>
273 void DefaultMultipliedLinearOp<Scalar>::applyImpl(
274  const EOpTransp M_trans,
275  const MultiVectorBase<Scalar> &X,
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
284  THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(
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 ) {
298  RCP<MultiVectorBase<Scalar> > Y_k;
299  RCP<const MultiVectorBase<Scalar> > X_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 ) {
317  RCP<MultiVectorBase<Scalar> > Y_k;
318  RCP<const MultiVectorBase<Scalar> > X_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 
334 template<class Scalar>
335 void DefaultMultipliedLinearOp<Scalar>::validateOps()
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 ) {
344  THYRA_ASSERT_LINEAR_OP_TIMES_LINEAR_OP_SPACES_NAMES(
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 
360 template<class Scalar>
361 void 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 
385 template<class Scalar>
386 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
387 Thyra::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 
402 template<class Scalar>
403 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
404 Thyra::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 
419 template<class Scalar>
420 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
421 Thyra::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
FancyOStream
basic_FancyOStream< char > FancyOStream
TEUCHOS_TEST_FOR_EXCEPT
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
Thyra::EOpTransp
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
Definition: Thyra_OperatorVectorTypes.hpp:160
Teuchos::Describable::description
virtual std::string description() const
Thyra::LinearOpBase::domain
virtual RCP< const VectorSpaceBase< Scalar > > domain() const =0
Return a smart pointer for the domain space for this operator.
VERB_DEFAULT
VERB_DEFAULT
Thyra::real_trans
EOpTransp real_trans(EOpTransp transp)
Return NOTRANS or TRANS for real scalar valued operators and this also is used for determining struct...
Definition: Thyra_OperatorVectorTypes.hpp:200
Thyra::DefaultMultipliedLinearOp::domain
RCP< const VectorSpaceBase< Scalar > > domain() const
Returns this->getOp(this->numOps()-1).domain() if <t>this->numOps() > 0 and returns Teuchos::null oth...
Definition: Thyra_DefaultMultipliedLinearOp_def.hpp:189
VERB_LOW
VERB_LOW
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::EVerbosityLevel
EVerbosityLevel
Teuchos::ArrayView
Thyra::DefaultMultipliedLinearOp
Concrete composite LinearOpBase subclass that creates an implicitly multiplied linear operator out of...
Definition: Thyra_DefaultMultipliedLinearOp_decl.hpp:120
Thyra::DefaultMultipliedLinearOp::clone
RCP< const LinearOpBase< Scalar > > clone() const
Definition: Thyra_DefaultMultipliedLinearOp_def.hpp:200
Teuchos::RCP
Thyra::DefaultMultipliedLinearOp::getOp
RCP< const LinearOpBase< Scalar > > getOp(const int k) const
Definition: Thyra_DefaultMultipliedLinearOp_def.hpp:164
VERB_HIGH
VERB_HIGH
Teuchos::Ptr
Thyra::Ordinal
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
Definition: Thyra_OperatorVectorTypes.hpp:127
Thyra::DefaultMultipliedLinearOp::initialize
void initialize(const ArrayView< const RCP< LinearOpBase< Scalar > > > &Ops)
Initialize given a list of non-const linear operators.
Definition: Thyra_DefaultMultipliedLinearOp_def.hpp:98
VERB_MEDIUM
VERB_MEDIUM
Teuchos::basic_FancyOStream
Thyra::DefaultMultipliedLinearOp::opSupportedImpl
bool opSupportedImpl(EOpTransp M_trans) const
Returns true only if all constituent operators support M_trans.
Definition: Thyra_DefaultMultipliedLinearOp_def.hpp:262
Thyra::MultiVectorBase
Interface for a collection of column vectors called a multi-vector.
Definition: Thyra_MultiVectorBase_decl.hpp:493
Thyra::LinearOpBase
Base class for all linear operators.
Definition: Thyra_LinearOpBase_decl.hpp:191
Thyra::DefaultMultipliedLinearOp::applyImpl
void applyImpl(const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha, const Scalar beta) const
Definition: Thyra_DefaultMultipliedLinearOp_def.hpp:273
Thyra::NOTRANS
Use the non-transposed operator.
Definition: Thyra_OperatorVectorTypes.hpp:162
Thyra::DefaultMultipliedLinearOp::DefaultMultipliedLinearOp
DefaultMultipliedLinearOp()
Constructs to uninitialized.
Definition: Thyra_DefaultMultipliedLinearOp_def.hpp:93
Thyra::DefaultMultipliedLinearOp::numOps
int numOps() const
Definition: Thyra_DefaultMultipliedLinearOp_def.hpp:135
toString
TEUCHOSCORE_LIB_DLL_EXPORT std::string toString(const EVerbosityLevel verbLevel)
Thyra::DefaultMultipliedLinearOp::range
RCP< const VectorSpaceBase< Scalar > > range() const
Returns this->getOp(0).range() if <t>this->numOps() > 0 and returns Teuchos::null otherwise.
Definition: Thyra_DefaultMultipliedLinearOp_def.hpp:178
Thyra::DefaultMultipliedLinearOp::describe
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Prints the details about the constituent linear operators.
Definition: Thyra_DefaultMultipliedLinearOp_def.hpp:220
Thyra::DefaultMultipliedLinearOp::uninitialize
void uninitialize()
Set to uninitialized.
Definition: Thyra_DefaultMultipliedLinearOp_def.hpp:124
Thyra::DefaultMultipliedLinearOp::opIsConst
bool opIsConst(const int k) const
Definition: Thyra_DefaultMultipliedLinearOp_def.hpp:142
Teuchos::RCP::ptr
Ptr< T > ptr() const
Thyra::DefaultMultipliedLinearOp::getNonconstOp
RCP< LinearOpBase< Scalar > > getNonconstOp(const int k)
Definition: Thyra_DefaultMultipliedLinearOp_def.hpp:153
VERB_EXTREME
VERB_EXTREME
THYRA_ASSERT_LINEAR_OP_TIMES_LINEAR_OP_SPACES_NAMES
#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.
Definition: Thyra_AssertOp.hpp:369
THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES
#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...
Definition: Thyra_AssertOp.hpp:250
OSTab
basic_OSTab< char > OSTab
Thyra::DefaultMultipliedLinearOp::description
std::string description() const
Prints just the name DefaultMultipliedLinearOp along with the overall dimensions and the number of co...
Definition: Thyra_DefaultMultipliedLinearOp_def.hpp:210
Generated by   doxygen 1.8.16