Thyra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Thyra_DefaultBlockedLinearOp_def.hpp
1 // @HEADER
2 // *****************************************************************************
3 // Thyra: Interfaces and Support for Abstract Numerical Algorithms
4 //
5 // Copyright 2004 NTESS and the Thyra contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef THYRA_DEFAULT_BLOCKED_LINEAR_OP_DEF_HPP
11 #define THYRA_DEFAULT_BLOCKED_LINEAR_OP_DEF_HPP
12 
13 
14 #include "Thyra_DefaultBlockedLinearOp_decl.hpp"
15 #include "Thyra_DefaultProductVectorSpace.hpp"
16 #include "Thyra_DefaultProductVector.hpp"
17 #include "Thyra_DefaultProductMultiVector.hpp"
18 #include "Thyra_MultiVectorStdOps.hpp"
19 #include "Thyra_VectorStdOps.hpp"
20 #include "Thyra_AssertOp.hpp"
21 #include "Thyra_ScaledAdjointLinearOpBase.hpp"
22 
23 
24 namespace Thyra {
25 
26 
27 // Constructors
28 
29 
30 template<class Scalar>
31 DefaultBlockedLinearOp<Scalar>::DefaultBlockedLinearOp()
32  :numRowBlocks_(0), numColBlocks_(0), blockFillIsActive_(false)
33 {}
34 
35 
36 // Overridden from PhysicallyBlockedLinearOpBase
37 
38 
39 template<class Scalar>
40 void DefaultBlockedLinearOp<Scalar>::beginBlockFill()
41 {
42  assertBlockFillIsActive(false);
43  uninitialize();
44  resetStorage(0,0);
45 }
46 
47 
48 template<class Scalar>
49 void DefaultBlockedLinearOp<Scalar>::beginBlockFill(
50  const int numRowBlocks, const int numColBlocks
51  )
52 {
53  assertBlockFillIsActive(false);
54  uninitialize();
55  resetStorage(numRowBlocks,numColBlocks);
56 }
57 
58 
59 template<class Scalar>
60 void DefaultBlockedLinearOp<Scalar>::beginBlockFill(
61  const RCP<const ProductVectorSpaceBase<Scalar> > &new_productRange
62  ,const RCP<const ProductVectorSpaceBase<Scalar> > &new_productDomain
63  )
64 {
65  using Teuchos::rcp_dynamic_cast;
66  assertBlockFillIsActive(false);
67  uninitialize();
68  productRange_ = new_productRange.assert_not_null();
69  productDomain_ = new_productDomain.assert_not_null();
70  defaultProductRange_ =
71  rcp_dynamic_cast<const DefaultProductVectorSpace<Scalar> >(productRange_);
72  defaultProductDomain_ =
73  rcp_dynamic_cast<const DefaultProductVectorSpace<Scalar> >(productDomain_);
74  // Note: the above spaces must be set before calling the next function!
75  resetStorage(productRange_->numBlocks(), productDomain_->numBlocks());
76 }
77 
78 
79 template<class Scalar>
80 bool DefaultBlockedLinearOp<Scalar>::blockFillIsActive() const
81 {
82  return blockFillIsActive_;
83 }
84 
85 
86 template<class Scalar>
87 bool DefaultBlockedLinearOp<Scalar>::acceptsBlock(
88  const int i, const int j
89  ) const
90 {
91  assertBlockFillIsActive(true);
92  assertBlockRowCol(i,j);
93  return true;
94 }
95 
96 
97 template<class Scalar>
98 void DefaultBlockedLinearOp<Scalar>::setNonconstBlock(
99  const int i, const int j
100  ,const RCP<LinearOpBase<Scalar> > &block
101  )
102 {
103  setBlockImpl(i, j, block);
104 }
105 
106 
107 template<class Scalar>
108 void DefaultBlockedLinearOp<Scalar>::setBlock(
109  const int i, const int j
110  ,const RCP<const LinearOpBase<Scalar> > &block
111  )
112 {
113  setBlockImpl(i, j, block);
114 }
115 
116 
117 template<class Scalar>
118 void DefaultBlockedLinearOp<Scalar>::endBlockFill()
119 {
120 
121  using Teuchos::as;
122 
123  assertBlockFillIsActive(true);
124 
125  // 2009/05/06: rabartl: ToDo: When doing a flexible block fill
126  // (Ops_stack_.size() > 0), we need to assert that all of the block rows and
127  // columns have been filled in. I don't think we do that here.
128 
129  // Get the number of block rows and columns
130  if (nonnull(productRange_)) {
131  numRowBlocks_ = productRange_->numBlocks();
132  numColBlocks_ = productDomain_->numBlocks();
133  }
134  else {
135  numRowBlocks_ = rangeBlocks_.size();
136  numColBlocks_ = domainBlocks_.size();
137  // NOTE: Above, whether doing a flexible fill or not, all of the blocks
138  // must be set in order to have a valid filled operator so this
139  // calculation should be correct.
140  }
141 
142  // Assert that all of the block rows and columns have at least one entry if
143  // the spaces were not given up front.
144 #ifdef TEUCHOS_DEBUG
145  if (is_null(productRange_)) {
146  for (int i = 0; i < numRowBlocks_; ++i) {
147  TEUCHOS_TEST_FOR_EXCEPTION(
148  !rangeBlocks_[i].get(), std::logic_error
149  ,"DefaultBlockedLinearOp<Scalar>::endBlockFill():"
150  " Error, no linear operator block for the i="<<i<<" block row was added"
151  " and we can not complete the block fill!"
152  );
153  }
154  for(int j = 0; j < numColBlocks_; ++j) {
155  TEUCHOS_TEST_FOR_EXCEPTION(
156  !domainBlocks_[j].get(), std::logic_error
157  ,"DefaultBlockedLinearOp<Scalar>::endBlockFill():"
158  " Error, no linear operator block for the j="
159  <<j<<" block column was added"
160  " and we can not complete the block fill!"
161  );
162  }
163  }
164 #endif
165 
166  // Insert the block LOB objects if doing a flexible fill.
167  if (Ops_stack_.size()) {
168  Ops_.resize(numRowBlocks_*numColBlocks_);
169  for ( int k = 0; k < as<int>(Ops_stack_.size()); ++k ) {
170  const BlockEntry<Scalar> &block_i_j = Ops_stack_[k];
171  Ops_[numRowBlocks_*block_i_j.j + block_i_j.i] = block_i_j.block;
172  }
173  Ops_stack_.resize(0);
174  }
175 
176  // Set the product range and domain spaces if not already set
177  if (is_null(productRange_)) {
178  adjustBlockSpaces();
179  defaultProductRange_ = productVectorSpace<Scalar>(rangeBlocks_());
180  defaultProductDomain_ = productVectorSpace<Scalar>(domainBlocks_());
181  productRange_ = defaultProductRange_;
182  productDomain_ = defaultProductDomain_;
183  }
184 
185  rangeBlocks_.resize(0);
186  domainBlocks_.resize(0);
187 
188  blockFillIsActive_ = false;
189 
190 }
191 
192 
193 template<class Scalar>
194 void DefaultBlockedLinearOp<Scalar>::uninitialize()
195 {
196  productRange_ = Teuchos::null;
197  productDomain_ = Teuchos::null;
198  numRowBlocks_ = 0;
199  numColBlocks_ = 0;
200  Ops_.resize(0);
201  Ops_stack_.resize(0);
202  rangeBlocks_.resize(0);
203  domainBlocks_.resize(0);
204  blockFillIsActive_ = false;
205 }
206 
207 
208 // Overridden from BlockedLinearOpBase
209 
210 
211 template<class Scalar>
212 RCP<const ProductVectorSpaceBase<Scalar> >
213 DefaultBlockedLinearOp<Scalar>::productRange() const
214 {
215  return productRange_;
216 }
217 
218 
219 template<class Scalar>
220 RCP<const ProductVectorSpaceBase<Scalar> >
221 DefaultBlockedLinearOp<Scalar>::productDomain() const
222 {
223  return productDomain_;
224 }
225 
226 
227 template<class Scalar>
228 bool DefaultBlockedLinearOp<Scalar>::blockExists(
229  const int i, const int j
230  ) const
231 {
232  assertBlockFillIsActive(false);
233  assertBlockRowCol(i,j);
234  return true;
235 }
236 
237 
238 template<class Scalar>
239 bool DefaultBlockedLinearOp<Scalar>::blockIsConst(
240  const int i, const int j
241  ) const
242 {
243 #ifdef TEUCHOS_DEBUG
244  TEUCHOS_TEST_FOR_EXCEPT(!blockExists(i,j));
245 #endif
246  assertBlockFillIsActive(false);
247  assertBlockRowCol(i,j);
248  return Ops_[numRowBlocks_*j+i].isConst();
249 }
250 
251 
252 template<class Scalar>
253 RCP<LinearOpBase<Scalar> >
254 DefaultBlockedLinearOp<Scalar>::getNonconstBlock(const int i, const int j)
255 {
256 #ifdef TEUCHOS_DEBUG
257  TEUCHOS_TEST_FOR_EXCEPT(!blockExists(i,j));
258 #endif
259  assertBlockFillIsActive(false);
260  assertBlockRowCol(i,j);
261  return Ops_[numRowBlocks_*j+i].getNonconstObj();
262 }
263 
264 
265 template<class Scalar>
266 RCP<const LinearOpBase<Scalar> >
267 DefaultBlockedLinearOp<Scalar>::getBlock(const int i, const int j) const
268 {
269 #ifdef TEUCHOS_DEBUG
270  TEUCHOS_TEST_FOR_EXCEPT(!blockExists(i,j));
271 #endif
272  assertBlockFillIsActive(false);
273  assertBlockRowCol(i,j);
274  return Ops_[numRowBlocks_*j+i];
275 }
276 
277 
278 // Overridden from LinearOpBase
279 
280 
281 template<class Scalar>
282 RCP< const VectorSpaceBase<Scalar> >
283 DefaultBlockedLinearOp<Scalar>::range() const
284 {
285  return productRange_;
286 }
287 
288 
289 template<class Scalar>
290 RCP< const VectorSpaceBase<Scalar> >
291 DefaultBlockedLinearOp<Scalar>::domain() const
292 {
293  return productDomain_;
294 }
295 
296 
297 template<class Scalar>
298 RCP<const LinearOpBase<Scalar> >
299 DefaultBlockedLinearOp<Scalar>::clone() const
300 {
301  return Teuchos::null; // ToDo: Implement this when needed!
302 }
303 
304 
305 // Overridden from Teuchos::Describable
306 
307 
308 template<class Scalar>
309 std::string DefaultBlockedLinearOp<Scalar>::description() const
310 {
311  assertBlockFillIsActive(false);
312  std::ostringstream oss;
313  oss
314  << Teuchos::Describable::description() << "{"
315  << "numRowBlocks="<<numRowBlocks_
316  << ",numColBlocks="<<numColBlocks_
317  << "}";
318  return oss.str();
319 }
320 
321 
322 template<class Scalar>
323 void DefaultBlockedLinearOp<Scalar>::describe(
324  Teuchos::FancyOStream &out_arg
325  ,const Teuchos::EVerbosityLevel verbLevel
326  ) const
327 {
328  using Teuchos::rcpFromRef;
329  using Teuchos::FancyOStream;
330  using Teuchos::OSTab;
331  assertBlockFillIsActive(false);
332  RCP<FancyOStream> out = rcpFromRef(out_arg);
333  OSTab tab1(out);
334  switch(verbLevel) {
335  case Teuchos::VERB_DEFAULT:
336  case Teuchos::VERB_LOW:
337  *out << this->description() << std::endl;
338  break;
339  case Teuchos::VERB_MEDIUM:
340  case Teuchos::VERB_HIGH:
341  case Teuchos::VERB_EXTREME:
342  {
343  *out
344  << Teuchos::Describable::description() << "{"
345  << "rangeDim=" << this->range()->dim()
346  << ",domainDim=" << this->domain()->dim()
347  << ",numRowBlocks=" << numRowBlocks_
348  << ",numColBlocks=" << numColBlocks_
349  << "}\n";
350  OSTab tab2(out);
351  *out
352  << "Constituent LinearOpBase objects for M = [ Op[0,0] ..."
353  << " ; ... ; ... Op[numRowBlocks-1,numColBlocks-1] ]:\n";
354  tab2.incrTab();
355  for( int i = 0; i < numRowBlocks_; ++i ) {
356  for( int j = 0; j < numColBlocks_; ++j ) {
357  *out << "Op["<<i<<","<<j<<"] = ";
358  RCP<const LinearOpBase<Scalar> >
359  block_i_j = getBlock(i,j);
360  if(block_i_j.get())
361  *out << Teuchos::describe(*getBlock(i,j),verbLevel);
362  else
363  *out << "NULL\n";
364  }
365  }
366  break;
367  }
368  default:
369  TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
370  }
371 }
372 
373 
374 // protected
375 
376 
377 // Overridden from LinearOpBase
378 
379 
380 template<class Scalar>
381 bool DefaultBlockedLinearOp<Scalar>::opSupportedImpl(EOpTransp M_trans) const
382 {
383  bool supported = true;
384  for( int i = 0; i < numRowBlocks_; ++i ) {
385  for( int j = 0; j < numColBlocks_; ++j ) {
386  RCP<const LinearOpBase<Scalar> >
387  block_i_j = getBlock(i,j);
388  if( block_i_j.get() && !Thyra::opSupported(*block_i_j,M_trans) )
389  supported = false;
390  }
391  }
392  return supported;
393 }
394 
395 
396 template<class Scalar>
397 void DefaultBlockedLinearOp<Scalar>::applyImpl(
398  const EOpTransp M_trans,
399  const MultiVectorBase<Scalar> &X_in,
400  const Ptr<MultiVectorBase<Scalar> > &Y_inout,
401  const Scalar alpha,
402  const Scalar beta
403  ) const
404 {
405 
406  using Teuchos::rcpFromRef;
407  typedef Teuchos::ScalarTraits<Scalar> ST;
408  typedef RCP<MultiVectorBase<Scalar> > MultiVectorPtr;
409  typedef RCP<const MultiVectorBase<Scalar> > ConstMultiVectorPtr;
410  typedef RCP<const LinearOpBase<Scalar> > ConstLinearOpPtr;
411 
412 #ifdef TEUCHOS_DEBUG
413  THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(
414  "DefaultBlockedLinearOp<Scalar>::apply(...)", *this, M_trans, X_in, &*Y_inout
415  );
416 #endif // TEUCHOS_DEBUG
417 
418  const bool
419  struct_transp = (real_trans(M_trans)!=NOTRANS); // Structural transpose?
420  const int
421  opNumRowBlocks = ( !struct_transp ? numRowBlocks_ : numColBlocks_ ),
422  opNumColBlocks = ( !struct_transp ? numColBlocks_ : numRowBlocks_ );
423 
424  //
425  // Y = alpha * op(M) * X + beta*Y
426  //
427  // =>
428  //
429  // Y[i] = beta+Y[i] + sum(alpha*op(Op)[i,j]*X[j],j=0...opNumColBlocks-1)
430  //
431  // , for i=0...opNumRowBlocks-1
432  //
433 
434  const RCP<const DefaultProductVectorSpace<Scalar> >
435  defaultProductRange_op = ( real_trans(M_trans)==NOTRANS
436  ? defaultProductRange_ : defaultProductDomain_ ),
437  defaultProductDomain_op = ( real_trans(M_trans)==NOTRANS
438  ? defaultProductDomain_ : defaultProductRange_ );
439 
440  const RCP<const ProductMultiVectorBase<Scalar> >
441  X = castOrCreateSingleBlockProductMultiVector<Scalar>(
442  defaultProductDomain_op, rcpFromRef(X_in));
443  const RCP<ProductMultiVectorBase<Scalar> >
444  Y = nonconstCastOrCreateSingleBlockProductMultiVector<Scalar>(
445  defaultProductRange_op, rcpFromPtr(Y_inout));
446 
447  for( int i = 0; i < opNumRowBlocks; ++i ) {
448  MultiVectorPtr Y_i = Y->getNonconstMultiVectorBlock(i);
449  for( int j = 0; j < opNumColBlocks; ++j ) {
450  ConstLinearOpPtr
451  Op_i_j = ( !struct_transp ? getBlock(i,j) : getBlock(j,i) );
452  ConstMultiVectorPtr
453  X_j = X->getMultiVectorBlock(j);
454  if(j==0) {
455  if (nonnull(Op_i_j))
456  Thyra::apply(*Op_i_j, M_trans,* X_j, Y_i.ptr(), alpha, beta);
457  else
458  scale(beta, Y_i.ptr());
459  }
460  else {
461  if (nonnull(Op_i_j))
462  Thyra::apply(*Op_i_j, M_trans, *X_j, Y_i.ptr(), alpha, ST::one());
463  }
464  }
465  }
466 }
467 
468 // Overridden from RowStatLinearOpBase
469 
470 template<class Scalar>
471 bool
472 DefaultBlockedLinearOp<Scalar>::rowStatIsSupportedImpl(
473  const RowStatLinearOpBaseUtils::ERowStat row_stat) const
474 {
475  using Teuchos::rcpFromRef;
476  using Teuchos::rcp_dynamic_cast;
477  using Teuchos::dyn_cast;
478  //typedef Teuchos::ScalarTraits<Scalar> ST; // unused
479  typedef RowStatLinearOpBase<Scalar> RowStatOp;
480  typedef RCP<const LinearOpBase<Scalar> > ConstLinearOpPtr;
481  typedef const LinearOpBase<Scalar> ConstLinearOp;
482 
483  // Figure out what needs to be check for each sub block to compute
484  // the require row stat
485  RowStatLinearOpBaseUtils::ERowStat
486  subblk_stat = RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM,
487  subblk_trans_stat = RowStatLinearOpBaseUtils::ROW_STAT_COL_SUM;
488  switch (row_stat) {
489  case RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM:
490  case RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM:
491  subblk_stat = RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM;
492  subblk_trans_stat = RowStatLinearOpBaseUtils::ROW_STAT_COL_SUM;
493  break;
494  case RowStatLinearOpBaseUtils::ROW_STAT_INV_COL_SUM:
495  case RowStatLinearOpBaseUtils::ROW_STAT_COL_SUM:
496  subblk_stat = RowStatLinearOpBaseUtils::ROW_STAT_COL_SUM;
497  subblk_trans_stat = RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM;
498  break;
499  default:
500  TEUCHOS_TEST_FOR_EXCEPT(true);
501  }
502 
503  // check each sub block for compatibility (null means zero,
504  // automatically compatible)
505  for( int i = 0; i < numRowBlocks_; ++i ) {
506  for( int j = 0; j < numColBlocks_; ++j ) {
507  ConstLinearOpPtr Op_i_j = getBlock(i,j);
508 
509  if (nonnull(Op_i_j)) {
510  // pull out adjoint, and scaling
511  ConstLinearOp * Orig_i_j = 0;
512  Scalar scalar = 1.0;
513  EOpTransp transp = NOTRANS;
514  Thyra::unwrap(*Op_i_j,&scalar,&transp,&Orig_i_j);
515 
516  const RowStatOp & row_stat_op = Teuchos::dyn_cast<const RowStatOp>(*Orig_i_j);
517 
518  // sub block must also support the required row stat operation
519  RowStatLinearOpBaseUtils::ERowStat stat = subblk_stat;
520  if(transp==NOTRANS || transp==CONJ)
521  stat = subblk_stat;
522  else if(transp==TRANS || transp==CONJTRANS)
523  stat = subblk_trans_stat;
524 
525  if(!row_stat_op.rowStatIsSupported(stat))
526  return false;
527  }
528  // else: sub block is null, "0" is good enough
529  }
530  }
531 
532  return true;
533 }
534 
535 template<class Scalar>
536 void
537 DefaultBlockedLinearOp<Scalar>::getRowStatImpl(
538  const RowStatLinearOpBaseUtils::ERowStat row_stat,
539  const Teuchos::Ptr<VectorBase< Scalar> > &rowStatVec) const
540 {
541  using Teuchos::rcpFromRef;
542  using Teuchos::rcpFromPtr;
543  using Teuchos::rcp_dynamic_cast;
544  using Teuchos::dyn_cast;
545  typedef Teuchos::ScalarTraits<Scalar> ST;
546  typedef RowStatLinearOpBase<Scalar> RowStatOp;
547  typedef RCP<const LinearOpBase<Scalar> > ConstLinearOpPtr;
548  typedef const LinearOpBase<Scalar> ConstLinearOp;
549  typedef RCP<VectorBase<Scalar> > VectorPtr;
550 
551  // Figure out what needs to be check for each sub block to compute
552  // the require row stat
553  RowStatLinearOpBaseUtils::ERowStat
554  subblk_stat = RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM,
555  subblk_trans_stat = RowStatLinearOpBaseUtils::ROW_STAT_COL_SUM;
556  switch (row_stat) {
557  case RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM:
558  case RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM:
559  subblk_stat = RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM;
560  subblk_trans_stat = RowStatLinearOpBaseUtils::ROW_STAT_COL_SUM;
561  break;
562  case RowStatLinearOpBaseUtils::ROW_STAT_INV_COL_SUM:
563  case RowStatLinearOpBaseUtils::ROW_STAT_COL_SUM:
564  subblk_stat = RowStatLinearOpBaseUtils::ROW_STAT_COL_SUM;
565  subblk_trans_stat = RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM;
566  break;
567  default:
568  TEUCHOS_TEST_FOR_EXCEPT(true);
569  }
570 
571  const RCP<ProductVectorBase<Scalar> >
572  Y = rcp_dynamic_cast<ProductVectorBase<Scalar> >(rcpFromPtr(rowStatVec));
573 
574  // using sub block row stat capability interrogate each for
575  // its constribution
576  for( int i = 0; i < numRowBlocks_; ++i ) {
577  VectorPtr blk_vec = Y->getNonconstVectorBlock(i);
578  put_scalar (ST::zero (), blk_vec.ptr ()); // clear vector just in case
579 
580  for( int j = 0; j < numColBlocks_; ++j ) {
581  ConstLinearOpPtr Op_i_j = getBlock(i,j);
582 
583  VectorPtr tmp_vec = createMember(Op_i_j->range());
584 
585  put_scalar (ST::zero (), tmp_vec.ptr ()); // clear vector just in case
586 
587  if (nonnull(Op_i_j)) {
588  // pull out adjoint, and scaling
589  ConstLinearOp * Orig_i_j = 0;
590  Scalar scalar = 1.0;
591  EOpTransp transp = NOTRANS;
592  Thyra::unwrap(*Op_i_j,&scalar,&transp,&Orig_i_j);
593 
594  const RowStatOp & row_stat_op = Teuchos::dyn_cast<const RowStatOp>(*Orig_i_j);
595 
596  // sub block must also support the required row stat operation
597  if(transp==NOTRANS || transp==CONJ)
598  row_stat_op.getRowStat(subblk_stat,tmp_vec.ptr());
599  else if(transp==TRANS || transp==CONJTRANS)
600  row_stat_op.getRowStat(subblk_trans_stat,tmp_vec.ptr());
601  else
602  { TEUCHOS_ASSERT(false); }
603 
604  // some the temporary into the block vector
605  Vp_V(blk_vec.ptr(),*tmp_vec);
606  }
607  }
608  }
609 
610  // do post processing as needed (take the inverse currently
611  switch (row_stat) {
612  case RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM:
613  case RowStatLinearOpBaseUtils::ROW_STAT_INV_COL_SUM:
614  reciprocal(*rowStatVec,rowStatVec.ptr());
615  case RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM:
616  case RowStatLinearOpBaseUtils::ROW_STAT_COL_SUM:
617  break;
618  default:
619  TEUCHOS_TEST_FOR_EXCEPT(true);
620  }
621 }
622 
623 // Overridden from ScaledLinearOpBase
624 
625 template<class Scalar>
626 bool
627 DefaultBlockedLinearOp<Scalar>::supportsScaleLeftImpl() const
628 {
629  using Teuchos::rcp_dynamic_cast;
630  using Teuchos::dyn_cast;
631  typedef RCP<const LinearOpBase<Scalar> > LinearOpPtr;
632  typedef const LinearOpBase<Scalar> LinearOp;
633  typedef const ScaledLinearOpBase<Scalar> ScaledOp;
634 
635  bool supported = true;
636  for( int i = 0; i < numRowBlocks_; ++i ) {
637  for( int j = 0; j < numColBlocks_; ++j ) {
638  LinearOpPtr Op_i_j = getBlock(i,j);
639 
640  EOpTransp transp = NOTRANS;
641 
642  if (nonnull(Op_i_j)) {
643  // pull out adjoint, and scaling
644  LinearOp * Orig_i_j = 0;
645  Scalar scalar = 1.0;
646  Thyra::unwrap(*Op_i_j,&scalar,&transp,&Orig_i_j);
647 
648  ScaledOp & scaled_op = Teuchos::dyn_cast<ScaledOp>(*Orig_i_j);
649 
650  if(transp==NOTRANS || transp==CONJ)
651  supported &= scaled_op.supportsScaleLeft();
652  else if(transp==TRANS || transp==CONJTRANS)
653  supported &= scaled_op.supportsScaleRight();
654  }
655  }
656  }
657 
658  return supported;
659 }
660 
661 template<class Scalar>
662 bool
663 DefaultBlockedLinearOp<Scalar>::supportsScaleRightImpl() const
664 {
665  using Teuchos::rcp_dynamic_cast;
666  using Teuchos::dyn_cast;
667  typedef RCP<const LinearOpBase<Scalar> > LinearOpPtr;
668  typedef const LinearOpBase<Scalar> LinearOp;
669  typedef const ScaledLinearOpBase<Scalar> ScaledOp;
670 
671  bool supported = true;
672  for( int i = 0; i < numRowBlocks_; ++i ) {
673  for( int j = 0; j < numColBlocks_; ++j ) {
674  LinearOpPtr Op_i_j = getBlock(i,j);
675 
676  if (nonnull(Op_i_j)) {
677  // pull out adjoint, and scaling
678  LinearOp * Orig_i_j = 0;
679  Scalar scalar = 1.0;
680  EOpTransp transp = NOTRANS;
681  Thyra::unwrap(*Op_i_j,&scalar,&transp,&Orig_i_j);
682 
683  ScaledOp & scaled_op = Teuchos::dyn_cast<ScaledOp>(*Orig_i_j);
684 
685  if(transp==NOTRANS || transp==CONJ)
686  supported &= scaled_op.supportsScaleRight();
687  else if(transp==TRANS || transp==CONJTRANS)
688  supported &= scaled_op.supportsScaleLeft();
689  }
690  }
691  }
692 
693  return supported;
694 }
695 
696 template<class Scalar>
697 void
698 DefaultBlockedLinearOp<Scalar>::scaleLeftImpl(
699  const VectorBase< Scalar > &row_scaling)
700 {
701  using Teuchos::dyn_cast;
702  using Teuchos::rcp_dynamic_cast;
703  typedef ScaledLinearOpBase<Scalar> ScaledOp;
704  typedef RCP<LinearOpBase<Scalar> > LinearOpPtr;
705  typedef RCP<const VectorBase<Scalar> > VectorPtr;
706 
707  const ProductVectorBase<Scalar> &
708  Y = dyn_cast<const ProductVectorBase<Scalar> >(row_scaling);
709 
710  // using sub block row stat capability interrogate each for
711  // its constribution
712  for( int i = 0; i < numRowBlocks_; ++i ) {
713  VectorPtr blk_vec = Y.getVectorBlock(i);
714 
715  for( int j = 0; j < numColBlocks_; ++j ) {
716  LinearOpPtr Op_i_j = getNonconstBlock(i,j);
717 
718  if (nonnull(Op_i_j)) {
719  // pull out adjoint, and scaling
720  LinearOpPtr Orig_i_j;
721  EOpTransp transp = NOTRANS;
722 
723  {
724  RCP<ScaledAdjointLinearOpBase<Scalar> > saOp
725  = rcp_dynamic_cast<ScaledAdjointLinearOpBase<Scalar> >(Op_i_j);
726  if(nonnull(saOp)) {
727  transp = saOp->overallTransp();
728  Orig_i_j = saOp->getNonconstOrigOp();
729  }
730  else
731  Orig_i_j = Op_i_j; // stick with the original
732  }
733 
734  RCP<ScaledOp> scaled_op = rcp_dynamic_cast<ScaledOp>(Orig_i_j);
735 
736  // sub block must support a row stat operation
737  TEUCHOS_ASSERT(scaled_op!=Teuchos::null); // guranteed by compatibility check
738 
739  // sub block must also support the required row stat operation
740  if(transp==NOTRANS || transp==CONJ)
741  scaled_op->scaleLeft(*blk_vec);
742  else if(transp==TRANS || transp==CONJTRANS)
743  scaled_op->scaleRight(*blk_vec);
744  }
745  }
746  }
747 }
748 
749 template<class Scalar>
750 void
751 DefaultBlockedLinearOp<Scalar>::scaleRightImpl(
752  const VectorBase< Scalar > &col_scaling)
753 {
754  using Teuchos::dyn_cast;
755  using Teuchos::rcp_dynamic_cast;
756  typedef ScaledLinearOpBase<Scalar> ScaledOp;
757  typedef RCP<LinearOpBase<Scalar> > LinearOpPtr;
758  typedef RCP<const VectorBase<Scalar> > VectorPtr;
759 
760  const ProductVectorBase<Scalar> &
761  Y = dyn_cast<const ProductVectorBase<Scalar> >(col_scaling);
762 
763  // using sub block row stat capability interrogate each for
764  // its constribution
765  for( int j = 0; j < numColBlocks_; ++j ) {
766  VectorPtr blk_vec = Y.getVectorBlock(j);
767 
768  for( int i = 0; i < numRowBlocks_; ++i ) {
769  LinearOpPtr Op_i_j = getNonconstBlock(i,j);
770 
771  if (nonnull(Op_i_j)) {
772  // pull out adjoint, and scaling
773  LinearOpPtr Orig_i_j;
774  EOpTransp transp = NOTRANS;
775 
776  {
777  RCP<ScaledAdjointLinearOpBase<Scalar> > saOp
778  = rcp_dynamic_cast<ScaledAdjointLinearOpBase<Scalar> >(Op_i_j);
779  if(nonnull(saOp)) {
780  transp = saOp->overallTransp();
781  Orig_i_j = saOp->getNonconstOrigOp();
782  }
783  else
784  Orig_i_j = Op_i_j; // stick with the original
785  }
786 
787  RCP<ScaledOp> scaled_op = rcp_dynamic_cast<ScaledOp>(Orig_i_j);
788 
789  // sub block must support a row stat operation
790  TEUCHOS_ASSERT(scaled_op!=Teuchos::null); // guranteed by compatibility check
791 
792  // sub block must also support the required row stat operation
793  if(transp==NOTRANS || transp==CONJ)
794  scaled_op->scaleRight(*blk_vec);
795  else if(transp==TRANS || transp==CONJTRANS)
796  scaled_op->scaleLeft(*blk_vec);
797  }
798  }
799  }
800 }
801 
802 // private
803 
804 
805 template<class Scalar>
806 void DefaultBlockedLinearOp<Scalar>::resetStorage(
807  const int numRowBlocks, const int numColBlocks
808  )
809 {
810  numRowBlocks_ = numRowBlocks;
811  numColBlocks_ = numColBlocks;
812  Ops_.resize(numRowBlocks_*numColBlocks_);
813  if (is_null(productRange_)) {
814  rangeBlocks_.resize(numRowBlocks);
815  domainBlocks_.resize(numColBlocks);
816  }
817  blockFillIsActive_ = true;
818 }
819 
820 
821 template<class Scalar>
822 void DefaultBlockedLinearOp<Scalar>::assertBlockFillIsActive(
823  bool wantedValue
824  ) const
825 {
826 #ifdef TEUCHOS_DEBUG
827  TEUCHOS_TEST_FOR_EXCEPT(!(blockFillIsActive_==wantedValue));
828 #else
829  (void)wantedValue;
830 #endif
831 }
832 
833 
834 template<class Scalar>
835 void DefaultBlockedLinearOp<Scalar>::assertBlockRowCol(
836  const int i, const int j
837  ) const
838 {
839 #ifdef TEUCHOS_DEBUG
840  TEUCHOS_TEST_FOR_EXCEPTION(
841  !( 0 <= i ), std::logic_error
842  ,"Error, i="<<i<<" is invalid!"
843  );
844  TEUCHOS_TEST_FOR_EXCEPTION(
845  !( 0 <= j ), std::logic_error
846  ,"Error, j="<<j<<" is invalid!"
847  );
848  // Only validate upper range if the number of row and column blocks is
849  // fixed!
850  if(Ops_.size()) {
851  TEUCHOS_TEST_FOR_EXCEPTION(
852  !( 0 <= i && i < numRowBlocks_ ), std::logic_error
853  ,"Error, i="<<i<<" does not fall in the range [0,"<<numRowBlocks_-1<<"]!"
854  );
855  TEUCHOS_TEST_FOR_EXCEPTION(
856  !( 0 <= j && j < numColBlocks_ ), std::logic_error
857  ,"Error, j="<<j<<" does not fall in the range [0,"<<numColBlocks_-1<<"]!"
858  );
859  }
860 #else
861  (void)i;
862  (void)j;
863 #endif
864 }
865 
866 
867 template<class Scalar>
868 void DefaultBlockedLinearOp<Scalar>::setBlockSpaces(
869  const int i, const int j, const LinearOpBase<Scalar> &block
870  )
871 {
872  using Teuchos::toString;
873  assertBlockFillIsActive(true);
874  assertBlockRowCol(i,j);
875 
876  // Validate that if the vector space block is already set that it is
877  // compatible with the block that is being set.
878  if( i < numRowBlocks_ && j < numColBlocks_ ) {
879 #ifdef TEUCHOS_DEBUG
880  RCP<const VectorSpaceBase<Scalar> >
881  rangeBlock = (
882  productRange_.get()
883  ? productRange_->getBlock(i)
884  : rangeBlocks_[i]
885  ),
886  domainBlock = (
887  productDomain_.get()
888  ? productDomain_->getBlock(j)
889  : domainBlocks_[j]
890  );
891  if(rangeBlock.get()) {
892  THYRA_ASSERT_VEC_SPACES_NAMES(
893  "DefaultBlockedLinearOp<Scalar>::setBlockSpaces(i,j,block):\n\n"
894  "Adding block: " + block.description(),
895  *rangeBlock,("(*productRange->getBlock("+toString(i)+"))"),
896  *block.range(),("(*block["+toString(i)+","+toString(j)+"].range())")
897  );
898  }
899  if(domainBlock.get()) {
900  THYRA_ASSERT_VEC_SPACES_NAMES(
901  "DefaultBlockedLinearOp<Scalar>::setBlockSpaces(i,j,block):\n\n"
902  "Adding block: " + block.description(),
903  *domainBlock,("(*productDomain->getBlock("+toString(j)+"))"),
904  *block.domain(),("(*block["+toString(i)+","+toString(j)+"].domain())")
905  );
906  }
907 #endif // TEUCHOS_DEBUG
908  }
909 
910  // Add spaces missing range and domain space blocks if we are doing a
911  // flexible fill (otherwise these loops will not be executed)
912  for( int k = numRowBlocks_; k <= i; ++k )
913  rangeBlocks_.push_back(Teuchos::null);
914  for( int k = numColBlocks_; k <= j; ++k )
915  domainBlocks_.push_back(Teuchos::null);
916 
917  // Set the incoming range and domain blocks if not already set
918  if(!productRange_.get()) {
919  if(!rangeBlocks_[i].get())
920  rangeBlocks_[i] = block.range().assert_not_null();
921  if(!domainBlocks_[j].get()) {
922  domainBlocks_[j] = block.domain().assert_not_null();
923  }
924  }
925 
926  // Update the current number of row and columns blocks if doing a flexible
927  // fill.
928  if(!Ops_.size()) {
929  numRowBlocks_ = rangeBlocks_.size();
930  numColBlocks_ = domainBlocks_.size();
931  }
932 
933 }
934 
935 
936 template<class Scalar>
937 template<class LinearOpType>
938 void DefaultBlockedLinearOp<Scalar>::setBlockImpl(
939  const int i, const int j,
940  const RCP<LinearOpType> &block
941  )
942 {
943  setBlockSpaces(i, j, *block);
944  if (Ops_.size()) {
945  // We are doing a fill with a fixed number of row and column blocks so we
946  // can just set this.
947  Ops_[numRowBlocks_*j+i] = block;
948  }
949  else {
950  // We are doing a flexible fill so add the block to the stack of blocks or
951  // replace a block that already exists.
952  bool foundBlock = false;
953  for( unsigned int k = 0; k < Ops_stack_.size(); ++k ) {
954  BlockEntry<Scalar> &block_i_j = Ops_stack_[k];
955  if( block_i_j.i == i && block_i_j.j == j ) {
956  block_i_j.block = block;
957  foundBlock = true;
958  break;
959  }
960  }
961  if(!foundBlock)
962  Ops_stack_.push_back(BlockEntry<Scalar>(i,j,block));
963  }
964 }
965 
966 
967 template<class Scalar>
968 void DefaultBlockedLinearOp<Scalar>::adjustBlockSpaces()
969 {
970 
971 #ifdef TEUCHOS_DEBUG
972  TEUCHOS_ASSERT_INEQUALITY(Ops_.size(), !=, 0);
973 #endif
974 
975  //
976  // Loop through the rows and columns looking for rows with mixed
977  // single-space range and/or domain spaces on operators and set the single
978  // spaces as the block space if it exists.
979  //
980  // NOTE: Once we get here, we can safely assume that all of the operators
981  // are compatible w.r.t. their spaces so if there are rows and/or columns
982  // with mixed product and single vector spaces that we can just pick the
983  // single vector space for the whole row and/or column.
984  //
985 
986  // Adjust blocks in the range space
987  for (int i = 0; i < numRowBlocks_; ++i) {
988  for (int j = 0; j < numColBlocks_; ++j) {
989  const RCP<const LinearOpBase<Scalar> >
990  op_i_j = Ops_[numRowBlocks_*j+i];
991  if (is_null(op_i_j))
992  continue;
993  const RCP<const VectorSpaceBase<Scalar> > range_i_j = op_i_j->range();
994  if (is_null(productVectorSpaceBase<Scalar>(range_i_j, false))) {
995  rangeBlocks_[i] = range_i_j;
996  break;
997  }
998  }
999  }
1000 
1001  // Adjust blocks in the domain space
1002  for (int j = 0; j < numColBlocks_; ++j) {
1003  for (int i = 0; i < numRowBlocks_; ++i) {
1004  const RCP<const LinearOpBase<Scalar> >
1005  op_i_j = Ops_[numRowBlocks_*j+i];
1006  if (is_null(op_i_j))
1007  continue;
1008  const RCP<const VectorSpaceBase<Scalar> >
1009  domain_i_j = op_i_j->domain();
1010  if (is_null(productVectorSpaceBase<Scalar>(domain_i_j, false))) {
1011  domainBlocks_[j] = domain_i_j;
1012  break;
1013  }
1014  }
1015  }
1016 
1017 }
1018 
1019 
1020 } // namespace Thyra
1021 
1022 
1023 template<class Scalar>
1024 Teuchos::RCP<Thyra::DefaultBlockedLinearOp<Scalar> >
1025 Thyra::defaultBlockedLinearOp()
1026 {
1027  return Teuchos::rcp(new DefaultBlockedLinearOp<Scalar>());
1028 }
1029 
1030 
1031 template<class Scalar>
1032 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
1033 Thyra::block1x1(
1034  const RCP<const LinearOpBase<Scalar> > &A00,
1035  const std::string &label
1036  )
1037 {
1038  RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1039  M = defaultBlockedLinearOp<Scalar>();
1040  M->beginBlockFill(1,1);
1041  M->setBlock(0, 0, A00);
1042  M->endBlockFill();
1043  if (label.length())
1044  M->setObjectLabel(label);
1045  return M;
1046 }
1047 
1048 
1049 template<class Scalar>
1050 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
1051 Thyra::block1x2(
1052  const RCP<const LinearOpBase<Scalar> > &A00,
1053  const RCP<const LinearOpBase<Scalar> > &A01,
1054  const std::string &label
1055  )
1056 {
1057  RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1058  M = defaultBlockedLinearOp<Scalar>();
1059  M->beginBlockFill(1,2);
1060  if(A00.get()) M->setBlock(0,0,A00);
1061  if(A01.get()) M->setBlock(0,1,A01);
1062  M->endBlockFill();
1063  if (label.length())
1064  M->setObjectLabel(label);
1065  return M;
1066 }
1067 
1068 
1069 template<class Scalar>
1070 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
1071 Thyra::block2x1(
1072  const RCP<const LinearOpBase<Scalar> > &A00,
1073  const RCP<const LinearOpBase<Scalar> > &A10,
1074  const std::string &label
1075  )
1076 {
1077  RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1078  M = defaultBlockedLinearOp<Scalar>();
1079  M->beginBlockFill(2,1);
1080  if(A00.get()) M->setBlock(0,0,A00);
1081  if(A10.get()) M->setBlock(1,0,A10);
1082  M->endBlockFill();
1083  if (label.length())
1084  M->setObjectLabel(label);
1085  return M;
1086 }
1087 
1088 
1089 template<class Scalar>
1090 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
1091 Thyra::block2x2(
1092  const RCP<const LinearOpBase<Scalar> > &A00,
1093  const RCP<const LinearOpBase<Scalar> > &A01,
1094  const RCP<const LinearOpBase<Scalar> > &A10,
1095  const RCP<const LinearOpBase<Scalar> > &A11,
1096  const std::string &label
1097  )
1098 {
1099  RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1100  M = defaultBlockedLinearOp<Scalar>();
1101  M->beginBlockFill(2,2);
1102  if(A00.get()) M->setBlock(0,0,A00);
1103  if(A01.get()) M->setBlock(0,1,A01);
1104  if(A10.get()) M->setBlock(1,0,A10);
1105  if(A11.get()) M->setBlock(1,1,A11);
1106  M->endBlockFill();
1107  if (label.length())
1108  M->setObjectLabel(label);
1109  return M;
1110 }
1111 
1112 
1113 template<class Scalar>
1114 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
1115 Thyra::nonconstBlock1x1(
1116  const RCP<LinearOpBase<Scalar> > &A00,
1117  const std::string &label
1118  )
1119 {
1120  RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1121  M = defaultBlockedLinearOp<Scalar>();
1122  M->beginBlockFill(1, 1);
1123  M->setNonconstBlock(0, 0, A00);
1124  M->endBlockFill();
1125  if (label.length())
1126  M->setObjectLabel(label);
1127  return M;
1128 }
1129 
1130 
1131 template<class Scalar>
1132 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
1133 Thyra::nonconstBlock1x2(
1134  const RCP<LinearOpBase<Scalar> > &A00,
1135  const RCP<LinearOpBase<Scalar> > &A01,
1136  const std::string &label
1137  )
1138 {
1139  RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1140  M = defaultBlockedLinearOp<Scalar>();
1141  M->beginBlockFill(1,2);
1142  if(A00.get()) M->setNonconstBlock(0,0,A00);
1143  if(A01.get()) M->setNonconstBlock(0,1,A01);
1144  M->endBlockFill();
1145  if (label.length())
1146  M->setObjectLabel(label);
1147  return M;
1148 }
1149 
1150 
1151 template<class Scalar>
1152 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
1153 Thyra::nonconstBlock2x1(
1154  const RCP<LinearOpBase<Scalar> > &A00,
1155  const RCP<LinearOpBase<Scalar> > &A10,
1156  const std::string &label
1157  )
1158 {
1159  RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1160  M = defaultBlockedLinearOp<Scalar>();
1161  M->beginBlockFill(2,1);
1162  if(A00.get()) M->setNonconstBlock(0,0,A00);
1163  if(A10.get()) M->setNonconstBlock(1,0,A10);
1164  M->endBlockFill();
1165  if (label.length())
1166  M->setObjectLabel(label);
1167  return M;
1168 }
1169 
1170 
1171 template<class Scalar>
1172 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
1173 Thyra::nonconstBlock2x2(
1174  const RCP<LinearOpBase<Scalar> > &A00,
1175  const RCP<LinearOpBase<Scalar> > &A01,
1176  const RCP<LinearOpBase<Scalar> > &A10,
1177  const RCP<LinearOpBase<Scalar> > &A11,
1178  const std::string &label
1179  )
1180 {
1181  RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1182  M = defaultBlockedLinearOp<Scalar>();
1183  M->beginBlockFill(2,2);
1184  if(A00.get()) M->setNonconstBlock(0,0,A00);
1185  if(A01.get()) M->setNonconstBlock(0,1,A01);
1186  if(A10.get()) M->setNonconstBlock(1,0,A10);
1187  if(A11.get()) M->setNonconstBlock(1,1,A11);
1188  M->endBlockFill();
1189  if (label.length())
1190  M->setObjectLabel(label);
1191  return M;
1192 }
1193 
1194 
1195 //
1196 // Explicit instantiation macro
1197 //
1198 // Must be expanded from within the Thyra namespace!
1199 //
1200 
1201 
1202 #define THYRA_DEFAULT_BLOCKED_LINEAR_OP_INSTANT(SCALAR) \
1203  \
1204  template class DefaultBlockedLinearOp<SCALAR >; \
1205  \
1206  template RCP<DefaultBlockedLinearOp<SCALAR > > \
1207  defaultBlockedLinearOp<SCALAR >(); \
1208  \
1209  template RCP<const LinearOpBase<SCALAR > > \
1210  block1x1( \
1211  const RCP<const LinearOpBase<SCALAR > > &A00, \
1212  const std::string &label \
1213  ); \
1214  \
1215  template RCP<const LinearOpBase<SCALAR > > \
1216  block1x2( \
1217  const RCP<const LinearOpBase<SCALAR > > &A00, \
1218  const RCP<const LinearOpBase<SCALAR > > &A01, \
1219  const std::string &label \
1220  ); \
1221  \
1222  template RCP<const LinearOpBase<SCALAR > > \
1223  block2x1( \
1224  const RCP<const LinearOpBase<SCALAR > > &A00, \
1225  const RCP<const LinearOpBase<SCALAR > > &A10, \
1226  const std::string &label \
1227  ); \
1228  \
1229  template RCP<const LinearOpBase<SCALAR > > \
1230  block2x2( \
1231  const RCP<const LinearOpBase<SCALAR > > &A00, \
1232  const RCP<const LinearOpBase<SCALAR > > &A01, \
1233  const RCP<const LinearOpBase<SCALAR > > &A10, \
1234  const RCP<const LinearOpBase<SCALAR > > &A11, \
1235  const std::string &label \
1236  ); \
1237  \
1238  template RCP<LinearOpBase<SCALAR > > \
1239  nonconstBlock1x1( \
1240  const RCP<LinearOpBase<SCALAR > > &A00, \
1241  const std::string &label \
1242  ); \
1243  \
1244  template RCP<LinearOpBase<SCALAR > > \
1245  nonconstBlock1x2( \
1246  const RCP<LinearOpBase<SCALAR > > &A00, \
1247  const RCP<LinearOpBase<SCALAR > > &A01, \
1248  const std::string &label \
1249  ); \
1250  \
1251  template RCP<LinearOpBase<SCALAR > > \
1252  nonconstBlock2x1( \
1253  const RCP<LinearOpBase<SCALAR > > &A00, \
1254  const RCP<LinearOpBase<SCALAR > > &A10, \
1255  const std::string &label \
1256  ); \
1257  \
1258  template RCP<LinearOpBase<SCALAR > > \
1259  nonconstBlock2x2( \
1260  const RCP<LinearOpBase<SCALAR > > &A00, \
1261  const RCP<LinearOpBase<SCALAR > > &A01, \
1262  const RCP<LinearOpBase<SCALAR > > &A10, \
1263  const RCP<LinearOpBase<SCALAR > > &A11, \
1264  const std::string &label \
1265  );
1266 
1267 
1268 #endif // THYRA_DEFAULT_BLOCKED_LINEAR_OP_DEF_HPP
Thyra::DefaultBlockedLinearOp::blockExists
bool blockExists(const int i, const int j) const
Definition: Thyra_DefaultBlockedLinearOp_def.hpp:228
Teuchos::VERB_DEFAULT
is_null
bool is_null(const boost::shared_ptr< T > &p)
Thyra::EOpTransp
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
Definition: Thyra_OperatorVectorTypes.hpp:128
Thyra::ProductVectorBase
Base interface for product vectors.
Definition: Thyra_ProductVectorBase.hpp:55
Teuchos::OSTab
basic_OSTab< char > OSTab
Thyra::DefaultBlockedLinearOp::productRange
Teuchos::RCP< const ProductVectorSpaceBase< Scalar > > productRange() const
Definition: Thyra_DefaultBlockedLinearOp_def.hpp:213
Thyra::unwrap
void unwrap(const LinearOpBase< Scalar > &Op, Scalar *scalar, EOpTransp *transp, const LinearOpBase< Scalar > **origOp)
Extract the overallScalar, overallTransp and const origOp from a const LinearOpBase object...
Definition: Thyra_ScaledAdjointLinearOpBase_def.hpp:18
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:218
Thyra::ProductVectorBase::getVectorBlock
virtual RCP< const VectorBase< Scalar > > getVectorBlock(const int k) const =0
Returns a non-persisting const view of the (zero-based) kth block vector.
Teuchos::VERB_LOW
Thyra::DefaultBlockedLinearOp::acceptsBlock
bool acceptsBlock(const int i, const int j) const
Definition: Thyra_DefaultBlockedLinearOp_def.hpp:87
Thyra::DefaultBlockedLinearOp::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_DefaultBlockedLinearOp_def.hpp:397
TEUCHOS_ASSERT_INEQUALITY
#define TEUCHOS_ASSERT_INEQUALITY(val1, comp, val2)
Teuchos::FancyOStream
basic_FancyOStream< char > FancyOStream
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Thyra::NOTRANS
Use the non-transposed operator.
Definition: Thyra_OperatorVectorTypes.hpp:130
Teuchos::VERB_EXTREME
Teuchos::dyn_cast
T_To & dyn_cast(T_From &from)
Thyra::DefaultBlockedLinearOp::describe
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Prints the details about the constituent linear operators.
Definition: Thyra_DefaultBlockedLinearOp_def.hpp:323
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:168
Thyra::CONJTRANS
Use the transposed operator with complex-conjugate clements (same as TRANS for real scalar types)...
Definition: Thyra_OperatorVectorTypes.hpp:140
Teuchos::RCP::get
T * get() const
Teuchos::EVerbosityLevel
EVerbosityLevel
Thyra::CONJ
Use the non-transposed operator with complex-conjugate elements (same as NOTRANS for real scalar type...
Definition: Thyra_OperatorVectorTypes.hpp:134
Teuchos::VERB_MEDIUM
Thyra::TRANS
Use the transposed operator.
Definition: Thyra_OperatorVectorTypes.hpp:136
Thyra::DefaultBlockedLinearOp::getRowStatImpl
virtual void getRowStatImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat, const Teuchos::Ptr< VectorBase< Scalar > > &rowStatVec) const
Definition: Thyra_DefaultBlockedLinearOp_def.hpp:537
Thyra::DefaultBlockedLinearOp::blockIsConst
bool blockIsConst(const int i, const int j) const
Definition: Thyra_DefaultBlockedLinearOp_def.hpp:239
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::VERB_HIGH
Thyra::DefaultBlockedLinearOp::setNonconstBlock
void setNonconstBlock(const int i, const int j, const Teuchos::RCP< LinearOpBase< Scalar > > &block)
Definition: Thyra_DefaultBlockedLinearOp_def.hpp:98
Thyra::MultiVectorBase
Interface for a collection of column vectors called a multi-vector.
Definition: Thyra_MultiVectorBase_decl.hpp:461
Teuchos::basic_FancyOStream
Thyra::ScaledAdjointLinearOpBase
Base class for LinearOpBase decorator subclasses that wrap a LinearOpBase object and adds on an extra...
Definition: Thyra_ScaledAdjointLinearOpBase_decl.hpp:41
Thyra::DefaultBlockedLinearOp::description
std::string description() const
Prints just the name DefaultBlockedLinearOp along with the overall dimensions and the number of const...
Definition: Thyra_DefaultBlockedLinearOp_def.hpp:309
Teuchos::Describable::description
virtual std::string description() const
Thyra::VectorBase
Abstract interface for finite-dimensional dense vectors.
Definition: Thyra_OperatorVectorTypes.hpp:338
toString
TEUCHOSCORE_LIB_DLL_EXPORT std::string toString(const EVerbosityLevel verbLevel)
Thyra::DefaultBlockedLinearOp::clone
Teuchos::RCP< const LinearOpBase< Scalar > > clone() const
Definition: Thyra_DefaultBlockedLinearOp_def.hpp:299
Thyra::DefaultBlockedLinearOp::domain
Teuchos::RCP< const VectorSpaceBase< Scalar > > domain() const
Definition: Thyra_DefaultBlockedLinearOp_def.hpp:291
Thyra::LinearOpBase
Base class for all linear operators.
Definition: Thyra_LinearOpBase_decl.hpp:159
Thyra::ScaledLinearOpBase
Applies left or right sclaing to the linear operator.
Definition: Thyra_ScaledLinearOpBase.hpp:50
Thyra::DefaultBlockedLinearOp::productDomain
Teuchos::RCP< const ProductVectorSpaceBase< Scalar > > productDomain() const
Definition: Thyra_DefaultBlockedLinearOp_def.hpp:221
Thyra::DefaultBlockedLinearOp::range
Teuchos::RCP< const VectorSpaceBase< Scalar > > range() const
Definition: Thyra_DefaultBlockedLinearOp_def.hpp:283
nonnull
bool nonnull(const boost::shared_ptr< T > &p)
Thyra::DefaultBlockedLinearOp::scaleLeftImpl
virtual void scaleLeftImpl(const VectorBase< Scalar > &row_scaling)
Definition: Thyra_DefaultBlockedLinearOp_def.hpp:698
Teuchos::as
TypeTo as(const TypeFrom &t)
Thyra::DefaultBlockedLinearOp::scaleRightImpl
virtual void scaleRightImpl(const VectorBase< Scalar > &col_scaling)
Definition: Thyra_DefaultBlockedLinearOp_def.hpp:751
Thyra::DefaultBlockedLinearOp::opSupportedImpl
bool opSupportedImpl(EOpTransp M_trans) const
Returns true only if all constituent operators support M_trans.
Definition: Thyra_DefaultBlockedLinearOp_def.hpp:381
Thyra::RowStatLinearOpBase
Interface for exxtracting row statistics as a VectorBase from a supporting LinearOpBase object...
Definition: Thyra_RowStatLinearOpBase.hpp:44
THYRA_ASSERT_VEC_SPACES_NAMES
#define THYRA_ASSERT_VEC_SPACES_NAMES(FUNC_NAME, VS1, VS1_NAME, VS2, VS2_NAME)
Helper assertion macro.
Definition: Thyra_AssertOp.hpp:134
Thyra::DefaultBlockedLinearOp::rowStatIsSupportedImpl
virtual bool rowStatIsSupportedImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat) const
Definition: Thyra_DefaultBlockedLinearOp_def.hpp:472
Thyra::DefaultBlockedLinearOp::getNonconstBlock
Teuchos::RCP< LinearOpBase< Scalar > > getNonconstBlock(const int i, const int j)
Definition: Thyra_DefaultBlockedLinearOp_def.hpp:254
TEUCHOS_ASSERT
#define TEUCHOS_ASSERT(assertion_test)
Thyra::DefaultBlockedLinearOp::setBlock
void setBlock(const int i, const int j, const Teuchos::RCP< const LinearOpBase< Scalar > > &block)
Definition: Thyra_DefaultBlockedLinearOp_def.hpp:108
Thyra::DefaultBlockedLinearOp
Concrete composite LinearOpBase subclass that creates single linear operator object out of a set of c...
Definition: Thyra_DefaultBlockedLinearOp_decl.hpp:54
TEUCHOS_TEST_FOR_EXCEPT
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
Teuchos::toString
std::string toString(const T &t)
Thyra::DefaultBlockedLinearOp::getBlock
Teuchos::RCP< const LinearOpBase< Scalar > > getBlock(const int i, const int j) const
Definition: Thyra_DefaultBlockedLinearOp_def.hpp:267
Generated on Thu Nov 21 2024 09:20:50 for Thyra by   doxygen 1.8.5