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