Thyra Version of the Day
Loading...
Searching...
No Matches
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
24namespace Thyra {
25
26
27// Constructors
28
29
30template<class Scalar>
32 :numRowBlocks_(0), numColBlocks_(0), blockFillIsActive_(false)
33{}
34
35
36// Overridden from PhysicallyBlockedLinearOpBase
37
38
39template<class Scalar>
41{
42 assertBlockFillIsActive(false);
43 uninitialize();
44 resetStorage(0,0);
45}
46
47
48template<class Scalar>
50 const int numRowBlocks, const int numColBlocks
51 )
52{
53 assertBlockFillIsActive(false);
54 uninitialize();
55 resetStorage(numRowBlocks,numColBlocks);
56}
57
58
59template<class Scalar>
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
79template<class Scalar>
81{
82 return blockFillIsActive_;
83}
84
85
86template<class Scalar>
88 const int i, const int j
89 ) const
90{
91 assertBlockFillIsActive(true);
92 assertBlockRowCol(i,j);
93 return true;
94}
95
96
97template<class Scalar>
99 const int i, const int j
100 ,const RCP<LinearOpBase<Scalar> > &block
101 )
102{
103 setBlockImpl(i, j, block);
104}
105
106
107template<class Scalar>
109 const int i, const int j
110 ,const RCP<const LinearOpBase<Scalar> > &block
111 )
112{
113 setBlockImpl(i, j, block);
114}
115
116
117template<class Scalar>
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) {
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) {
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
193template<class Scalar>
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
211template<class Scalar>
214{
215 return productRange_;
216}
217
218
219template<class Scalar>
222{
223 return productDomain_;
224}
225
226
227template<class Scalar>
229 const int i, const int j
230 ) const
231{
232 assertBlockFillIsActive(false);
233 assertBlockRowCol(i,j);
234 return true;
235}
236
237
238template<class Scalar>
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
252template<class Scalar>
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
265template<class Scalar>
267DefaultBlockedLinearOp<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
281template<class Scalar>
284{
285 return productRange_;
286}
287
288
289template<class Scalar>
292{
293 return productDomain_;
294}
295
296
297template<class Scalar>
300{
301 return Teuchos::null; // ToDo: Implement this when needed!
302}
303
304
305// Overridden from Teuchos::Describable
306
307
308template<class Scalar>
310{
311 assertBlockFillIsActive(false);
312 std::ostringstream oss;
313 oss
315 << "numRowBlocks="<<numRowBlocks_
316 << ",numColBlocks="<<numColBlocks_
317 << "}";
318 return oss.str();
319}
320
321
322template<class Scalar>
324 Teuchos::FancyOStream &out_arg
325 ,const Teuchos::EVerbosityLevel verbLevel
326 ) const
327{
328 using Teuchos::rcpFromRef;
330 using Teuchos::OSTab;
331 assertBlockFillIsActive(false);
332 RCP<FancyOStream> out = rcpFromRef(out_arg);
333 OSTab tab1(out);
334 switch(verbLevel) {
337 *out << this->description() << std::endl;
338 break;
342 {
343 *out
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<<"] = ";
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
380template<class Scalar>
382{
383 bool supported = true;
384 for( int i = 0; i < numRowBlocks_; ++i ) {
385 for( int j = 0; j < numColBlocks_; ++j ) {
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
396template<class Scalar>
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;
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
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
435 defaultProductRange_op = ( real_trans(M_trans)==NOTRANS
436 ? defaultProductRange_ : defaultProductDomain_ ),
437 defaultProductDomain_op = ( real_trans(M_trans)==NOTRANS
438 ? defaultProductDomain_ : defaultProductRange_ );
439
441 X = castOrCreateSingleBlockProductMultiVector<Scalar>(
442 defaultProductDomain_op, rcpFromRef(X_in));
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
470template<class Scalar>
471bool
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:
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
535template<class Scalar>
536void
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;
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:
569 }
570
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:
620 }
621}
622
623// Overridden from ScaledLinearOpBase
624
625template<class Scalar>
626bool
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
661template<class Scalar>
662bool
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
696template<class Scalar>
697void
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
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 {
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
749template<class Scalar>
750void
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
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 {
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
805template<class Scalar>
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
821template<class Scalar>
822void 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
834template<class Scalar>
835void DefaultBlockedLinearOp<Scalar>::assertBlockRowCol(
836 const int i, const int j
837 ) const
838{
839#ifdef TEUCHOS_DEBUG
841 !( 0 <= i ), std::logic_error
842 ,"Error, i="<<i<<" is invalid!"
843 );
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()) {
852 !( 0 <= i && i < numRowBlocks_ ), std::logic_error
853 ,"Error, i="<<i<<" does not fall in the range [0,"<<numRowBlocks_-1<<"]!"
854 );
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
867template<class Scalar>
868void 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()) {
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()) {
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
936template<class Scalar>
937template<class LinearOpType>
938void 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
967template<class Scalar>
968void 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
1023template<class Scalar>
1025Thyra::defaultBlockedLinearOp()
1026{
1027 return Teuchos::rcp(new DefaultBlockedLinearOp<Scalar>());
1028}
1029
1030
1031template<class Scalar>
1033Thyra::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
1049template<class Scalar>
1051Thyra::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
1069template<class Scalar>
1071Thyra::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
1089template<class Scalar>
1091Thyra::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
1113template<class Scalar>
1115Thyra::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
1131template<class Scalar>
1133Thyra::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
1151template<class Scalar>
1153Thyra::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
1171template<class Scalar>
1173Thyra::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
virtual std::string description() const
T * get() const
Concrete composite LinearOpBase subclass that creates single linear operator object out of a set of c...
Teuchos::RCP< const VectorSpaceBase< Scalar > > domain() const
bool blockIsConst(const int i, const int j) const
virtual void scaleRightImpl(const VectorBase< Scalar > &col_scaling)
void applyImpl(const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha, const Scalar beta) const
bool acceptsBlock(const int i, const int j) const
Teuchos::RCP< const VectorSpaceBase< Scalar > > range() const
virtual bool rowStatIsSupportedImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat) const
Teuchos::RCP< LinearOpBase< Scalar > > getNonconstBlock(const int i, const int j)
Teuchos::RCP< const LinearOpBase< Scalar > > getBlock(const int i, const int j) const
bool blockExists(const int i, const int j) const
virtual void scaleLeftImpl(const VectorBase< Scalar > &row_scaling)
void setNonconstBlock(const int i, const int j, const Teuchos::RCP< LinearOpBase< Scalar > > &block)
Teuchos::RCP< const ProductVectorSpaceBase< Scalar > > productDomain() const
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Prints the details about the constituent linear operators.
Teuchos::RCP< const ProductVectorSpaceBase< Scalar > > productRange() const
virtual void getRowStatImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat, const Teuchos::Ptr< VectorBase< Scalar > > &rowStatVec) const
void setBlock(const int i, const int j, const Teuchos::RCP< const LinearOpBase< Scalar > > &block)
std::string description() const
Prints just the name DefaultBlockedLinearOp along with the overall dimensions and the number of const...
Teuchos::RCP< const LinearOpBase< Scalar > > clone() const
bool opSupportedImpl(EOpTransp M_trans) const
Returns true only if all constituent operators support M_trans.
Base class for all linear operators.
Interface for a collection of column vectors called a multi-vector.
Base interface for product vectors.
virtual RCP< const VectorBase< Scalar > > getVectorBlock(const int k) const =0
Returns a non-persisting const view of the (zero-based) kth block vector.
Interface for exxtracting row statistics as a VectorBase from a supporting LinearOpBase object.
Applies left or right sclaing to the linear operator.
Abstract interface for finite-dimensional dense vectors.
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
#define TEUCHOS_ASSERT_INEQUALITY(val1, comp, val2)
bool is_null(const std::shared_ptr< T > &p)
#define THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(FUNC_NAME, M, M_T, X, Y)
This is a very useful macro that should be used to validate that the spaces for the multi-vector vers...
#define THYRA_ASSERT_VEC_SPACES_NAMES(FUNC_NAME, VS1, VS1_NAME, VS2, VS2_NAME)
Helper assertion macro.
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.
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
const char * toString(EConj conj)
Return a string name for a EOpTransp value. `*.
EOpTransp real_trans(EOpTransp transp)
Return NOTRANS or TRANS for real scalar valued operators and this also is used for determining struct...
@ TRANS
Use the transposed operator.
@ NOTRANS
Use the non-transposed operator.
@ CONJTRANS
Use the transposed operator with complex-conjugate clements (same as TRANS for real scalar types).
@ CONJ
Use the non-transposed operator with complex-conjugate elements (same as NOTRANS for real scalar type...
TypeTo as(const TypeFrom &t)
T_To & dyn_cast(T_From &from)
std::string toString(const T &t)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)